Transkingdom Analysis of the Female Reproductive Tract Reveals Bacteriophages form Communities

Abstract

The female reproductive tract (FRT) microbiome plays a vital role in maintaining vaginal health. Viruses are key regulators of other microbial ecosystems, but little is known about how the FRT viruses (virome), particularly bacteriophages that comprise the phageome, impact FRT health and dysbiosis. We hypothesize that bacterial vaginosis (BV) is associated with altered FRT phageome diversity, transkingdom interplay, and bacteriophage discriminate taxa. Here, we conducted a retrospective, longitudinal analysis of vaginal swabs collected from 54 BV-positive and 46 BV-negative South African women. Bacteriome analysis revealed samples clustered into five distinct bacterial community groups (CGs), and further, bacterial alpha diversity was significantly associated with BV. Virome analysis on a subset of baseline samples showed FRT bacteriophages clustering into novel viral state types (VSTs), a viral community clustering system based on virome composition and abundance. Distinct BV bacteriophage signatures included increased alpha diversity along with discriminant Bacillus, Burkholderia, and Escherichia bacteriophages. Bacteriophage-bacteria transkingdom associations were also identified between Bacillus and Burkholderia viruses and BV-associated bacteria, providing key insights for future studies elucidating the transkingdom interactions driving BV-associated microbiome perturbations. In this cohort, bacteriophage-bacterial associations suggest complex interactions, which may play a role in the establishment and maintenance of BV.

Keywords: HIV; bacterial vaginosis; bacteriophage; female reproductive tract; microbiome; transkingdom associations; virome.

Figure 1

Bacteriome profiling by community group of self-collected vaginal swabs from South African women. (A) Relative abundance of the 16 most frequently identified bacterial taxa (y-axis) by sample (x-axis), grouped by community group (CG), bacterial vaginosis (BV) Status, visit number, highly active antiretroviral therapy (HAART) status, and human immunodeficiency virus (HIV) status (color key shown). Percent abundance is indicated by gradient key. Using Ward’s linkage hierarchical clustering, samples clustered into five distinct bacterial community profiles termed CG. (B) Bacterial composition (color key shown) for each of the five CGs (x-axis) expressed as relative abundance (y-axis). (C) Bar plot showing the relative abundance of 16S rRNA copies per 10 ng total DNA (y-axis) of L. iners (blue), L. crispatus (purple), L. gasseri (pink), and L. jensenii (green) bacterial species as determined by qPCR of FRT samples (x-axis) that clustered into CG1.

Figure 2

FRT bacteriome clusters into distinct community groups (CGs) that differ by alpha and beta diversity. (A) Bacterial Shannon diversity (y-axis) by CG (x-axis) as determined by linear mixed-effects model. Center bar represents median; gray box bounded by upper/lower interquartile ranges (IQR); whiskers represent range; dots represent outliers; color-filled areas are representative of density/distribution of diversity values. **, p < 0.01; ***, p < 0.001. (B) Principal coordinate analysis (PCoA) plot of the weighted UniFrac distances colored by CG (color key shown).

Figure 3

FRT DNA bacteriophages cluster into two unique community groups. Self-collected vaginal swabs were processed for DNA virome analysis by enriching for viral nucleic acid, libraries built and sequenced. (A) Relative abundance of the 32 most frequent bacteriophage species (y-axis) by sample (x-axis). Ward’s linkage hierarchical clustering analysis was used to cluster samples into distinct bacteriophage community profiles called viral state types (VSTs). VST, bacterial vaginosis (BV) status, highly active antiretroviral therapy (HAART) status, and human immunodeficiency virus (HIV) status (color key) are shown. Percent abundance is indicated by gradient key. (B) Bacteriophage Shannon diversity (y-axis) by VST (x-axis) as determined by linear regression model. Center bar represents median; gray box is bounded by upper/lower interquartile ranges (IQR); whiskers represent range; dots represent outliers; color-filled areas are representative of density/distribution of diversity values. ***, p < 0.001. (C) Principal coordinate analysis (PCoA) plots of beta diversity distances, as determined by permutational multivariate analysis of variance, colored by VST.

Figure 4

Transkingdom associations within the FRT of South African women. Heatmap of estimated Kendall’s correlation coefficients between FRT bacterial taxa (x-axis) and sequences assigned to bacteriophage (y-axis). Asterix indicate significant correlations after multiple comparisons correction by the Benjamini–Hochberg procedure, *, p < 0.05; **, p < 0.01. Magnitude and sign of the Kendall’s rank correlation coefficient are indicated by gradient key. Red indicates positive correlations; blue indicates negative correlations.

Figure 5

Discriminant FRT bacterial and bacteriophage species associated with bacterial vaginosis. (A) Significantly discriminant bacterial taxa by BV status were determined by univariate analysis using mixed-effects models. Relative abundance of each taxon (y-axis) is shown in BV-negative (orange) and BV-positive subjects (green; x-axis). All bacterial taxa presented are significant with an FDR-adjusted p < 0.05. (B) Bacteriophage Shannon diversity (y-axis) by BV status (x-axis) as determined by a linear regression model. *, p < 0.05. (C) Grouped bar plot showing the number of identified lytic and lysogenic bacteriophage contigs per sample by BV status. Lytic genes are represented as black triangles and lysogenic genes as black circles. Significance assessed using Mann–Whitney test with multiple comparisons correction by the Benjamini–Hochberg procedure. ***, p < 0.001; ****, p < 0.0001. (D) Significantly discriminant bacteriophage species by clinical BV diagnosis was determined by univariate analysis using a linear regression model. All bacteriophage taxa presented are significant with an FDR-adjusted p < 0.05.