Evolutionary Diversity and Structural Dynamics of the Outer Membrane Protein Ail in Yersinia

1. Introduction

Yersinia pestis is a category A select agent that can be found in populations of more than 200 species of wild rodents and was responsible for the deaths of more than 200 million people during multiple plague pandemics [1,2,3]. During evolution, different phylogenetic groups of this bacterial pathogen adapted to the effective circulation in the ecosystems of different natural plague foci due to the allele polymorphisms of the genes coding for some pathogenicity factors and housekeeping genes. Polymorphisms of genes encoding the synthesis of capsular antigen F1 (CaF1) [4,5], V antigen [6], plasminogen activator Pla [7,8], and lipopolysaccharide [9] have been studied in quite detail. However, much less attention has been paid to the Y. pestis Ail protein, which has been found to be an important virulence factor during plague infection [10]. The Y. pestis Ail is involved in bacterial autoaggregation [11], adhesion to and internalization into host cells [12,13,14], Yersinia outer proteins (Yop) delivery [11], serum resistance [15,16], and inhibition of inflammatory response [17].

It has been believed that the Ail protein is highly conserved within Y. pestis species [10], but it was recently shown that different isoforms of Ail are characteristic of different phylogenetic groups of the microbe [18] circulating in populations of different species of voles in natural plague foci located in different countries of Eurasia [18,19]. The straightforward way to perform their comparative study would be to collect strains with the all variants of the ail gene in one laboratory, construct their knockout mutants, and then cross-complement the mutations. However, while scientists can freely exchange genetic information with each other, exchange with strain specimens of category A select agent is limited [20]. Another way to conduct a comparative assessment of Ail isoforms is to synthesize allelic genes in vitro and use these synthesized genes to complement Δail mutations. Unfortunately, to construct representative sets of isogenic strains, including wild-type isolates with different ail alleles, their Δail derivatives, as well as cross-complemented ones, it is also necessary to overcome almost insurmountable export control restrictions. Thus, today in silico study of the properties of Ail isoforms offers the most accessible means of their comparative assessment.

In the current work, the primary sequences of Ail proteins from 84 strains of Y. pestis were aligned and compared in order to characterize their natural polymorphisms. Analysis of polymorphisms of Y. pestis ail gene with the help of multiple sequence alignment allowed us to identify nine natural isoforms of the ail gene in different Y. pestis strains. Then, these nine Ail isoforms were subjected to extensive molecular dynamics simulations and structural-dynamic characterization to understand the effect of sequence polymorphisms on intrinsic disorder propensity, structural and dynamic aspects of the Y. pestis Ail variants.

2. Results and Discussion

2.1. Evolutionary Aspects of Ail Polymorphism

Analyzing Y. pestis Ail polymorphism via multiple sequence alignment nine natural isoforms of Ail were found in different strains of Y. pestis showing a very high degree of sequence homology (98.32%). Figure 1 illustrates this observation by showing results of the multiple sequence alignment of these proteins using the MEGA11 [21].

Analysis of data shown in Figure 1 illustrates that Yersinia Ail can exist in nine isoforms:

1—The ancestral isoform (found in Y. pseudotuberculosis IP 32953 as well as in Y. pestis 0.PE7 and 0.PE4A);

2—The microti isoform caucasica (found in Y. pestis subsp. microti bv. caucasica, 0.PE2);

3—The microti isoform angola (found in Y. pestis strain Angola, 0.PE3);

4—The microti isoform 0.PE4B/ 0.PE4С var1 (found in Y. pestis subsp. microti bv. talassica, altaica and hissarica, 0.PE4B and bv. qinghaiensis (0.PE4C);

5—The microti isoform 0.PE4С var2 (found in some strains Y. pestis subsp. microti bv. xilingolensis, 0.PE4C);

6—The microti isoform 0.PE4С var3 (found in some strains Y. pestis subsp. microti bv. xilingolensis, 0.PE4C);

7—The pestis-microti isoform (found in Y. pestis subsp. microti bv. ulegeica, 0.PE5 and Y. pestis subsp. pestis except SNP- types 0.ANT1, 2.ANT);

8—The pestis isoform 0.ANT1 (found in Y. pestis subsp. pestis 0.ANT1);

9—The pestis isoform 2.ANT (found in Y. pestis subsp. pestis 2.ANT).

To further illustrate the peculiarities of Ail microevolution, the phylogenetic relationships of different Yersinia Ail isoforms are shown in Figure 2. This Y. pestis phylogeny scheme was adapted from the SNP-based phylogenetic tree generated using SplitsTree4 (http://www.splitstree.org/) [22]. Analysis of the Y. pestis whole-genome sequences available from GenBank/EMBL/DDBJ was performed using Wombac 2.0 (Bacterial core Genome SNPs for Phylogenomic Trees from NGS Reads and / or draft genomes; http://www.mybiosoftware.com/wombac-1–1-bacterial-core-genome-snps-phy logenomic-trees-ngs-reads-andor-draft-genomes.html).

Figure 2 shows that Y. pestis 0.PE7 has the identical Ail of Y. pseudotuberculosis IP 32953. Y. pseudotuberculosis IP32953 (O:1b) is traditionally used in Y. pestis phylogenetic studies as the closest out-group representative of the mutual progenitor belonging to O:1b [23]. Y. pestis 0.PE7 is the oldest currently existing branch of Yersinia pestis bacteria. Probably, the same Ail isoform had the ancestor of Y. pestis. The second branch seceded from Y. pestis tree is subsp. microti bv. caucasica, all of its representatives have Ail of the microti isoform caucasica characterized by additional S insertion at position 134. This change is in the extracellular loop 3 (EL 3) that is involved in binding vitronectin, factor H, C4BP, as well as maintaining the Y. pestis serum resistance phenotype (Figure 1) [11,24,25]. Third Ail isoform belongs to unique strain Angola of Y. pestis (0.PE3) and has changes F80L and R136I, K138E in extracellular loop 2 and 3. SNP-type 0.PE4 is the most polymorphic group of Y. pestis by the presence of Ail isoforms. Distal branch 0.PE4A has the identical Ail of Y. pseudotuberculosis IP 32953. SNP-groups 0.PE4B and 0.PE4C has polymorphism F80V. Y. pestis subsp. microti bv. xilingolensis acquired additional polymorphism. Part of them has D85G, another has D85Y change. All of these polymorphisms are in extracellular loop 2 of Ail (Figure 1 and Figure 2).

All of Y. pestis subsp. pestis and the most distant Y. pestis subsp. microti of the root bv. ulegeica has V126F change in the extracellular loop 3 of Ail. SNP-groups 0.ANT1 and 2.ANT have additional polymorphism. The pestis isoform 0.ANT1 has Y46D in extracellular loop 1 of Ail. Thes Y. pestis isoform 2.ANT has S128P in extracellular loop 3 of Ail.

2.2. Structural and Dynamic Aspects of Ail Variants

We now employ MD simulations to explore structural and dynamic aspects of the Ail variants. Each construct is simulated (in triplicate) in a bilayer membrane patch with lipid composition identical to the Gram-negative outer membrane model composition used in the experimental setup. Figure 3a shows the experimentally determined structure of the Y. pestis Ail protein [26], whereas Figure 3b shows the MD simulation box utilized in this study.

Figure 4a shows the root-mean-square deviation (RMSD) values for all the Ail protein constructs analyzed in this study. RMSD values derived from the individual MD simulation replicas are seen to exhibit some variation, for some constructs more than others. We use the RMSD values to quantify the extent to which the protein structure deviates from its initial starting point, with higher values reflecting greater deviation. The two strains that appear consistently deviating more than the rest based on mean RMSD values (across replicas) are 0.ANT1 and 0.PE4_Microtus_91001. Curiously, in comparison with the reference strain IP_32953, both of these are strain variants with multiple rather than single substitutions. To explore whether the mean RMSD values correlate with sequence deviation from the reference strain, we computed the percent-identity value after multiple sequence alignment against the reference strain, IP_32953. The results indicate correlation between mean RMSD and percent identity to the reference (Figure 4b). The Pearson’s r value for the correlation is -0.52, indicating moderate negative correlation (Figure 4b). Therefore, more severe mutations from the reference might lead to increased flexibility or alternate conformational states. Such structural flexibility is likely functionally significant as it can affect the protein’s ability to interact with extracellular molecules.

To further elucidate the possibility of varying structural flexibility among Ail variants, we computed the root-mean-square fluctuations (RMSF) from the MD simulations. Our analysis revealed that the core of the β-barrel structure is stable across all variants, with backbone RMSF values consistently near or below 1Å (Figure 5). This is noteworthy, since engineered variants of Ail usually aims to introduce mutations in the core of the β-barrel [27], in contrast to where natural variation is seen (Figure 1 and Figure 2). The extracellular loops exhibits significantly greater flexibility, as is evident from the RMSF line plots (Figure 5). The regions corresponding to the extracellular loops ECL1-4 are highlighted in the RMSF plots, allowing for comparative analyses.

We then asked whether the maximum RMSF values for individual ECLs can be connected to the evolutionary variations in their sequences. Discerning such connections directly from the RMSF line plots (Figure 5) is challenging. To address this, we performed a correlation analysis (Figure 6). When comparing the maximum RMSF values of each loop to the percent sequence identity relative to the reference strain, we see significant variation among strains. Yet, only some of this variation correlated with the sequence deviation. Specifically, we see that there is negligent or weak correlation between the maximum RMSF values for ECL1, ECL2, and ECL3, which have Pearson’s r values of 0.00, -0.20, 0.25, respectively.

Surprisingly, ECL4 displayed a strong negative correlation with the sequence identity, with a Pearson’s r of -0.77. This is unexpected, since the sequence of ECL4 is conserved across all variants tested (Figure 1). A plausible explanation is that the flexibilities and conformational dynamics of the other loops (ECL1-3) can influence the behavior of the ECL4. The smaller size of ECL4 further makes it more sensitive to such effects. This point is further emphasized by comparing the numerics of the maximum RMSF values e.g., by inspecting the heat-map plot (Figure 6b).

The text continues here (Figure 2 and Table 2).

We also explored the radius-of-gyration (R_GYR) for the Ail constructs to assess global structural differences induced by sequence polymorphism (Figure 7a). Despite the restricted nature of mutations in the Ail variants, a moderate positive correlation was observed between the sequence identity and average R_GYR values, with a Pearson’s r of 0.43 (Figure 7b). This suggests that even subtle sequence variations can influence the overall compactness of the Ail structure.

2.3. Intrinsic Disorder Aspects of Ail Variants

Intrinsic disorder and structural flexibility are known to have important functional implications, and all proteomes analyzed so far invariantly contain noticeable proportions of intrinsically disordered proteins (IDPs) and hybrid proteins containing ordered domains and functional intrinsically disordered regions (IDRs) [28,29,30,31]. IDPs/IDRs have multiple crucial biological functions that complement the order-based protein functionality [32,33]. Although due to their membrane surroundings, transmembrane proteins are typically considered as mostly ordered proteins, many of them are known to contain functionally important IDRs [34,35,36,37,38].

To understand the effect of sequence polymorphism on intrinsic disorder propensity of Yersinia Ail protein, a set of established disorder predictors was used in this study, such as members of the PONDR^® family [39,40,41] and IUPred [42]. Since for the individual Yersinia Ail isoforms the outputs of these disorder predictors generally agree (see Supplementary Materials, Figure S1), and since PONDR^® VSL2 [43] is one of the more accurate stand-alone disorder predictors [43,44,45], in the subsequent studies we utilized the outputs of the PONDR^® VSL2 predictor.

Figure 8a shows that the Yersinia Ail protein is predicted to have several IDRs (i.e., regions with the disorder scores above 0.5) and flexible regions (i.e., regions with disorder cores between 0.15 and 0.5) positions of which typically correspond to loops of the protein (cf. Figures 1, 5, and 8a). Amino acid substitutions found in different isoforms of the Yersinia Ail protein have some effect on the local intrinsic disorder propensity of this protein as shown by Figure 8a representing overlaid disorder profiles of nine isoforms. To better see the polymorphism-induced changes in intrinsic disorder propensity, Figure 8b represents “disorder difference spectra” calculated using the per-residue disorder outputs of the PONDR^® VSL2 as a simple difference between the per-residue disorder propensities calculated for individual Ail isoforms and the per-residue disorder propensities calculated for the ancestral Yersinia Ail isoform found in Y. pseudotuberculosis IP 32953 as well as in Y. pestis 0.PE7 and 0.PE4A. In such representation, positive peaks show regions, for which local intrinsic disorder propensity is increased relative to the ancestral isoform, whereas negative peaks shows regions with reduced local disorder propensity. Polymorphisms caused some increase in the overall disorder predisposition of 0.ANT1, 0.PE2, 0.PE3_Angola, 0.PE4B, and 0.PE4_Xilingolensis, whereas isoforms co92, 0.ANT, and 0.PE4_microtus_91001 are somewhat less disordered than the ancestral Yersinia Ail isoforms. In most cases, individual isoforms have major peaks of same sign. The exceptions are given by 0.PE2 and 0.PE3_Angola, which have noticeable positive and negative peaks.

Figure S1 shows that there is generally a good agreement between the disorder propensities predicted by PONDR^® VSL2 and structural flexibility evaluated in terms of RMSF. Here, regions with high RMSF are typically located within regions with high disorder score. Therefore, it seems that the structural flexibility of many Ail regions is somehow encoded in their intrinsic disorder propensity. Furthermore, major “dips” in the PONDR^® VSL2 plot typically correspond to the segments of predicted secondary structure, whereas predicted IDRs and flexible regions mostly coincide with the Ail loops. This is further evidenced by comparing the PONDR^® VSL2 disorder profiles (Figure 8a) and per-residue RMSF plots (Figure 8c). Curiously, Figure 8d shows that polymorphisms have more pronounced and diversified effects of local structural flexibility than on intrinsic disorder propensity. This plot contains “difference RMSF spectra” calculated by subtracting RMSF profile generated for the ancestral Yersinia Ail isoform from the RMSF profiles of individual isoforms. It is clearly seen that mutations cause both increase and decrease in the local flexibility.

Figure 8e represents the global correlation of the mean RMSF values calculated for all Yersinia Ail isoforms and their corresponding mean PONDR^® VSL2 scores, whereas Figure 8f shows correlation between the total RMSF differences of the individual Ail isoforms relative to the ancestral Yersinia Ail isoform and their corresponding total PONDR^® VSL2 score differences. These plots show a bit counterintuitive negative correlations, suggesting that mutations causing increase in the overall disorder propensities of isoforms are linked to the decreased overall structural flexibility.

Finally, to illustrate experimentally observed conformational dynamics of the of the Yersinia pestis outer membrane protein Ail, Figure 9 shows NMR solution structure of this protein in lipid bilayer nanodiscs and demonstrates that several its regions are characterized by high structural flexibility. This is evidenced by the structural fuzziness of the extracellular loops, where each of the included structural models have own structure non-overlapping with the loop structures of other models (see bottom side of Figure 9a,c). Therefore, these data indicate that the results of our computational analyses of the Ail structural dynamics are supported by the results of previously published NMR experiments [26].

3. Materials and Methods

Data on genomes of 84 Y. pestis strains representing main phylogenetic groups were used in this study as a source of amino acid sequences of Ail isoforms. They were SNP-groups, 0.ANT-4.ANT (17), 2.MED (5), 1.ORI (3), 1.IN (3) belonging to subsp. pestis; and 0.PE2 (18), 0.PE3 (1), 0.PE4 (29), 0.PE5 (5), 0.PE7 (2) from subsp. microti. The geographic position and background data for the phylogenetic groups were documented in previous publications [46,47,48]. First, sequence of the outer membrane protein adhesion invasion locus (Ail) from 84 strains, representatives of 17 Y. pestis phylogenetic groups, were aligned using the MEGA11 [21] to characterize their natural polymorphism. We also used Ail isoform of Y. pseudotuberculosis IP 32953. According to recent studies Y. pseudotuberculosis О:1b is the evolutionary ancestral of Y. pestis [23].

3.1. Multiple Sequence Alignment

Based on the sequence similarity, it was found that the Y. pestis Ail naturally exists in nine isoforms. First form is the identical to Ail of Y. pseudotuberculosis IP 32953 and are found in genomes of strains from 0.PE7 and 0.PE4A SNP-groups. Second isoform found in Y. pestis subsp. microti bv. caucasica (0.PE2). Third Ail isoform belongs to unique Y. pestis strain Angola (0.PE3). The strains of bv. hissarica / bv. altaica / bv. talassica (0.PE4B) and bv. qinghaiensis (0.PE4C) have Ail fourth isoform. Fifth and sixth forms of Ail have the representatives of Y. pestis subsp. microti bv. xilingolensis (0.PE4C). Seventh Ail isoform is the most represented among sequenced Y. pestis genomes. The strains belonged to Y. pestis subsp. pestis bv. antiqua (0.ANT2, 0.ANT3, 1.ANT, 3.ANT, 4.ANT), medievalis (2.MED), orientalis (1.ORI), and SNP-group 1.IN as well as Y. pestis subsp. microti bv. ulegeica have the identical isoform of Ail. Unlike others Y. pestis subsp. pestis SNP-groups 0.ANT1 and 2.ANT have different Ail isoforms, eighth and ninth forms respectively.

3.2. Molecular Dynamics Simulations

Constructs corresponding to the nine unique isoforms were characterized using molecular dynamics simulations. The Yersinia pestis outer membrane protein Ail resolved by NMR in phospholipid bilayer nanodiscs was used as structural template (PDB ID: 5VJ8, residues 27-182) [26] with the mutations corresponding to the various isoforms introduced in silico. The proteins were embedded in a bilayer membrane with a lipid ratio equivalent to that of used in the NMR study [26], namely a DMPC:DMPG ratio of 75:25. Protonation states were determined using PROPKA-3.0 [49] and the relative orientation of the proteins were found by using PPM 3.0 [50]. The systems were subsequently solvated by TIP3P water [51] and 150 mM KCl salt using the CHARMM-GUI builder [52,53]. The CHARMM36m force field [54] was using to describe proteins, lipids and ions. The standard CHARMM-GUI protocol for minimization and relaxation were employed. Each of the nine constructs were then subjected to molecular dynamic (MD) simulation for 3 × 600 ns (triplicate) in the NPT ensemble by initializing independent velocity distributions, totaling 27 systems and 16.2 μs. The target pressure (1 atm) and temperature (T=303.15 °K) were maintained using a Monte Carlo barostat [55] (coupling frequency of 100 steps) and the Langevin thermostat [56] (friction coefficient of 1 ps^-1), respectively. Long-range electrostatic forces were calculated using the particle mesh Ewald (PME) method with a grid spacing of approximately 1 Å and 5th-order B-splines for interpolation [57]. Electrostatic forces were switched starting at 10 Å and switched entirely off by 12 Å. Periodic boundary conditions were applied in the X-,Y-, and Z-directions. All MD simulations were performed using OpenMM 7 [58]. An integration time step of 2.0 fs was used for the velocity Verlet algorithm. The lengths of all bonds that involve a hydrogen atom were constrained. Figure S2 in Supplementary Materials represents Root Mean Square Deviation (RMSD) for the various Ail constructs across three replicates of MD simulation, whereas Figure S3 shows corresponding Root Mean Square Fluctuation (RMSF) for the various Ail constructs.

3.3. Evaluation of the Effect of Polymorphism on the Intrinsic Disorder Predisposition

The effect of polymorphism on the intrinsic disorder predisposition of the Yersinia outer membrane protein Ail was evaluated out using a set of predictors from the PONDR^® family, PONDR-FIT, PONDR^® VLXT, and PONDR^® VSL2 servers [39,40,41] as well as by the IUPred algorithm [42]. The outputs of these tools were combined using the RIDAO web platform that generates an integral disorder profile of an individual query protein [59]. Since the results of all predictors were generally in a good agreement with each other, the effect of natural polymorphism on the local intrinsic propensity of the Yersinia outer membrane protein Ail was investigated using PONDR^® VSL2. We also analyzed the effect of sequence polymorphism on intrinsic disorder propensity of this protein using PONDR^® VSL2 in the form “difference disorder spectra”, which were calculated by subtracting the per-residue disorder propensities calculated for the ancestral Yersinia Ail isoform found in Y. pseudotuberculosis IP 32953 as well as in Y. pestis 0.PE7 and 0.PE4A (see below) from the per-residue disorder propensities calculated by the PONDR^® VSL2 for eight different Yersinia Ail isoforms described below.

4. Conclusions

The currently existing strains of Y. pestis can be divided into two groups that differ fundamentally in their pathogenicity for humans [19]. The phylogenetically more ancient and conditionally pathogenic for men representatives of the subspecies microti are characterized by a significant diversity of amino acid sequences of their pathogenicity factors such as capsular antigen F1 [4,5], V antigen [6], plasminogen activator Pla [7,8], and lipopolysaccharide [9]. A similar situation is observed in the case of Ail, where comparative analysis of sequence variants across 81 Y. pestis strains revealed polymorphisms located in the extracellular loops. These polymorphisms were correlated to measurable differences in dynamic and structural properties as revealed by molecular dynamics simulations and intrinsic order predictions. It could be speculated that the evolution of Ail into the pestis-microti isoform, as well as into the pestis isoforms 0.ANT1 and 2.ANT, is the main reason for the change in the host specificity of this phylogenetic group of Y. pestis towards pathogenicity for humans. However, it is more likely that the high human pathogenicity of the pestis subspecies strains is the result of a change in the structural and functional organization of not just one protein, but a whole complex of pathogenicity factors. Thus, Ail variation likely represents but one side of a broader evolutionary adaptation shaping host specificity in Y. pestis.