Whole-genome sequencing of the clinical isolate of Legionella pneumophila ALAW1 from the West Bank allows high-resolution typing and determination of pathogenicity mechanisms

ABSTRACT

Background

Legionella pneumophila is water-based bacterium causing Legionnaires’ disease (LD). We describe the first documented case of nosocomial LD caused by L. pneumophila sequence type (ST) 461 and serogroup 6. The etiology of LD was confirmed by culturing the bronchoalveolar lavage sample retrieving L. pneumophila strain ALAW1. A 7-days treatment of the LD patient with Azithromycin and Levofloxacin allowed complete recovery.

Methods

In details, we sequenced the whole genome of the L. pneumophila ALAW1 using Illumina HiSeq platform. The sequence of ALAW1 was aligned with the genome sequence from the closely related reference strain Alcoy 2300/99 and a whole-genome phylogeny based on single nucleotide polymorphisms (SNPs) was created using Parsnp software. Also, the TYGS web-server was used in order to compare the genome with type strain.

Results

An analysis of the population structure by SNP and TYGS comparison clustered ALAW1 with the reference genome Alcoy 2300/99. Blastp analysis of the type IV secretion Dot/Icm system genes showed that these genes were highly conserved with (≤25%) structural differences at the protein level.

Conclusions

Overall, this study provides insights into detailed genome structure and demonstrated the value of whole-genome sequencing as the ultimate typing tool for Legionella.

KEYWORDS:

• Legionella pneumophila
• ST 461
• VACC11
• L pneumophila str. Alcoy
• WGS

Introduction

L.pneumophila is an opportunistic bacterial pathogen with widespread distribution in freshwater environments. This bacterial species is the main etiology of legionellosis worldwide except Australia and New Zealand [1–3]. The term ‘legionellosis’ describes Legionnaires’ disease (LD), a severe form of atypical pneumonia, and a non-pneumonic febrile illness called Pontiac Fever. L. pneumophila has 15 serogroups (Sgs); Sg1 is the most common causative agent of LD, followed by Sg6 [4,5]. Many studies have demonstrated that an important source for LD is the drinking water distribution systems (DWDS) in large buildings like hospitals and hotels [6–8]. The contamination of hospital water systems with L. pneumophila is considered to pose a high risk for patients, especially for vulnerable and immunocompromised people. To this end, it is well known that LD is an important cause of hospital-acquired pneumonia [9]. The presence of L. pneumophila in DWDS could be a serious health risk to hospital staff and patients, but the magnitude of the problem is often unrecognized [1,9].

For the identification of possible sources of contamination/infection, high-resolution genotyping of new isolates is needed to correlate environmental isolates with clinical isolates. This is currently done using the standard molecular approach Sequence-Based Typing (SBT). SBT of L. pneumophila is done by sequencing a set of seven reference genes per isolate, providing a specific Sequence Type (ST) that can be matched with an International database [10–12]. Another typing method is Multi Locus Variable numbered of tandem repeats Analysis (MLVA) allowing a somewhat higher resolution than SBT for most L. pneumophila strains [13]. Both SBT and MLVA were used for typing environmental and clinical strains of the West Bank in previous studies [7,14,15].

Studies on this ubiquitous water-based pathogen should be directed towards both sides: i) the clinical side from pneumonia patients, and ii) the environmental side from DWDS of hospitals and public buildings [16]. The problem is apparent in the West Bank because awareness about the prevalence of L. pneumophila or LD is lacking, and data scarcity for such arid regions. The study by Sharaby et al [17] described a set of genotyped clinical and environmental L. pneumophila strains retrieved from the closely geographic area Israel. They studied the pathogenicity potential (Ex: hemolytic activity and cytotoxicity toward amoebae and macrophages.) of both clinical and environmental strains. The authors concluded that understanding the virulence characteristics of L. pneumophila genotypes may improve the assessment of public health risks of L. pneumophila in DWDS [17]. Furthermore, in the West Bank, there are no specific guidelines for L. pneumophila surveillance or protection from exposure in hospitals or public buildings [18]. Thus, the main aim of studying entensively the genomics of ALAW1 clinical isolate from this LD case is to focus into genome details, to compare among L. pneumophila intra-clonal genomic characteristics portraying the importance of hospital water habitat.

Case presentation

In January 2014, a 66-year-old woman was admitted to hospital F in Bethlehem – West Bank (Figure S1) (day 0) after evaluation in the emergency ward for high blood pressure and severe gastric pain (Figure 1). The woman had been prescribed Arcoxia (NSAIDs) and Prednisolone (Corticosteroids) 6 days prior to admission (day −6) as treatment of joint inflammation. On day 4 after admission, the patient complained of dyspnea, chest pain, mild fever, and reproductive cough. A chest X-ray showed fluffy infiltrates in both lungs. Atypical pneumonia was suspected. The incubation period of LD is usually 2–10 days, but it has been recorded in some cases to be of up to 16 days. The severity of the disease ranges from a mild cough to a rapidly fatal pneumonia. Initial symptoms include fever, loss of appetite, headache, malaise, and lethargy. Some patients may also experience myalgia, diarrhea, and confusion [19,20]. On day 5, the patient had difficulty breathing with low oxygen and was transferred to the Intensive Care Unit (ICU) of a specialized hospital in East Jerusalem (hospital E) (Figure S1) where bronchoalveolar lavage (BAL) and sputum) were collected to be analyzed by cultivation for bacterial pathogens. On the same day, treatment with Azithromycin followed by Levofloxacin was administered. On day 9, bacteriological culture was negative, and the diagnosis was confirmed as an atypical pneumonia caused by a fastidious pathogen. On day 12, the patient recovered and was discharged (Figure 1).

Figure 1. Timeline of events during investigation of a hospital acquired Legionnaires’ disease — in the West Bank, 2014. Abbreviations: NSAIDs: Non-steroidal anti-inflammatory drugs, ICU: Intensive Care Unit, BAL: Bronchoalveolar Lavage, Sg: Serogroup, ST: Sequence Type and AQU: Al-Quds University. Timeline represents days and the data presented are from one hospitalized patient. BAL sample was analyzed twice, at day 5 and day 10.

On day 17, the clinical isolate ALAW1 was retrieved from the BAL sample, and was identified as L. pneumophila Sg6 Dresden and ST 461 [14]. According to Zayed et al [7] ST 461 was endemic in the DWDS of hospital F in Bethlehem. During the patient’s stay in hospital F, she had two warm showers that may have caused the infection. On day 24 and at the request of the patient, a researcher from the Microbiology Research Laboratory collected and cultured samples from bath shower and faucets in the patients’ home, but the culture was negative for L. pneumophila (Figure 1). The patient’s two stories family house appeared to be in good condition and was connected to a municipal water network rendering a nosocomial infection during her stay in hospital F in Bethlehem most likely.

Materials and methods

L. pneumophila clinical isolate ALAW1 and reference genomes

L. pneumophila clinical isolate ALAW1 was selected for this study. For a comparison, genome sequences of 15 L. pneumophila reference strains retrieved from NCBI GenBank database (https://www.ncbi.nlm.nih.gov/genome/browse#!/prokaryotes/416/) were used as reference (Table S1). Details for all L. pneumophila isolates from the West Bank are in [21].

Short read sequencing, genome assembly, and annotation

Whole-genome sequencing based on short reads was applied to this study. Therefore, the Illumina Hiseq 2500 was used, which is well established and produces low error rates [22,23]. DNA of the clinical isolate was provided for whole-genome sequencing to the Genome Analytics unit (GMAK) of the Helmholtz Center for Infection Research (HZI).

In detail, a DNA library for the Illumina sequencer was prepared using the NEBNext Ultra kit according to the manufacturer’s instructions (NEB, Ipswitch, MA, USA), and sequenced using 100 bp paired end runs on the Illumina Hiseq 2500. Processed sequence reads of the clinical isolate were assembled using Velvet version 1.2.10 with k-mer size of 61.

Assembled contigs of isolate ALAW1 were annotated using Prokka v1.11 [24]. Prokka connects several software tools in order to predict the coordinates of CDS, tRNAs, rRNAs, CRISPRs, and other genomic features encoded on contigs and chromosomes. In this case, Prokka was applied using a genus database (–usegenus) based on the high-quality, manually annotated L. pneumophila strain Corby [25]. With that, 2997 Coding DNA Sequence (CDS) annotations could be successfully transferred to isolate ALAW1.

The Genome Sequence of isolate ALAW1 was uploaded to NCBI GenBank, Acc. No. JALDWK000000000. Furthermore, short read sequencing data is provided in NCBI Sequence Read Archive under BioProject ID PRJNA817377.

Genome analysis

For phylogenomic identification and phylogenomic tree construction, GenBank files of L. pneumophila clinical isolate ALAW1 and 15 reference genomes were submitted to the Type Strain Genome Server (TYGS) (https://tygs.dsmz.de) (Table S1). Hereby, a comparison and clustering of genome data was performed based on the established digital DNA – DNA hybridization (dDDH) threshold of 70% [26]. Subspecies clustering was based on a 93% dDDH threshold as previously introduced [27].

The software Parsnp was used for SNP analysis in order to construct a phylogenomic tree for the core genome of already published L. pneumophila genomes [28]. Parsnp results were visualized using Gingr [28]. Gingr provides an interactive display of multi-alignment variants and phylogenetic trees estimated from the core genome alignment.

Finally, BLASTp was applied to align the amino acid sequences against the Virulence Factors of Pathogenic Bacteria (VFDB) database [29,30]. Amino acid sequences with 75% match identity was chosen and the description of the best hit was assigned as the annotation of predicted gene compared to L. pneumophila str. Philadelphia1 as default bacteria on the webpage.

Results and discussion

Herein, we report, for the first time, that a strain representing L. pneumophila ST461 was isolated from a BAL sample from a patient hospitalized with LD in January 2014.

General features of the genome of L. pneumophila clinical isolate ALAW1

Table 1 summarizes the main genome features of the L. pneumophila clinical isolate ALAW1. The genome of the L. pneumophila clinical isolate ALAW1 was sequenced with an approximately 105× coverage using the Illumina Hiseq platform. 35 contigs greater than 200 bp could be assembled in total summed up to a total genome size of 3,368,497 bp, 3041 genes and 2997 CDS were identified. It contained genomic islands of a total length of 161,380 bp representing 4.8% of the genome.

Table 1. Main features of L. pneumophila clinical isolate used in the study.

Isolate designation	ALAW1_cl_Ps
Sg (mAb)^1	6 Dresden
ST^2	ST461
VACC^3	VACC11
Year of Isolation	2014
Source of isolation	Clinical sample (BAL^)
Genome Size (bp)	3,368,497
No. of Genes	3,040
Genomic Islands (bp)	161,380
Genomic Islands (%)	4.8
CDS*	3,018
GC(%)	38.2
Contigs	40

* Coding DNA Sequence.

¹Serogruop (monoclonal antibody)

²Sequence Type

³´ VNTR-Analysis -Clonal Complex

^Bronchoalveolar Lavage

Identification of L. pneumophila clinical isolate ALAW1 based on SBT and MLVA

The clinical isolate belonged to Sg 6 Dresden and ST 461. By MLVA, ALAW1 is identified as MLVA-genotype (Gt) 9(92). Isolates of ST 461 of the West Bank were split up into three MLVA-genotypes, i.e. Gt 9(92), Gt 10(93) and Gt 10(141). These MLVA genotypes form the Clonal Complex 11 (VACC11), a clonal complex first described for the West Bank [7]. Detailed ST and MLVA patterns and genes are seen in (Table S2).

Comparison of taxonomic resolution of STs and SNP analysis

Identification of L. pneumophila clinical isolate ALAW1 based on Illumina Hiseq Sequencing was performed using the TYGS webserver. Phylogenomic tree of L. pneumophila clinical isolate revealed correct taxonomic task to L. pneumophila. Hereby, the recommended dDDH values of 91.3–94.7% were computed against the type strain L. pneumophila str. Philadelphia1 as default reference strain the web page and thereby fulfill the criteria for bacterial species identification [26]. A phylogenomic tree based on Illumina Hiseq sequencing was constructed using the TYGS web server4 and ALAW1 isolate cluster in the same clade with the reference strains Alcoy 2300/99 and Corby (Figure 2a, for the full TYGS phylogenomic tree see Figure S2). According to the TYGS recommended results of similarity according to the webpage. ALAW1 is most similar to L. pneumophila str. Alcoy 2300/99 with a dDDH value of 93.2% and a G+C content difference of 0.18% (Figure 2a and Table S3)

Figure 2. A) Zoomed in Phylogenomic tree of L.pneumophila clinical isolate (Alaw1_cl_ps) and 15 L. pneumophila reference genomes and subspecies delineation based on the GBDP phylogenetic analyses retrieved from the TYGS website. the dark green squares showed ALAW1 clustered together with the two reference genomes Alcoy 2300/99 and Corby. The branch lengths are scaled in terms of GBDP distance formula d4 (the whole phylogenomic tree of Legionella species and subspecies in supplementary materials (Figure S2). B) Gingr visualization of 15 L. pneumophila reference genomes and L.pneumophila clinical isolate aligned with Parsnp. The leaves of the reconstructed phylogenetic tree (left) are paired with their corresponding rows in the multi-alignment. L. pneumophila str. Alcoy 2300/99 is the reference genome. The constructed tree share 72% of the core genome. The gray region isn’t a part of core genome. The white region is a part of core genome but no SNPs accumulation. The heatmap on the button of the figure indicates more reddish meaning more SNPs accumulation..

Additionally, a phylogenomic tree was constructed using SNP analysis of 15 reference genomes from GeneBank (Table S1) and our clinical L. pneumophila genome sequenced strain. The main branches or clusters were observed in the constructed phylogenomic tree (Figure 2b). ALAW1 is clustered with the two complete reference genomes (L. pneumophila strain Alcoy 2300/99 and Corby) (Figure 2b).

By ParSNP analysis L. pneumophila str. Alcoy 2300/99 was shown to be the closest reference strain so far with a complete genome available in the GenBank database. We found 11,674 SNPs compared to L. pneumophila clinical isolate. Indeed, L. pneumophila str. Alcoy 2300/99 showed 1,303 SNPs less than performing the complete ParSNP analysis with L. pneumophila str. Corby as a reference genome. Thus, the best result for our SNP-based phylogenomic tree was retrieved with L. pneumophila str. Alcoy 2300/99 as reference and documented in Figure 2b and Table 2. Overall, we observed a typical concordance between phylogenomic tree based on SNPs and phylogenomic tree based on TYGS webpage.

Table 2. Whole-genome SNP comparison of L. pneumophila isolates and reference strains.

L. pneumophila Clinical Isolate	ALAW1_cl_Ps
Lpn-LPE509	55,387
Lpn-Phildelphia1	58,320
Lpn-Thunderbay	56,257
Lpn-ATCC43290	56,577
Lpn-lpm7613	56,771
Lpn-Lens	69,093
Lpn-Lorraine	52,304
Lpn-HL06041035	47,447
Lpn-D7630	46,594
Lpn-Paris	45,200
Lpn-OLDA	45,262
Lpn-Pontiac	47,414
Lpn-Toronto	47,230
Lpn-Alcoy 2300/99	11,674
Lpn-Corby	12,977

Bold = Least No. of SNPs.

Results of TYGS phylogenomic tree of L. pneumophila clinical isolate are identical with the constructed SNP phylogenomic tree. Comparing the phylogenomic tree obtained with these two molecular methods and mathematical algorithms, we found that ST is in concordance with the SNP genotyping, but it cannot cover the needed high resolution at intra-clonal level provided by WGS.

Actually, a threshold for the number of SNPs necessary to identify outbreak-associated L. pneumophila strains still has to be established. For example, more than 200 SNPs have been reported in phylogenetically closely related L. pneumophila outbreak strains [31,32]. Mercante et al [33] showed that up to 20 core SNPs were identified in comparison of Philadelphia clade L. pneumophila isolates. In the present case study, 11,674 SNPs were identified in the same sub-species cluster with the reference strain Alcoy 2300/99 from Valencia, Spain. While SBT from the same sub-lineage in a phylogenomic cluster is different as L. pneumophila str. Alcoy 2300/99 is ST578 and ALAW1 is ST461. These results are in concordance with Khodr et al [34]. They sequenced six ST1 genomes (four clinical and environmental isolates from a hospital and the other two were unrelated) and observed that geographically unrelated isolates differed by more than 1,500 SNPs (Table 2).

Strains of VACC11 are affiliated with ST461 and were classified as Sg 6 Dresden [7]. ST461 was previously reported by the European Working Group for Legionella infections (EWGLI) [7,8,35,36], to be found in hospitals water in Poland [35]. Moreover, ST461 was identified in Michigan, (USA) water systems and showed high capability to efficiently infect THP-1 macrophages [36]. More recently, ST461 was identified in hotel water in southern Israel [8] and in the West Bank hospital water systems [7].

Identification of pore-forming genes mediating cytotoxicity in L. pneumophila clinical isolate ALAW1

Central to the pore-forming mediated cytotoxicity of L. pneumophila are the Dot/Icm loci, which taken together directly assemble to a type IV secretion system (T4SS) [37,38]. Also, the toxin rtxA plays an important role in the pore-mediated cytotoxicity [39–41]. Although all L. pneumophila strains examined until today contain the complete Dot/Icm loci, sequence variations among the Dot/Icm genes among different L. pneumophila strains have been reported [37,42]. Eleven Dot/Icm T4SS genes (icmT [43,44], icmS [45], icmR [45], icmQ [45], icmL/dotI [46], icmK/dotH [46], icmE/dotG [46], icmC/dotE [46], dotB [46], dotA [46] and icmW [45]) and rtxA gene [41] are responsible for the pore-forming mediated cytotoxicity of L. pneumophila. The 11 Dot/Icm T4SS pore-forming mediated cytotoxicity genes and the rtxA gene were identified by a BLASTp search against the Virulence Factor Database (VFDB) using L. pneumophila strain Philadelphia1 as default reference genome for the reference strain L. pneumophila str. Alcoy 2300/99 and L. pneumophila clinical isolate ALAW1, which had been genome sequenced (Table 3).

Table 3. Percentage of nucleotide identity of orthologous pore-forming activity genes with respect to the BLASTp search against the VFDB using L. pneumophila strain Philadelphia1 as default reference genome.

	Secretion system
	Dot/Icm type IVB secretion system
L. pneumophila strain	icmT		icmS		icmR		icmQ
	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)
Lpn-Alcoy 2300/99	82%	260	100%	344	95%	362	100%	575
ALAW1_cl_Ps	82%	260	100%	344	95%	362	100%	575
L. pneumophila strain	icmL/dotI		icmK/dotH		icmE/dotG		icmC/dotE
	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)
Lpn-Alcoy 2300/99	84%	638	84%	1,085	88%	3,146	99%	584
ALAW1_cl_Ps	84%	638	83%	1,085	88%	3,146	99%	548
							Toxin
							RtxA
L. pneumophila strain	dotB		dotA		icmW		rtxA
	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)	Id	Length (bp)
Lpn-Alcoy 2300/99	100%	1,133	81%	3,059	98%	455	85%	14,009
ALAW1_cl_Ps	100%	1,133	78%	3,119	98%	455	75%	4,877

Table 3 shows 75% to 100% similarity between the pore-forming mediated cytotoxicity genes. Morozova et al [47] showed that the Dot/Icm genes are highly conserved in L. pneumophila strains. After WGS technology was available on the market, Gomez-Valero et al [48,49] confirmed the previous study showing high conservation (98%) among orthologs of the reference strains Corby, Paris, Philadelphia, and Lens with few exceptions in the dotA gene. The dotA gene is an essential gene for virulence activity of L. pneumophila strains since it encodes an integral membrane protein with eight domains. This explains why a dotA mutant of L. pneumophila strain Corby is being used as a negative control for all virulence assays [50]. Costa et al [51] analyzed 300 dotA gene sequences from L. pneumophila strains and demonstrated that pathogenic L. pneumophila strains belong to a subset of the genotypes existing in the environment. Khodr et al [34] explained the high variation of the dotA gene of L. pneumophila by indicating that this gene is a target for host speciation and adaptive evolution to different hosts and environments. Dumenil et al [52] showed that IcmR is a regulator gene for the IcmQ gene that possesses pore-forming activity. In addition, Gomez-Valero et al [53] demonstrated that dotB, IcmS and IcmW are highly conserved genes. These facts are in accordance with our results (Table 3) showing that IcmR (95%), IcmS and dotB (100%) and IcmQ (100%) were highly conserved gene while dotA had only a 78% gene similarity for our L. pneumophila clinical isolate ALAW1. More recently, Zayed et al studied the infectivity and cytotoxicity to all clonal complexes and their affiliated genotypes of the West Bank [54]. Overall, the Dot/Icm system is a highly conserved and complex molecular system (Table 3).

Taken together, this detailed analysis suggests that the 12 genes studied out of more than 200 of the Dot/Icm components of L. pneumophila are representing a large repertoire of effectors, which are necessary for virulence [48]. In general, both of the L. pneumophila strains studied shared the same Dot/Icm T4SS with less than 25% structural differences at the protein level.

According to David et al [55] a certain local micro-evolution could be observed if isolates from the same site of isolation have been obtained at different times. Such a case of micro-evolution could have occurred in the environmental isolates obtained from the DWDS of hospital F. Environmental strains of ST461 has had time to diversify by genetic drift since it was endemic in Bethlehem area (hospital F) and the clinical isolate originated from environmental source. Also, conditions in the water systems such as disinfection, temperature, and amoeba might have increased the virulence of the ST461 strains. A thorough genome-based comparison of the environmental strains and the clinical isolate ALAW1 might give insight into this development and the respective mechanisms.

Conclusions and future plans

This study of L. pneumophila clinical isolate ALAW1 highlights conserved genetical features that may be critical for pathogenesis in human lungs. Future plans aim to compare the clinical isolate genome with genome sequenced environmental isolates belonging to ST461 from hospital F hotspots. Moreover, complete genome sequencing of L. pneumophila ST461 strains is highly recommended to better understand the evolution and environmental adaptation within the VACC11 clade.

Patient consent

This study has consent form PHRC/HC/211/17

Acknowledgements

We thank Josefin Koch, Verena Maiberg and Suha Butmeh for their skilful technical support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Data are available at GenBank accession numbers JALDWK000000000. Furthermore, short read sequencing data is provided in NCBI Sequence Read Archive under BioProject ID PRJNA817377.Full WGS methods and results are provided in the Supplementary Data.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/20018525.2023.2168346.

Additional information

Funding

This work was supported by funds from the Deutsche Forschungsgemeinschaft (DFG-project Number HO-930/5-1/2) and the Hildegard-Elisabeth Foundation, Germany