The core genome evolution of Lactobacillus crispatus as a driving force for niche competition in the human vaginal tract

Abstract

The lower female reproductive tract is notoriously dominated by Lactobacillus species, among which Lactobacillus crispatus emerges for its protective and health-promoting activities. Although previous comparative genome analyses highlighted genetic and phenotypic diversity within the L. crispatus species, most studies have focused on the presence/absence of accessory genes. Here, we investigated the variation at the single nucleotide level within protein-encoding genes shared across a human-derived L. crispatus strain selection, which includes 200 currently available human-derived L. crispatus genomes as well as 41 chromosome sequences of such taxon that have been decoded in the framework of this study. Such data clearly pointed out the presence of intra-species micro-diversities that could have evolutionary significance contributing to phenotypical diversification by affecting protein domains. Specifically, two single nucleotide variations in the type II pullulanase gene sequence led to specific amino acid substitutions, possibly explaining the substantial differences in the growth performances and competition abilities observed in a multi-strain bioreactor culture simulating the vaginal environment. Accordingly, L. crispatus strains display different growth performances, suggesting that the colonisation and stable persistence in the female reproductive tract between the members of this taxon is highly variable.

FIGURE 1

Comparative analysis between different Lactobacillus species and between 22 non‐identical L. crispatus strains. In panel (A), Venn diagram shows the eight Lactobacillus species sharing the 159 genes used to measure the magnitude of intra‐species genetic diversity. Below, the species‐specific number of SNPs identified within the common protein‐encoding genes is reported for each of the considered Lactobacillus species. Panel (B) depicts the L. crispatus pan‐ and core‐genome size. The number of discovered genes (vertical axe) is plotted as a function of the number of sequentially added genomes (horizontal axe). Panel (C) shows the phylogenomic tree based on the concatenated 959 core genes shared among the 22 non‐identical L. crispatus genomes. The tree was constructed by the neighbour‐joining method. Bootstrap percentages based on 1000 replicates above 50 are shown at node points. For each strain, the isolation source is highlighted with a coloured circle. On the right, an aligned portion of the L. crispatus core genome exemplifies the relationships between phylogenomic clusters and SNP patterns. In the top row, nucleotide positions showing variants are highlighted with an asterisk, while a dot highlights non‐variant sites.

FIGURE 2

Identification of the 52 HVGs. In panel (A), Box‐Whisker plot was used to represent the gene distribution based on the number of SNP sites obtained by comparing the nucleotide sequence of every 959 protein‐coding genes shared among all the 22 Lactobacillus crispatus chromosomes. For each gene, the number of SNP sites was expressed as the average of all the pairwise comparisons against the reference sequence (homologous gene sequence of L. crispatus GCF_015669875). The Q3 + 1.5IQR was used as a cut‐off to select the 52 HVGs. Panel (B) reports the functional annotation and the number of SNP sites of each HVG.

FIGURE 3

Phylogenetic analysis based on the 52 HVG sequences. Proteomic tree based on concatenating the 52 protein encoding core genes identified as highly variable across the 22 non‐identical Lactobacillus crispatus genomes. Phylogenetic groups are highlighted in different colours. For comparison, the phylogenomic tree resulting from the whole core genome (presented in Figure 1) is visualised using the redial layout. For each strain, the coloured circle represents the isolation source, while the diameter of the black circle is proportional to the number of SNPs identified within the core genome using the GCF_015669875.1 genome sequence as reference.

FIGURE 4

Differential growth and competitive abilities between Lactobacillus crispatus strains. Panel (A) shows the optical density (OD) registered after 48 h of anaerobic growth in different nutritive substrates. Panel (B) illustrates the design of the bioreactor‐based experiment simulating the vaginal environment. In panel (C), the bar chart reports the quantification of the metagenomic reads (using average RPKM measures) mapping marker genes unique to each L. crispatus strain throughout the 48 h of growth in the bioreactor. The standard deviations are plotted as error bars. Different lowercase letters indicate significant differences at p‐value < 0.05 according to the Bonferroni test. In panel (D), alignment of partial amino acid sequences corresponding to the type II pullulanase genes of L. crispatus LB97 and LMG11440 strains genes highlights two amino acid substitutions Gln (Q) to Lys (K) and Arg (R) to Thr (T).