Molecular evolution of virulence genes and non-virulence genes in clinical, natural and artificial environmental Legionella pneumophila isolates

L. pneumophila isolates

One hundred and thirty-nine strains of L. pneumophila were enrolled in this study. These isolates included 51 artificial isolates, 59 natural isolates, and 29 clinical isolates. The source natures, geographic locations, collection dates and sequence types (STs) based on the Sequence-Based Typing (SBT) scheme (Gaia et al., 2005; Ratzow et al., 2007) of these isolates are shown in Table S1. Briefly, the 29 clinical strains were isolated between 1947 and 2012, from the non-China regions, including different cities of USA, Germany, UK, France, etc. Their details were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov). The environmental isolates were from 14 different sites in two cities (Guangzhou and Jiangmen) of Guangdong Province, China between October 2003 and September 2007. All the environmental isolates were selected for sequencing partial trpA, cca , lssD, lspE, and icmK genes. We selected the most variable regions of these genes through a sequence alignment with the known sequences in the NCBI database, in order to achieve maximum genetic variability. The variability of the gene regions was measured by calculating the single-nucleotide polymorphisms of the sequence using Dnasp v5 (Librado & Rozas, 2009; Rozas, 2009) (Table S2).

Genomic DNA extraction, PCR, and DNA sequencing

Genomic DNA extraction of the artificial and natural strains was performed as shown in our previous report (Zhan, Hu & Zhu, 2016). PCR was employed to amplify fragments of DNA. The corresponding oligonucleotide primers are shown in Table S2. The PCR was performed using an EasyPfu PCR SuperMix (Transgene Biotech, Beijing, China) according to the manufacturer’s instructions and carried out using the GeneAmp PCR system (MJ Research PTC-200) with the following thermal conditions: 95 °C for 3 min followed by 35 cycles of 95 °C for 20 s, 60 °C for 20 s and 72 °C for 30 s (lspE, lssD, and icmK loci) or 70 s (cca and trpA loci), and a final extension at 72 °C for 5 min. For confirmation purposes, each PCR reaction was performed with a positive control (L. pneumophila strain ATCC33152 genomic DNA as the PCR template) and a negative control (sterile water as the PCR template). PCR products were purified using an EasyPure Quick Gel Extraction (Transgene Biotech, Beijing, China) and then transferred to Guangzhou IGE Biotechnology Ltd for sequencing.

Sequence analysis

The quality of DNA sequencing was manually checked by Chromas (http://technelysium.com.au). The corresponding gene sequences of 29 clinical isolates were obtained from the NCBI database. Multiple sequence alignments were performed using ClustalX 2.1 (Chenna et al., 2003; Thompson et al., 1997). Genetic variability analyses were performed using DnaSP v5 (Librado & Rozas, 2009; Rozas, 2009). Phylogenetic analyses were conducted by a MEGA7 package (Kumar, Stecher & Tamura, 2016). Neighbor-Joining (NJ) phylogenetic trees were obtained for each locus separately with the MEGA7 based on the Kimura 2-parameter model (Kimura, 1980; Saitou & Nei, 1987). The tree was drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. NJ tree nodes were evaluated by bootstrapping with 1,000 replicates. The ratios of synonymous and non-synonymous substitutions were calculated according to the Nei-Gojobori method with Jukes-Cantor correction as implemented in the MEGA7 (Kumar, Stecher & Tamura, 2016).

Molecular evolution analysis

The aligned sequences of the five loci were screened using RDP4 to detect intragenic recombination (Martin et al., 2015; Martin et al., 2017). Six methods implemented in the program RDP4 were utilized. These methods were RDP (Martin & Rybicki, 2000), GENECONV (Padidam, Sawyer & Fauquet, 1999), BootScan (Martin et al., 2005), MaxChi (Smith, 1992), Chimaera (Posada, 2002), and SiScan (Gibbs, Armstrong & Gibbs, 2000). Potential recombination event was considered as that identified by at least two methods according to Coscolla’s report (Coscolla & Gonzalez-Candelas, 2009). Common settings for all methods were to consider sequences as linear, statistical significance was set at the P < 0.05 level, with Bonferroni correction for multiple comparisons and requiring phylogenetic evidence and polishing of breakpoints.

Tajima’s D, Fu and Li’s D* and F* were calculated for testing the mutation neutrality hypothesis as previously described by Coscolla and colleagues (2007). These statistics were calculated with the program Dnasp v5 (Rozas, 2009) using a statistical significance level P < 0.05 and applying the false discovery rate to correct for multiple comparisons, and 1000 replicates in a coalescent simulation. Non-neutrality evolution was considered as that identified by at least two of the methods.

Population structure analysis

Hierarchical analysis of molecular variance (AMOVA) for clinical and environmental sequences using the alignment of the five loci was performed using Arlequin Ver3.5.2 (Excoffier & Lischer, 2010). This analysis provides estimates of variance components and F-statistics analogues speculating the correlation of haplotype diversity at different levels of the hierarchical subdivision. We defined the hierarchical subdivision of these isolates at three levels. At the upper level, the three groups considered were clinical, artificial and natural isolates. As populations within groups, the intermediate level, we reckoned the strains isolated from the same geographic location as subpopulations. Therefore, natural and artificial isolates were both split into two subgroups based on the cities where they were isolated (Guangzhou and Jiangmen subgroups), and clinical strains were also divided into two subgroups based on the continents where they were isolated (eg. America or non-America, including Europe and Australia). The third level corresponded to the different haplotypes which were found within the six subgroups considered in the previous level. The statistical significance of fixation indices was tested using a non-parametric permutation approach (Excoffier, Smouse & Quattro, 1992).

Nucleotide sequence accession numbers

The 550 sequences from L. pneumophila environmental isolates determined in this study were deposited in the GenBank Nucleotide Sequence Database with accession numbers KY708328 –KY708437 (cca), KY708438 –KY708547 (trpA), KY708768 –KY708877 (lssD), KY708658 –KY708767 (lspE) and KY708548 –KY708657 (icmK).

Sequence analysis and genetic variability of L. pneumophila isolates from different sources

In general, we obtained the gene sequences from 29 clinical, 59 natural and 51 artificial environmental isolates. Genetic diversity estimates in all the clinical, natural and artificial isolates are presented in Table 1. The highest nucleotide diversity (π) in the L .pneumophila isolates was found in the lssD locus, varied from 0.04213 to 0.05399, while the lowest nucleotide diversity was found in trpA locus, varied from 0.01041 to 0.01246. Both the most haplotype (h) and highest haplotype diversity (Hd) of the five gene loci was found in natural isolates. For the trpA locus, the nucleotide diversity, number of polymorphic nucleotide sites (S), population mutation ration (θ), average number of pairwise nucleotide differences (k), and the total number of mutations (η) were all higher in clinical isolates. In contrast, another NV locus, cca did not show higher nucleotide diversity and nucleotide differences in clinical isolates, but the number of polymorphic nucleotide sites, population mutation ratio and the total number of mutations were higher in clinical isolates. The highest nucleotide diversity and the average number of pairwise nucleotide differences of cca locus were found in artificial isolates. Different results were found in the three virulence loci: the haplotype, haplotype diversity, nucleotide diversity and the average number of pairwise nucleotide differences were higher in natural isolates.

The rates of non-synonymous substitutions per non-synonymous site (dN) were extremely low and different between these genes, ranging from 0.00243 (for lspE in clinical isolates) to 0.0295 (for icmK in natural isolates), despite the relatively similar values of polymorphic sites (37 vs. 77). Synonymous substitutions (dS) ranged from 0.0338 for trpA in natural isolates to 0.3091 for lssD in natural isolates (Table S3). An obviously different dN/dS ratio was observed, ranging from 0.0165 (for lssD in clinical isolates) to 0.2543 (for icmK in natural isolates) (Table 1). We found different dN/dS ratios among clinical, natural, and artificial isolates. dN/dS ratios of the cca, trpA, lspE, and icmK locus were higher in natural isolates than in artificial and clinical isolates (Table 1).

Phylogenetic analysis of the five gene loci

Phylogenetic trees were derived separately for each locus to test the phylogenetic relationships between these isolates. The NJ trees showed three to six main groups (Figs. 1 to 5). The number of haplotypes in these loci ranged from 15 (lssD locus) to 89 (icmK locus). The NV loci displayed a similar number of haplotypes (23 for cca and 22 for trpA), while the virulence loci displayed large variable range of haplotypes (16 for lssD, 19 for lspE and 89 for icmK). More than half of the clinical isolates (16/29 to 18/29, 55.17% to 62.07%) presented a single dominant allelic profile in the five loci (Figs. 1 to 5). Similarly, more than half (27/51 to 38/51, 52.94% to 74.51%) of the artificial isolates presented a distinct dominant allelic profile of the four loci (Figs. 1 to 4). Some of the allelic profiles were unique to clinical, artificial or natural isolates. Six allelic profiles of cca were unique to clinical isolates, comprising 26.09% of all profiles. Two to six allelic profiles of trpA, lssD, lspE and icmK loci were unique to clinical isolates, constituting 6.74% to 25% of the select profiles. A significantly different proportion of unique alleles of the five loci between natural and artificial strains was found (P = 0.0008, paired t-test). 30.43% (7/23), 40.91% (9/22), 40% (6/15), 42.11% (8/19), and 50.56% (45/89) allelic profiles of the cca, trpA, lssD, lspE and icmK loci were unique to natural isolates, respectively. By contrast, only 4.35% (1/23), 13.64% (3/22), 6.67% (1/15), 15.79% (3/19), and 35.96% (32/89) profiles of the cca, trpA, lssD, lspE and icmK loci were unique to artificial isolates. Most of the clinical isolates distributed in group A cluster of cca locus (86.21%, 25/29) and mostly in the subgroup 2 (cca-A-S2, 68.97%, 20/29); group B cluster of trpA locus (68.97%, 20/29) and mostly in the subgroup 2 (trpA-B-S2, 55.17%, 16/29); group B cluster of lssD locus (58.62%, 17/29), group A cluster of lspE locus (79.31%, 23/29) and mostly in the subgroup 2 (lspE-A-S2, 62.07%, 18/29); group B cluster of icmK locus (58.62%, 17/29). Artificial environmental isolates mainly distributed in group B cluster of cca locus (64.71%, 33/51), group A cluster of trpA locus (70.59%, 36/51), group A clusters of lssD locus (76.47%, 39/51), group A clusters of lspE (82.35%, 42/51) and group A cluster of icmK (84.31%, 43/51) (Figs. 1 to 5). The natural isolates in the NJ trees of cca, trpA, lssD, and lspE were more dispersed, while in the NJ tree of icmK, there mainly distributed in the group A, subgroup 3 and 4 clusters (Fig. 5).

Relationships between clinical, natural and artificial isolates

We performed a hierarchical analysis of molecular variance (AMOVA) for the 139 isolates based on five loci considered above. The largest proportion of the genetic variation in each locus was found within populations, as this level accounted for 76.39% to 90.88% of the total variation (Table 2). The proportion of the total genetic variation explained by differences between clinical, artificial and natural isolates was small, ranging from −5.04% for lssD to 18.87% for cca, and it was not significant (Table 2). In contrast, geographical difference contributed to a part of genetic variation, especially for the virulence loci (accounting for 10.38% to 21.78% variation), and was all significant.

Evolution and recombination analysis

Tajima’s D, Fu and Li’s D* and F* statistics were calculated for testing the mutation neutrality hypothesis. The results showed most of the genes were in accord with the neutral hypothesis, except the NV genes of clinical isolates (Fig. 6, Table 3).

The intragenic recombination events of each gene locus in different types of isolates were detected using RDP4 individually. For clinical isolates, RDP reported recombination events happened in the cca locus of C17, C21, C22, and C23, in the lssD locus of C17, and in the lspE locus of C16. No recombination event was identified on these loci of artificial environmental isolates. However, recombination events were detected only in virulence loci of natural isolates including the lssD locus of N71 and N220, the lspE locus of N52 and N53 and the icmK locus of N115 and N152. Details are shown in Table 4. These results might indicate that horizontal exchange of genetic material of virulence genes in clinical and natural isolates was more widespread than that in artificial isolates, and it was more prevalent in natural isolates.