Virulence effector SidJ evolution in Legionella pneumophila is driven by positive selection and intragenic recombination

Introduction

Legionella pneumophila (L. pneumophila) is the most common causative agent of legionellosis, which manifests as atypical pneumonia or non-pneumonia type illnesses, e.g., Pontiac fever (Fields, Benson & Besser, 2002). As a pathogenic Gram-negative bacterium, L. pneumophila is widely present in natural environments including natural water and soil sources, in which free-living amoebae is its natural host (Albert-Weissenberger, Cazalet & Buchrieser, 2007; Van Heijnsbergen et al., 2014). From natural environments, L. pneumophila can colonize artificial environments (e.g., cooling towers and hot-water systems), and then spread in aerosols, infecting the susceptible person. As an intracellular parasite for protozoa, L. pneumophila infects mammalian host cells using similar mechanisms: phagocytosis or macropinocytosis (De Carvalho, Barrias & De Souza, 2015; Horwitz, 1984; Peracino, Balest & Bozzaro, 2010). After entering the host cells, L. pneumophila creates an intracellular niche named Legionella-containing vacuoles (LCVs), which are resistant to the acidification and LCV-lysosome fusion, and permissive for its replication (Isberg, OConnor & Heidtman, 2009). In this process, host-cell functions are modulated by hundreds of effector proteins encoded by the L. pneumophila genome and delivered by L. pneumophila type IV Dot/Icm secretion system (Hubber & Roy, 2010; Segal, Feldman & Zusman, 2005). The genomes of L. pneumophila consist of a single circular chromosome of about 3.4 Mb (Gomez-Valero et al., 2011). Some L. pneumophila strains (e.g., Paris, Lens and Lorraine) also contain a plasmid (Gomez-Valero et al., 2011). The number of genes in L. pneumophila chromosome is about 3,000, of which, 98%–99% are protein-coding genes and about 300 are the type IV Dot/Icm effectors (Gomez-Valero et al., 2011). As a facultatively pathogenic bacterium interacting with free-living amoebae, L. pneumophila exhibit a genome larger than their close relatives such as Coxiella burnetii and Francisella tularensis due to gene conservation and acquisition (Moliner, Fournier & Raoult, 2010). This was also proved by Gomez-Valero et al. (2019) that the number of gene gain events in 2,837 representative proteins of genus Legionella considerably exceeded the number of loss events.

Functional redundancy among groups of effector proteins is required for L. pneumophila to survive in different host cells (Newton et al., 2010; Richards et al., 2013; Shames et al., 2017). However, only a few of these proteins are necessary for intracellular replication, and elimination of numerous effector genes rarely leads to detectable defects in intracellular growth (O’Connor et al., 2011). There were some critical components for both intracellular growth and disease within animals that have been identified in L. pneumophila, including SdhA, SidJ, and AnkB (Al-Khodor et al., 2008; Anand, Choi & Isberg, 2020; Harding et al., 2013; Jeong, Sexton & Vogel, 2015; Liu & Luo, 2007). SidJ, encoded by a 2,622 to 2,628 bp length gene, is a member of the Dot/Icm effector and plays a key role in regulating several host cellular processes and pathways through another effector member named SidE family (SidEs), including SdeA (1,499 aa), SdeB (1,920 aa), SdeC (1,533 aa) and SidE (1,495 aa) (Liu & Luo, 2007). These proteins are encoded by genes with lengths of 4,497 bp, 5,760 bp, 4,599 bp, and 4,485 bp respectively. SidEs localize to the cytoplasmic face of the LCV in the early stages of L. pneumophila infection. They are required for the mono-ADP-ribosyltransferase activity involved in ubiquitin activation, which is regulated by SidJ glutamylase activity. Such activity is subsequently modulated by the eukaryote-specific protein calmodulin (CaM) via binding (Gan et al., 2019; Liu & Luo, 2007). SidJ catalyzes glutamylation of SidEs (Bhogaraju et al., 2019; Black et al., 2019), and in turn, inhibits their unrestrained ubiquitin (Ub) ligase activity; which is shown to be harmful to the host by poisoning the cellular Ub pool and possibly blocking the action of other L. pneumophila effectors that manipulate the host Ub machinery (Bhogaraju et al., 2016). Three functional domains of SidJ were identified in a previous report, including a kinase domain (KD) in the center of protein spanning residues 336 to 593, which forms a catalytic structure and an N-terminal (NTD) and C-terminal domains (CTD) (Black et al., 2019). Given the crucial role of SidJ in this regulation network, studying the nature of evolution of sidJ is of great importance to understand the virulence in L. pneumophila and the interaction between the bacteria and its hosts. A study on 32 unrelated strains of L. pneumophila revealed that recombination was an important strategy in the evolutionary adaptive process and played an active role in sidJ genetic plasticity (Costa et al., 2014). Recombination of sidJ might also provide a broad-host-range for L. pneumophila by preventing host specialization and contributing to the resilience of the species (Costa et al., 2014). Besides recombination, selection was another fundamental evolutionary force that shaped DNA sequence variation. Interaction between recombination and natural selection within a gene can either increase or decrease sequence diversity. Moreover, recombination can generate genetic variation, which is tested by natural selection, and as such, it plays an important role in fueling adaptive evolution (Jouet, McMullan & Van Oosterhout, 2015).

The ultimate goal of this work is to understand the underlying patterns in the evolution of the sidJ gene of L. pneumophila during its lifecycle (e.g., infect the amoebae via amoebal attack and present a sympatric lifestyle) through identifying individual codons under positive selection. We utilized various algorithms to identify the molecular evolution patterns of sidJ in a relatively large number of L. pneumophila strains. It is shown here both intragenic recombination and positive selection drove the adaptive evolution of sidJ and shaped its genetic plasticity. Codons of sidJ that were identified to experience positive selection might play key roles in regulating the binding affinity of SidJ to CaM; and thus, change the glutamylases activation of SidJ, which might, in turn, manipulate the host Ub machinery balance. This reticular regulation network might be an important strategy for the survival and adaptability of L. pneumophila to variable host cells.

Materials and Methods

Mapping of positively selected sites to structure models of proteins

The three-dimensional structure of SidJ and CaM was modeled using the Phyre server (Kelley et al., 2015), and the SWISS-MODEL (http://swissmodel.expasy.org) (Waterhouse et al., 2018). The positive selection sites were mapped onto the structure and visualized by PyMol (http://www.pymol.org/) (Lilkova et al., 2015).

Results and Discussion

Intragenic recombination drives the diversification of the sidJ

To explore whether intragenic recombination drives the diversification of the sidJ, we categorized the 116 strains into two groups: the recombinant group and the non-recombinant group. Then, we utilized DnaSP to study the difference in genetic diversity between the two groups. It showed that most of the parameters that represent genetic diversity were higher in the recombinant group. These parameters included nucleotide diversity, polymorphic sites, and the average number of nucleotide differences (Table 2), suggesting that recombination added a high density of polymorphisms in sidJ. The phylogeny of sidJ alleles also showed that recombinant alleles roughly formed an outside subclade (sidJ-3 and sidJ-4, Fig. 1B) compared with those non-recombinants. Recombination also introduced an excess of non-synonymous and synonymous diversity for sidJ (0.02454 Vs. 0.01371 and 0.1660 Vs.0.1063, Table 2). These results were partly consistent with the research with another intracellular bacteria, Mycobacterium tuberculosis, of which the genetic diversification was partly driven by recombination and led to a high genetic diversity and genomic plasticity (Namouchi et al., 2012). Furthermore, we did not find evidence for recombinant sidJ alleles experiencing demographic expansion based on the current pools of strains. The mismatch distributions for the total sidJ set, recombinant sidJ, and non-recombinant sidJ genes were roughly multimodal with P > 0.05 for the SSD. Also, P-values for Harpending’s Raggedness index was lower than 0.05 in each group, indicating that no demographic expansion exists (Figs. 2A–2C). However, a potential reduction of recombinant sidJ alleles in the L. pneumophila community was found. Fu and Li’s D* and F* tests showed significantly positive values (Table 2), indicating an excess of intermediate-frequency alleles which might result from bottleneck populations, thus causing demographic reduction (Rowbotham, 1980). Parsimony (TCS) network of sidJ alleles showed no central allele (node). Many mutations were found among the alleles, and these alleles did not comprise a scattered star structure (Fig. 3), suggesting that the expansion of the L. pneumophila population with a specific mutation in the sidJ gene has not taken place (Leigh & Bryant, 2015).

Table 2:

Summary of genetic diversity parameters for sidJfromL. pneumophila strains.

Parameters	Over all	Non-recombinant alleles	Recombinant alleles
Sequences, n	116	83	33
Haplotypes, h	39	25	14
Haplotype diversity, Hd	0.909	0.843	0.866
Nucleotide diversity, π	0.05317	0.03178	0.05123
(standard deviation)	0.00278	0.00342	0.00193
Polymorphic sites, S	544	390	361
Theta per site (from S)	0.03900	0.03178	0.03396
(standard deviation)	0.00925	0.00755	0.01029
Average number of nucleotide differences, k	139.246	83.255	134.163
Total number of mutations, Eta	600	414	381
dN	0.02344	0.01371	0.02454
dS	0.1833	0.1063	0.1660
dN/dS	0.1278	0.1290	0.1478
Tajima’s D	0.7917 (P > 0.10)	0.01120 (P > 0.10)	1.64769 (P > 0.10)
Fu and Li’s D*	1.46609 (0.10 > P > 0.05)	1.30041( P > 0.10)	1.62736 (P < 0.02)
Fu and Li’s F*	1.38900 (P > 0.10)	0.90431 (P > 0.10)	1.94078 (P < 0.02)

DOI: 10.7717/peerj.12000/table-2

Table 3:

Log-likelihood values and parameter estimates for the sidJgene ofL. pneumophilausing modifiedtopology tree of the alleles.

Model	nP	lnL	Estimates of parameters	LRT P-value	Positively sites
M3 (discrete)	81	−9330.3176	p0 = 0.8126, p1 = 0.1727, p2 = 0.01473, ω0 = 0.01711, ω1 = 0.7739, ω2 = 2.5605	P <10−9	200N*, 868T*, 869S**
M0 (one ratio)	77	−9567.0414	ω0 = 0.1636		Not Allowed
M2a (selection)	80	−9331.9663	p0 = 0.8393, p1 = 0.1533, p2 = 0.00744, ω0= 0.02473, ω1 = 1.00000, ω2= 3.1041	P=0.0803	N/A
M1a (neutral)	78	−9334.4879	p0 = 0.8392, p1 = 0.1608 ω0= 0.02390, ω1= 1.0000		Not Allowed
M8a (beta& ω)	80	−9333.3117	p0 = 0.8654, p = 0.04277, q = 0.3594 p1 = 0.01345,ω = 1.1155	P <10−9	58G**, 200N**, 820A*, 867R*, 868T**, 869S**
M7 (beta)	78	−9358.5621	p = 0.03809, q = 0.17020		Not Allowed

DOI: 10.7717/peerj.12000/table-3

Notes:

P is the number of parameters in the ω distribution; lnL is the log likelihood; ω is ratio of dN/dS, LRT P-value indicates the value of chi-square test; Parameters indicating positive selection are presented in bold; positive selection sites were identified by the Bayes empirical Bayes (BEB) methods under M8 model or by Naive Empirical Bayes (NEB) methods under M3 and M2a models.

The posterior probabilities (p) ≥0.90, (p) ≥0.95 and p ≥ 0.99 are indicated by *, ** and ***, respectively. (Yang, 2007; Yang & Bielawski, 2000) recommended that results from M8 model were preferred to find sites under positive selection pressure, and it is more robust to recombination which was proved by Anisimova, Nielsen & Yang (2003)

Recombination and positive selection shape the population structure of SidJ in L. pneumophila

The evolutionary history of the SidJ proteins corresponding to 39 alleles was studied by using MEGA X. Similar topology of the trees was found when compared to that of the sidJ alleles (Figs. 1B and 4). Three paired alleles (Lorraine and ATCC35096, Ice27 and Aco13, and D5265 and D4954) encoded SidJ with the same protein sequences. The properties of these alleles and the information of their representative strains were studied. Of the 39 sidJ alleles, some were distributed both in clinical and environmental strains, while some were only distributed in environmental strains (Table S1). We defined those alleles that could be found in clinical strains as clinically associated alleles. Among the 14 recombinant alleles, 12 (85.71%) were clinically associated. In contrast, among the 25 non-recombinant alleles, only 15 (60%) were clinically associated. This result suggested that recombinant sidJ alleles were more likely to be clinically associated alleles, although not significant (P =0.093, Chi-Square test). (Fig. 5A). A larger pool of L. pneumophila strains is required to sufficiently explain the association of recombinant sidJ with clinical strain. Based on this result, we propose that recombination is an important strategy for L. pneumophila to survive in different environments, and for infecting human cells. We did find specific mutation profiles of positive selection sites in clinically associated alleles (Fig. 5B) or recombinant alleles (Fig. 5C), for example G58M mutation happened less frequently in clinically associated alleles, while T868P, T869P happen more frequently in recombinant alleles of sidJ (Table S5 and Figs. 5B–5C). Considering that the positive selection is usually beneficial to the survival of the individual bacteria carrying the mutation, these results indicated that mutations in positive selection sites increased SidJ variability, and made more extensive the adaptability in environmental hosts for L. pneumophila. This might lead to broad coevolution of L. pneumophila genes (e.g., sidJ) with viable environmental hosts before it could infect human cells and thus be of crucial importance in the virulence of this bacteria. The fact that the finding of alleles harboring particular mutations in positive selection sites were more likely recombinants further demonstrated that recombination in sidJ enhanced the environmental adaptability of some L. pneumophila strains (Fig. 5C).

Mutation of positive selection sites of SidJ might influence the binding of CaM to SidJ

Here, the four definitive positive selection sites are located in the NTD and CTD, but not in the KD. CaM binding is required by SidJ glutamylase activity. An IQ (I841 and Q842) motif located in the C-terminal domain of SidJ is involved in CaM binding by burying in a hydrophobic cleft of the CaM C lobe. The CTD of CaM semi-encircles the C-terminal helix of SidJ and the NTD domain of CaM makes extensive contacts with the N-terminus of SidJ (Bhogaraju et al., 2019). Residues in these domains also play roles in mediating the formation of the SidJ-CaM complex, and in turn, stabilize the position of the N-lobe of the KD, and thereby leading to the formation of a stable catalytic pocket for SidES (Bhogaraju et al., 2019). A previous site-directed mutagenesis study verified the importance of I841 and Q842 for optimal binding. In addition, some other sites including Q830, S808, E812, R796, R660, R804 that engaged in hydrogen-bonding interactions with corresponding CaM residues also showed the importance of the binding affinity of SidJ with CaM (Gan et al., 2019). Therefore, the mutation on the positive selection sites might change the level of interactions of SidJ with CaM through hydrogen bonds and salt bridges. Structural studies on a truncated SidJ–CaM complex indicated that some of the residues in the C- terminus of SidJ were crucial to the complex which is adjacent to the IQ motif (Fig. 6A). Our protein structural modeling showed a similar three-dimensional image as the truncated one (Fig. 6B). The definitive positive selection site on codon 58 of SidJ indicated that the mutation of this site was functional, although a previous study suggested that a truncated SidJ lacking the first 99 residues (SidJ(ΔN99)) showed activity indistinguishable from that of the full-length protein (Gan et al., 2019). This implied that a full-length protein of SidJ was required for a more detailed description of the function of some important amino acid sites of the whole protein. As shown in Fig. 6C, the G58E or G58M mutations significantly influence the three-dimensional structure of SidJ. The G58E mutation might be associated with the interaction with CaM because it was spatially closer to the I841 and Q842 than the original 58G (Fig. 6C). A similar explanation of the influence of codons 868 and 869 was determined as they were closer to those that could interact with CaM residues (Fig. 6B) and where the CaM docked and could mediate most of the interactions with CaM (Bhogaraju et al., 2019). Although codon 200 of SidJ had four mutation profiles, we found two of which might be functional. As shown in Fig. 6D, the N200T and N200V mutations significantly affected the hydrogen bonds of SidJ. Two hydrogen bonds of N200 could be formed with K196 and S204, while T200 could form three ones with K196, P197, and S204. In contrast, V200 could only form one with K196 (Fig. 6D). The more hydrogen bonds indicated a more stable α-helix, in turn, stabilizes the whole protein structure. Thus, the mutations in codon 200 may also mediate the interaction between SidJ and CaM and all these mutations in the positive selection sites could only adjust, but not abolish the affinity of CaM binding to SidJ. Based on protein structure modeling, we propose a potential explanation of the influence of the mutation on the four positive selection sites. Experimental data was still required to understand the exact function of these mutations. Mutation in positive selection sites might facilitate the survival of the lifeform containing mutated alleles. As an intracellular bacterium, entering the host and establishing infection were crucial for L. pneumophila lifecycle, and in which SidJ was dedicated to balance the host Ub ligase activity and was important for successful infection. L. pneumophila was shown to survive as an intracellular parasite of free-living protozoa in aquatic and moist soil environments (Bhogaraju et al., 2019; Fields, Benson & Besser, 2002; Gan et al., 2019). Protozoa provide a specific environment for gene exchange between L. pneumophila and other microorganisms invading them as pathogens or symbionts, also protozoa might be act as donors and transfer their own DNA to L. pneumophila (Mondino, Schmidt & Buchrieser, 2020). Many potential environmental hosts of L. pneumophila have been identified, including Dictyostelium discoideum (soil amoeba), Acanthamoebae castellanii, Entamoeba histolytica, Naegleria, and Tetrahymena, etc. (Hagele et al., 2000; Steinert et al., 1994; Rowbotham, 1980; Swart et al., 2018). Figure 7A showed the phylogenetic relationship of CaM among humans and potential hosts of L. pneumophila. Despite that CaM was relatively conserved, the protein sequences of CaM in these eukaryotes are variable (Fig. 7B, Table S6). And this might lead to a slightly structural difference among human CaM and those of protozoa hosts of L. pneumophila (Fig. 7C). Given that most of the time L. pneumophila is inhabiting in environmental hosts, we also propose that the variability in SidJ, especially that in positive selection sites might be important towards L. pneumophila survival in different environmental hosts and has an adaptive function to a broad selection of environments. The exact functions of these mutations to L. pneumophila living in the environmental host are worthy of further study. These results also suggested that more natural variants of a protein from a broad-host bacterium were required to discover the mechanism of how this protein and its variants are involved in the infection.

Conclusions

We presented molecular evolution analyses on a large and comparative set of sidJ alleles, derived from a collection of L. pneumophila strains. We found that among the 39 recognized sidJ alleles, about one-third were recombinants generated by eight PREs. Intragenic recombination also drove the sidJ diversification manifested by a higher genetic diversity in the recombinants as compared with that in non-recombinants. In addition, we found definitive positive Darwinian selection of SidJ at the codon level. Four codons in the NTD and CTD domains of SidJ were identified, and their mutation profiles were also determined. Protein structural modeling of SidJ provided possible functional explanations for the mutations in positive selection sites. It might influence the binding affinity of CaM to SidJ, thus regulate SidJ glutamylase activity to SidEs, and in turn balance the Ub ligase activity in different hosts. This study gave us a deeper understanding of the adaptive mechanisms of this intracellular bacterium to different hosts and highlighted the importance of the NTD and CTD domains in SidJ kinase activity that is activated by the binding of CaM. Further studies should focus on experimental evidence to illustrate the function and mechanism of natural sidJ mutants (with a mutation in positive selection sites) in regulating balanced Ub ligase activity in different L. pneumophila hosts.

Supplemental Information