Common cervicovaginal sequencing methods result in discordant molecular diagnoses of bacterial vaginosis and reveal strain level effects of Gardnerella vaginalis.

Samples and clinical data were collected between January 2019 and January 2024 in Miami-Dade County, Florida as part of a larger longitudinal study to investigate the behavioral and biological factors influencing recurrent BV and HIV risk among cisgender women of reproductive age as part of the “Women, HIV, Immunology, Microbiome and Sexual Health” (WHIMS) cohort. HIV-negative, cisgender woman aged 18–45 years reporting sexual activity within the previous three months were enrolled. Exclusion criteria were pregnancy, antibiotic use, or STI diagnosis in the last 2 months, immunosuppressed status, or previous cervical surgery. Participants were recruited through consent-to-contact registries, health fairs, and public flyers in community clinics and bus stops, in collaboration with the Miami Center for AIDS Research (CFAR). Candidates were screened over the phone before being invited to complete an in-person enrollment visit. The study was conducted at the University of Miami Infectious Diseases and CFAR clinical Core research unit. During the in-person visit, participants completed a written informed consent, questionnaires, and underwent a vaginal exam. During the vaginal exam, vaginal swabs were collected in the following sequence (1. for Gram stain and Nugent scoring, 2. for Amsel criteria including clue cell count and amine production test (“whiff test”), 3. for microbiome analysis). Participants who reported menstrual bleeding or spotting, or who had blood noted during the vaginal examination were requested to reschedule their enrollment clinical exam and testing. This study was conducted in accordance with the ethical principles of the Declaration of Helsinki. Ethical approval was obtained from the University of Miami Institutional Review Board (protocol number 20180758), and all participants provided written informed consent prior to any study related assessments. Follow-up assessments included a questionnaire, rapid HIV and pregnancy testing, and a gynecological examination to evaluate clinical BV criteria. Diagnosis of BV for enrollment was determined by Amsel’s criteria. BV was diagnosed by Amsel’s criteria when three out of four Amsel’s criteria labs were positive (1. Thin milky discharge present; 2. Vaginal pH > 4.5; 3. Greater than 20% clue cells proportion on a vaginal swab slide preparation; 4. Production of amine odors after application of KOH to a vaginal swab preparation, the “whiff test”). For Nugent scoring, a vaginal discharge slide preparation was Gram stained and evaluated by a laboratory technician. The Nugent score of each participant and the grouping of the score (0–3 = BV negative, 4–6 = BV intermediate, and 7–10 = BV positive) was determined. Participants who met Amsel’s criteria for BV were prescribed 500mg metronidazole perorally twice daily for 7 days. Vaginal swabs and cervical cytobrush samples were collected during enrollment and during at least two follow-up clinical visits at 1 and 6-months post-enrollment. Participants were followed at 1 month and 6 months and tested for BV by Amsel criteria at each time point. Participants who tested positive for BV by Amsel’s criteria at their 1-month clinical visit were prescribed an additional 7-day course of metronidazole (500mg perorally twice daily) and assessed at an additional clinical visit 2 months after initial enrollment ( Table 1 ). Out of 269 samples (n = 121 participants) in the parent study, our final analysis included 88 samples from 44 participants from all visits ( Table 1 ). For this analysis, we included only samples which had high quality sequencing results for metataxonomic, metagenomics, and metatranscriptomics each recapitulates the microbiome in a different way. We assessed all three to determine how they relate to each other each of the three approaches addressed. Total DNA from vaginal swabs was extracted using the PowerSoil Pro kit HT (Qiagen, Valencia, CA) per manufactures instructions. Extracted DNA was amplified using a dual-indexing qPCR method for metataxonomics utilizing the 341F-806R primers to target the V3-V4 region of the 16S SSU rRNA. This method is a two-step qPCR method where the 16S-specific primers and adaptor tails for adding indices are added in the first amplification and the Illumina flow cell adaptors are added in a secondary amplification 58 . Amplified samples were then normalized, pooled, and sequenced on an Illumina MiSeq platform with 15% phiX with a read length of 300 base pairs. A water control was extracted simultaneously with study specimens to monitor for contamination during the library preparation. To remove the adaptor and primers from resulting sequencing files the Cutadapt tool was utilized. The DADA2 R package pipeline was then used to filter, trim, check quality of sequences, remove chimeras and produce an amplicon sequence variants (ASV) table using a default trimming parameters and Phred quality threshold of 30 59 . The Silva database version 138 was used for taxonomic classification 60 . DNA from the same total DNA vaginal swab extractions used for metataxonomics was used for metagenome sequencing. Quarter reactions using the Illumina Nextera XT Library System were used for DNA amplification, quantification, and normalization (Illumina, San Diego, CA). The AmpPure clean up kit (Beckman Coulter, Brea, CA) was used on all samples to increase sequencing efficiency. Sequencing was carried out on the Illumina NovaSeq platform targeting a mean depth of least 7 million 150bp paired end reads per sample. Human reads were removed using BMtagger 61 , and quality control and trimming were done with fastp 62 using default quality metrics and a Phred cutoff of 30. Paired-end reads were mapped to the VIRGO2 vaginal microbiome gene catalogue. Cervical cytobrush samples were used for RNA extractions to increase RNA yield, and previous research has shown correspondence between the vaginal and cervical microbiome 63 . RNA was extracted from cervical cytobrush dry cell pellets stored in RNAlater using the Qiagen QIAamp HT RNeasy kit (Qiagen, Valencia, CA) per manufactures instructions. Ribosomal RNA depletion using the RiboZero Plus (Illumina, San Diego, CA) kit preceded library preparation using the Illumina NovaSeq platform with a S4 2×150-bp flow cell. Trimmomatic was used to remove adaptors and filter low-quality reads with a Phred score below 30. We targeted 2,250 million pass filter reads for each lane and 40 million reads per sample. Samples with fewer than 3 million reads that passed the quality filtering were removed from the analysis 64 . Filtered reads were mapped to the vaginal microbiome gene catalog VIRGO2 (an update to VIRGO 16 provided by the Jacques Ravel group in advance of publication), and mapped reads were imported into R using tximport 61 to generate a length-scaled transcripts per million (TPM) data matrix for downstream analysis. For both the metagenome and transcriptome, a python script was used to convert gene count data into a composition matrix for CST classification using Valencia. CST and subCST classification for each sequencing approach was performed using Valencia 17 . SubCSTs split the original CST III and CST IV into lettered subgroups to accommodate the greater variation observed in those CSTs. All statistical analyses were performed in R (version 4.3.0, Already Tomorrow ). The vegan 65 package was used for PERMANOVA, Shannon, and Simpson diversity calculations. The rstatix 66 package was used for the Wilcoxon, Kruskal-Wallis, and Dunn statistical tests used in analyzing Shannon and Simpson diversity differences between BV status groups by Amsel criteria and Nugent score groups within each subCST. The umap 67 (Uniform Manifold package was used for dimensional reduction plots. Relationships between CSTs, individual bacterial species relative abundance (metagenome and metataxonomic) or relative contribution (metatranscriptome) to the microbiome, and individual Amsel criteria were evaluated by fitting Bayesian logistic regressions using the brm function from the brms 40 package in R. Models contained either CST and Amsel criteria data, or bacterial species abundance and Amsel criteria data. The Bernoulli distribution was used in models assessing the binary Amsel’s criteria clinical test results, and the categorical family was for CST models. Beta distribution with a logit link was used for models assessing the relationship between Nugent score and bacterial abundance with Nugent score scaled to fall between 0 and 1 to conform with the distribution assumptions. Each model was stratified by participant ID as a fixed effect to control for multiple sampling and demographic effects by allowing both intercept and slope to be determined for each participant. Model robustness was assessed using R-hat and effective samples size (ESS). Twelve Markov Chain Monte Carlo (MCMC) iterations were employed in each model with 500 warm-up samples and 2500 additional iterations. The Adapt delta parameter was set to 0.95 and max tree depth was set to 25. To avoid overfitting, all Amsel and CST models were run with heavy tailed student-t priors with 3 degrees of freedom, a location of 0, and a scale of 2.5. Similarly, the normal priors were used for Nugent score models with an intercept of 0 and scale of 5. For all non-Bayesian statistical analyses, alpha was set to 0.05 and the Benjamini-Hochberg 68 correction method was used for false discovery rate (FDR) for all multiple comparison testing. All plots, model summaries and tables were created in R utilizing the ggplot2 69 and gt 70 packages. dplyr 71 and tidyverse 72 were used for data formatting and utility functions. Nearly all Bayesian models had an Rhat (a measure of model convergence with 1.0 being perfect convergence) less than or equal 1.001 indicating excellent convergence. Two models assessing the effect microbes on the incidence of positive whiff tests had Rhat statistics above 1.001; the metagenomic based model assessing the effect of Fannyhessia on the whiff test, and the metagenomic model assessing the effect of all Lactobacillus on the whiff test. The Rhat for the metagenomic model assessing the effect of the CST II community structure on the discharge lab was also above 1.001 (1.01). For models comparing individual species relative abundance in the metagenome or relative contribution to the metatranscriptome with individual Amsel’s criteria, the estimated effect of a one-unit increase (1% relative abundance or contribution) in the predictor on the odds of a positive lab outcome was reported. For the models assessing the relationship between individual species and Nugent score, the estimated effect of a one-unit increase (1% relative abundance or contribution) in the predictor on expected Nugent score was reported. For the models assessing the relationship between CSTs and Amsel’s criteria results, the change in odds of a positive criteria compared to CST I were reported ( Figures 3 – 8 and Supplemental figures 3 – 8 ).

To compare the dominant bacteria identified by each sequencing approach, we calculated the top 20 species’ relative abundance for each sequencing method ( Figure 1 ). Only four species were shared among the top 20 species across each of the three sequencing approaches, and only six were found in both DNA-based sequencing approaches: metataxonomic and metagenomic sequencing. Thirteen species were shared between the top 20 species in the metagenomic and metatranscriptomics results. The highest abundance species identified across all sequencing methods was L. iners ( Figure 1 ). Our analysis of the intra-CSTs variation in participants’ BV status by Amsel and Nugent scores identified significant discordance between molecular and clinical assessments, with women characterized as BV negative by Amsel’s criteria and Nugent scoring assigned to the same BV-positive CSTs (CST IIIB, and IV). However, based on metataxonomics we did not identify significant differences in overall diversity within CSTs based on Amsel and Nugent status as measured by Wilcoxon and Kruskal-Wallis tests ( Supplemental figure 1 ). We did find that within metagenomic-defined CST IIIB, Amsel BV-positive participants showed significantly higher Shannon diversity than those who were Amsel BV-negative ( p adj = 0.0238 and p adj = 0.042, respectively; Supplemental figure 1 ). CST IIIB Amsel BV-positive participants also exhibited significantly lower Simpson diversity than those who were Amsel BV-positive ( p adj = 0.0303; Supplemental figure 2D ). Comparison of Nugent score groups for Shannon and Simpson metrics in the metagenomic data did not reveal significant differences between the CST groups. No significant differences in diversity within CSTs by Amsel’s criteria or Nugent score groups were identified in the metatranscriptomics data when examining Shannon or Simpson metrics. However, while not reaching significance, the same directional trend was observed in the metatranscriptome and metagenome; that substantial variability exists within molecular CSTs when separated by clinical diagnosis ( Supplemental figure 1 – 2 ). Our assessment of the interactions between microbiome structure and Amsel and Nugent criteria using PERMANOVA analysis based on UMAP distance matrices revealed significant microbiome variation by both Amsel’s criteria and Nugent score group identity ( Figure 2 ). For the metataxonomics data, Amsel’s criteria groupings were significantly different overall, explaining an estimated 19% of total ordinal variation ( p =0.001). Pairwise comparison of Nugent score groups showed significant differences in overall bacterial composition: Nugent groups 4–6 were 8% different compared to Nugent groups 0–3 (p = 0.001); Nugent groups 4–6 were 5% different compared to 7–10 (p = 0.007); and Nugent groups 0–3 were 22% different from groups 7–10 (p = 0.001) ( Figure 2 ). Amsel criteria and Nugent score groups were also modeled with CST as an interaction effect. No significant interaction effect was found between Amsel’s criteria and CST. With Nugent score groups, there was a significant interaction effect when comparing groups, except for Nugent score group 0–3 compared to 4–6. Metagenomic data PERMANOVA analysis showed similar results, with Amsel-BV status having a significantly effect and explaining and estimated 18% of the compositional variation (p=0.001). A significant interaction effect was found between Amsel criteria and subCST (p=0.022). Nugent groups 4–6, compared to 0–3, explained 7% of the variation (p = 0.001); Nugent groups 4–6, compared to 7–10, explained 5% of the variation (p = 0.004); and Nugent groups 0–3, compared to 7–10, explained 24% (p = 0.001) ( Figure 2 ) No significant interaction effect was found between any comparisons of Nugent groups and metagenomic subCSTs. Amsel-BV status samples (BV+ vs BV−) had significantly different metatranscriptomic compositions, with Amsel-BV status accounting for 23% of the total variation (p=0.001). Like the metataxonomic data, no significant interaction effect was found between Amsel-BV status and subCST. The only significant interaction effect with subCST when performing comparisons was for Nugent groups 0–3 compared to 7–10 (p=0.002). When accounting for subCST in the models, the p -values were all significant for subCST (p=0.001). SubCST in all models explained a larger portion of the variations, ranging from 49% to 72%, depending on the model ( Figure 2 ). Our analysis assessing the concordance of CST assignments across three sequencing methods are reported as the percentage of matched and mismatched CST assignments for each sample between sequencing strategies in a pairwise manner ( Supplemental table 1 ). The comparisons revealed that the percentage of matching assignments between metagenomic and metataxonomic CSTs and subCSTs were 84% and 66%, respectively. Comparing assignments based on the metatranscriptome to the metagenome or the metatranscriptome to metataxonomic CSTs demonstrated slightly lower concordance (65% and 59%; and 68% and 59%, respectively). Notably, the majority of subCST mismatches between the RNA-based and the DNA-based methods (metatranscriptomics vs. metataxonomics and metagenomics) were identified as CST IV-C classifications in the metatranscriptome shifting to subCST III in the two DNA-based assignments. In summary, we found that sequencing methodology affects CST and subCST assignments and implied molecular BV status. Unsurprisingly, the discordance is more pronounced when comparing cervical RNA and vaginal DNA-based sequencing, despite their similarities. Models assessing the relationship between individual microbes identified in the metagenome and metatranscriptome and the result of the Amsel’s criteria amine odor (“whiff test”) lab showed strain-specific relationships within the G. vaginalis species ( Figure 3 ). In each case, G. vaginalis strain B was one of the top three predictors of a positive whiff lab, substantially higher than G. vaginalis overall. In both the metagenome and metatranscriptome, the estimated effect of Dialister and Parvimonas was also substantial ( Figure 3 ). This is consistent with past research demonstrating substantial differences between metagenomic composition and metatranscriptome contributions. 20 Models assessing the effect of microbes on the likelihood of a participant having a pH above 4.5 showed that G. vaginalis strain B was a strong predictor of elevated pH in the metatranscriptome and the metagenome ( Figure 4 ). Prevotella amnii and Dialister similarly had high estimated effects on the likelihood of a pH >4.5. The estimated effect of P. micra on the pH measurement based on the metatranscriptome data was also substantial, indicating with moderately high likelihood that P. micra is associated with vaginal pH changes in our cohort ( Figure 4 ). Gardnerella vaginalis strain B also had the highest estimated impact on the likelihood of having clue cells based on the metagenome, and G. vaginalis strain B was the third strongest predictor of a positive clue cell test based on the metatranscriptome, behind P. micra and D. micraerophilus ( Figure 5 ). Additionally, P. bivia , D. micraerophilus , and P. micra had high estimate effects on the likelihood of a positive clue cells ( Figure 5 ). Nearly all 95% credible intervals estimating the effect of bacteria on the likelihood of abnormal vaginal discharge overlapped zero, resulting in a low confidence of the direction and specificity of the associations of each species. The exceptions to this trend were L. gasseri and P. bivia in both the metagenome and metatranscriptome, and L. jensenii in the metagenome, which were associated with abnormal discharge being reported in the Amsel’s criteria ( Figure 6 ). The models assessing the effect of the relative abundance and metatranscriptome contribution of individual bacterial species on Nugent score most closely mirrored the results of the modeling using the clue cell Amsel’s criteria. The association of Gardnerella on with Nugent score, both at the genus level and specifically G. vaginalis strain B, was less pronounced than the association of Gardnerella with individual Amsel’s criteria. The estimated effect had wider 95% confidence intervals that overlapped zero, suggesting relatively low confidence in the estimated effect Gardnerella specifically on Nugent score. The species with the highest and most consistent estimated effect on Nugent score were Dialister and P. micra ( Figure 7 ). Modeling the association between metagenomic, metatranscriptomic, metataxonomic CSTs with the likelihood of individual positive criteria demonstrated that sequencing type influenced how individual CSTs were associated with Amsel’s criteria ( Figure 8 ). Interestingly, compared to CST I, CST III was more associated with an elevated pH and positive whiff test, but less associated with the presence of clue cells and abnormal discharge. In the metagenome, CST IVA was the CST most associated with a positive whiff test and clue cells. CST IVC was the most associated with elevated pH and CST II was most associated with abnormal vaginal discharge. However, each estimate distribution for abnormal vaginal discharge overlapped zero, indicating very low confidence that abnormal discharge is strongly associated with any metagenomic-based CST. In the metatranscriptome, CST IVB was the CST with the highest association with a positive whiff test and an elevated pH, CST V was most associated with the presence of clue cells, and CST II was most associated with abnormal vaginal discharge. For CSTs established using bacterial composition generated by metataxonomics, CST II and IVB were most associated with a positive whiff test, CST V was most associated with the presence of clue cells, and CST II was most associated with abnormal vaginal discharge and an elevated pH ( Figure 8 ).

Despite decades of research, bacterial vaginosis affects an unacceptable number of women worldwide and poses a significant risk for morbidities including STI and HIV acquisition and transmission 11 . This study evaluated BV using a multi-omic dataset aiming to increase our understanding of how different sequencing methodologies can influence conclusions that are drawn from microbiome studies about the prevalence of clinical and molecular BV and their level of correspondence. Additionally, we used those data to elucidate how individual species or CSTs influence the likelihood of positive results for clinical BV using common clinical criteria. The etiology of BV remains poorly defined and BV research and clinical care continues to evaluate BV both as a condition defined by G. vaginalis as an etiological agent and as a condition defined by changes in the vaginal microbiota and the presence of non-optimal, anaerobic, polymicrobial community structures 14 , 21 , 28 , 32 , 33 . In this study, we sought to contribute data connecting and contrasting those two definitions. Comparing the Shannon and Simpson metrics we calculated provides insight into microbial community diversity within each CST ( Supplemental figures 1 – 2 ). While modest samples size limited the number of statistically significant associations we found, we identified intra-CST variation in CST III, IV, and V. This finding highlights the fact that even with the precision of molecular profiling, enough variation remains within CST groupings to encompass both Amsel and Nugent negative participants in classically BV positive molecular CSTs, namely CST IIIB, and IV. While this was not surprising, it demonstrates again the complex associations between clinical and molecular assessments of BV, reinforcing the idea that BV is not one discrete condition, but rather a collection of related conditions grouped together out of necessity. Likewise, Amsel’s criteria, the only assessment available for immediate clinical diagnosis, do not always align with CSTs or molecular microbial diversity. Recent studies have already begun to increase the granularity of molecular profiling in part to address this disparity 21 . Understanding the differences in microbial diversity within CSTs based on clinical diagnosis could aid in refining diagnostic criteria or and inform the development of more personalized therapeutic strategies for different versions of BV. Within the metagenomic data, participants with Amsel’s BV in CST III-B and IV-B showed significantly higher Shannon and Simpson diversity, suggesting increased microbial richness even within a CST could potentially lead to a higher risk of symptom onset and diagnosis of BV by Amsel criteria. These results align with previous studies indicating that microbial composition plays a crucial role in the development and progression to a non-optimal vaginal microbiome 34 . Metagenomic sequencing yielded the only statistically significant intra-CST variation, demonstrating the value of sequencing beyond metataxonomic resolution. While metatranscriptomics data provide increased granularity, they measure bacterial activity rather than simple abundance. The metatranscriptome data in this study were also generated from cervical cytobrush samples. Combined, these factors contribute to the lack of concordance between metatranscriptome and metagenome results. Indeed, differences between the vaginal metagenome and metatranscriptome have been previously demonstrated 20 . Thus, some discordance between the two methods was expected, and continuing to parse the specific ways that they differ will help elucidate the nuanced differences between different types of BV. We also noted discordance in subCSTs between different sequencing methods (metagenomics, metataxonomic, and metatranscriptomics, Supplemental table 1 ). Our results revealed that while there was a notable concordance between metagenomics and metataxonomics methods, the level of agreement was not as high as anticipated, especially considering both sets of sequencing data were generated from the same sample set and same DNA extractions. This discordance is likely attributable to the inherent imprecision of 16S taxa databases compared to the specificity of the vaginal shotgun database (VIRGO2) we were able to use for metagenomic and metatranscriptome classification. Previous studies have also shown inconsistencies in taxonomic assignments based on different databases 35 , 36 . Nonetheless, although concordance at the CST level between shotgun metagenomics and metataxonomic sequencing reached 84%, when examining the subCST level, this concordance dropped to 66%. This difference underscores the importance of considering finer taxonomic resolutions, as discrepancies at the subCST level suggest variations in microbial community structure that could have important functional implications. These results also demonstrate that the choice of sequencing technology alone could have substantial influence on molecular BV characterization and diagnosis. Metatranscriptomics exhibited lower concordance with both metagenomic and metataxonomic methods ( Supplemental table 1 ). While we acknowledge that the metatranscriptomics data came from a different sample type (cervical vs vaginal), the discrepancy in concordance nevertheless emphasizes the importance of considering the choice of sequencing approach in characterizing the vaginal microbiome, as different methods capture microbial community structures in different ways (community composition vs community contribution) and with varying degrees of fidelity and granularity. Most of the discordance between CSTs calculated from metatranscriptomics and CSTs calculated from metataxonomic sequencing were instances when the metatranscriptomics CST was IV, while the metataxonomic and metagenomic CSTs were classified as CST III. Of note, CST III is dominated by L. iners , which has been associated with transition to BV and makes the L-isoform of lactic acid, which lowers pH less robustly then D-isoform that the other vaginal Lactobacillus species produces 37 . Our results contribute to the body of evidence suggesting that metatranscriptomics, since it measures gene expression, infers bacterial “active abundance” and rather than the simple relative abundance captured by metagenomics sequencing 20 . The fact that CSTs based on bacterial activity differ from CSTs based on bacterial abundance alone highlights the potential that species with relatively low abundance could exert a significant influence on BV dynamics because they are disproportionately active. Despite limited past research, additional studies are needed that examine the vaginal microbiota at the expression level to fully elucidate this phenomenon 20 . Our modeling results were generally consistent with the consensus understanding of BV as a condition significantly influenced by the relative abundance of Lactobacillus species, particularly L. crispatus and L. iners , and the relative abundance of G. vaginalis and other opportunistic anaerobes 12 , 20 , 27 , 38 ( Figures 3 – 8 ). We chose to model the associations between individual bacteria and specific vaginal symptoms using a Bayesian approach to take advantage of the power and flexibility of the brms 39 package. Since the Bayesian statistical approach doesn’t produce p-values or explicitly calculate significance in the way a frequentist approach does, we reported credible intervals of each estimated posterior distribution 40 . For convenience we will refer to results with an estimated 95% credible interval that doesn’t overlap zero as “credible” effects hereafter. Our results showed that there is a credible association between G. vaginalis , particularly G. vaginalis strain B, and a positive result in each individual Amsel’s criteria lab apart from abnormal discharge. Gardnerella overall and G. vaginalis strain B were also credibly associated with an elevated Nugent score. That pattern was pronounced compared with the effect of other species and even other strains of G. vaginalis . Parvimonas micra was only identified in the metatranscriptome but was also credibly associated with a positive result in each of the individual Amsel’s criteria labs except discharge, which has previously been shown to have low diagnostic utility 41 . P. micra was also credibly associated with an elevated Nugent score, underscoring the similarity between the associations between P. micra and G. vaginalis and Amsel’s criteria in the metatranscriptome. Parvimonas has received minimal attention in the context of BV, but has been shown to be an important trigger of inflammation and has been reported in isolated cases of bacteremia 42 – 45 . These data suggest that relatively low abundance bacteria in the vaginal community could play important roles in the progression of inflammation ultimately leading to BV and vaginal symptoms. Dialister micraerophilus was only identified at the species level in the metatranscriptome, but the genus Dialister was identified in the metagenome as well. Given the similarity in the estimated effects of D. micraerophilus and the genus Dialister , it is possible that the taxa identified in the metagenome was composed primarily of D. micraerophilus . Both D. micraerophilus and Dialister were consistently and credibly associated with positive results in the individual elements of Amsel’s criteria and with an increased Nugent score. Dialister micraerophilus is a known BV-associated opportunistic pathobiont of the female reproductive tract 46 that has previously been implicated in increased risk of recurrent BV among people living with HIV 47 and pregnancy-associated vaginitis 48 . Our results add to these findings that Dialister could be important in determining the clinical manifestations of Amsel or Nugent-diagnosed BV and should be considered in the development of future molecular diagnostics and BV treatments. Surprisingly, of the three species credibly associated with vaginal discharge in the metagenomic data, two were Lactobacillus species: L. gasseri and L. jensenii. Prevotella bivia in the metagenome was also significantly associated with the presence of discharge. The Prevotella genus, including P. bivia has previously been shown to degrade mucins and has been associated with disruption of the vaginal epithelial glycocalyx 49 . Those previous findings combined with our results indicate that in context of BV, G. vaginalis may be the driving force of the condition as a binary state, but additional species play may play key mechanistic roles in the manifestations of singular Amsel’s criteria, in this case P. bivia and discharge. These findings are consistent with recent suggestions that the role Gardnerella plays in BV may be primarily colonization, and other species are opportunistic niche occupants that drive symptoms 50 , 51 .The associations between species and Nugent score were consistent with previous findings demonstrating that low Nugent scores have a higher association with Lactobacillus and high Nugent scores are associated with higher abundance of G. vaginalis and other aerobic species 15 , 38 . Notably, G. vaginalis strain B was also consistently associated with higher probabilities for individual Amsel’s criteria. This novel finding is, to our knowledge, the first indication that G. vaginalis has a strain level influence on BV. These data should help inform future studies on BV treatment and diagnostics 52 – 54 and suggest an opportunity to identify differences between strain B and other G. vaginalis strains that could help further illuminate BV symptomatology. This tend surfaced despite the two primary limitations of this study, a modest sample size, and the comparison of sequencing methods from both vaginal and cervical samples. However, while the Bayesian modeling results were compared across sample types, they were analyzed individually. There were differences between the cervical transcriptome and the vaginal metagenome, but the most interesting finding was apparent in both: the strain-specific effect of G. vaginalis. These results are particularly relevant in light of new data from probiotic and live-biotherapeutic BV interventions 55 – 57 . Those clinical trials found a benefit of seeding the vaginal microbiome with Lactobacillus crispatus strains, which improved symptoms and microbiome stability 55 – 57 . However, the effect appears largely reliant on continuous use of the probiotics, which may indicate that remnant populations of non-optimal microbes persist and opportunists. One trial of LACTIN-V also reported data showing that mild discharge continued in the treatment group at rates higher than the in the placebo group, while vaginal odor diverged 57 , a finding that aligns with our own, and also suggest a complex relationship between vaginal symptoms and the microbiome. In summary, our results highlight that sequencing methodology and granularity have a substantial effect on CST assignments and therefore the implied community structure and molecular BV status. Furthermore, we highlight that G. vaginalis has a strain dependent association with the individual elements of Amsel’s criteria and Nugent score. In addition to G. vaginalis , Dialister , P. micra , and P. bivia were consistently associated with BV-results for the individual elements of Amsel’s criteria. Our work particularly builds on a previous study which analyzed the symptomatology of different molecular variations of BV 13 . Srinivasan et al. similarly found that assessing the presence-absence patterns of bacteria identified via metataxonomics with Amsel’s criteria showed that some bacteria are related to a specific symptom rather than the entire diagnostic criteria of BV 13 . Like our study, they saw that P. bivia, G. vaginalis , and Dialister have important but distinct associations with pH, the whiff test, and clue cells. Our study differs and builds on those results by identifying strain level variation within G. vaginalis , and the associations we saw between Parvimonas and Amsel’s BV criteria. Srinivasan et al. also found that Parvimonas was uniquely associated with discharge 13 , while our results suggest that Parvimonas is negatively associated with discharge, but positively associated with the other three Amsel’s criteria. Taken together, our results highlight the complexity of BV and provide targets for continued exploration of the differences between BV types. Our findings are also consistent with the idea that no diagnostic consistently captures the nuances of BV or fully describes the underlying differences that lead to unacceptably high rates of recurrence.

Bacterial vaginosis (BV) affects 19–30% of women in the US and up to 60% in vulnerable populations with limited access to healthcare including lower-resourced women in Sub-Saharan Africa 1 – 4 . BV substantially impacts quality of life 5 , 6 with symptoms including vaginal discomfort, abnormal vaginal discharge, and odor changes 7 . BV is associated with pre-term birth 8 and increased risk of sexually transmitted infections 1 , including HIV 9 – 11 . Despite the prevalence of BV and the severity of the associated morbidities, BV remains understudied and without a universally accepted etiology. Bacterial vaginosis was differentiated from non-specific vaginitis in the 1950s, most notably in the work of Gardner and Dukes that first associated Gardnerella vaginalis (at the time named haemophilus vaginalis ) with the clinical manifestations now typically diagnosed as BV 12 . They claimed that G. vaginalis fulfilled Kochs postulates and should be considered the infectious cause of BV 12 . However, subsequent research has demonstrated that BV involves more complex changes and interactions of the vaginal microbiota 1 , 2 . One of the central dynamics in BV is the relative abundance of two taxa, Lactobacillus spp. and Gardnerella vaginalis . When Lactobacillus predominates the vaginal microbiota, it has a protective effect by producing lactic acid and maintaining a pH below 4.5 in the vaginal microenvironment, thereby restricting the growth of non-optimal microbes. In contrast, the presence of G. vaginalis is associated with a higher vaginal pH and increased polymicrobial vaginal microbiota. A relative decrease in vaginal Lactobacillus and accompanying increase in G. vaginalis is often associated with vaginal symptoms characteristic of BV 13 . Those observations informed the development of two clinical diagnostic tools for BV: Amsel’s criteria 14 and Nugent scoring 15 . Women are diagnosed with BV by Amsel’s criteria if clinical examination identifies three of four vaginal characteristics: abnormal vaginal discharge, a vaginal pH above 4.5, the production of amine odors when vaginal discharge is exposed to 10% potassium hydroxide (often termed the whiff test), and the presence of more than 20% clue cells in vaginal swab preparations evaluated under high power microscopy 14 . Nugent scoring relies on the microscopic evaluation of gram-stained vaginal smears and the presence of specific microbial morphotypes: large Gram-positive rods, small Gram-negative or variable rods, and Gram-negative curved rods. Low numbers of gram-positive rods and high numbers of Gram-negative morphotypes contribute to a high Nugent score and positive BV diagnosis. Applying sequencing approaches to characterize the composition and function of the vaginal microbiome has added another layer to our understanding of BV and vaginal microbiota. 16S rRNA gene metataxonomics (metataxonomics hereafter) studies of vaginal microbiota has demonstrated that there are five broad types of vaginal microbiota referred to as community state types (CSTs) 13 , 16 – 19 that differ in their microbial composition. Studies of vaginal microbiota composition that employ sequencing techniques have shown that most 11 , 19 women have vaginal microbiota predominated by one Lactobacillus species. These communities fall into CST I characteristically predominated by L. crispatus , CST II predominated by L. gasseri , CST III predominated by L. iners, and CST V predominated by L. jensenii . In contrast, some women have polymicrobial vaginal microbiome compositions that often include Gardnerella, Prevotella, Fannyhessea (previously Atopobium ), and Dialister among others 20 , 21 ; this community is designated CST IV. By measuring vaginal microbiota composition and structure, CST profiling represents a more precise and granular characterization and CST IV has been proposed as a molecular diagnosis for BV since it largely corresponds with Amsel’s and Nugent-defined BV 11 , 19 . Amsel’s and Nugent assessments measure different characteristics of the vaginal microenvironment, but both were heavily influenced by an understanding, now somewhat outdated, that BV is an infectious condition that can be cured by eliminating an overgrowth of anaerobic bacteria 14 , 15 . Considerable research provides evidence that BV is not a simple infection 1 , 11 , 21 , 22 , suggesting instead that BV is a collection of similar shifts from an optimal vaginal microbiota dominated by Lactobacillus. The continued lack of consensus on the definition and diagnostic criteria for BV has profound real-world clinical implications. Inconsistent diagnostic methods can contribute to misdiagnosis, treatment failures, and high recurrence rates. Amsel’s criteria rely on clinical symptoms, meaning asymptomatic individuals with a microbiome composition consistent with BV may go undiagnosed and untreated, while symptomatic individuals without a corresponding microbial profile may receive unnecessary antibiotic therapy. The clinical treatments for BV, with the antibiotics metronidazole and clindamycin being the most common 23 , are quite effective for initial clinical symptom resolution, but BV has a recurrence rate of up to 60% within a year of treatment 24 – 26 . While the clinical symptoms of BV are of primary concern to individuals who seek treatment, those symptoms do not always correspond to observed perturbations of the vaginal microbiome 21 . Those relationships are complicated further by the subjectivity of some symptom reporting. Nonetheless, research suggests associations between changes in the vaginal microbiota and numerous other conditions including endometriosis 27 – 29 , pelvic inflammatory disease 28 , STI risk 9 , 10 , 21 , and pregnancy complications 8 . Many of those associations appear independent of symptomatic BV, suggesting that even transient asymptomatic changes in the vaginal microbiome are of profound importance to advance women’s and neonatal health. This highlights that clinical symptom-based and microbial-based diagnoses should align if they describe the same discrete condition. Nonetheless, the etiology of BV, and even its definition, remains ambiguous. Data demonstrate that BV is not a simple condition caused by the presence of a single etiological agent and the complex transition to non-optimal vaginal microbiome states suggests instead that BV is a collection of several distinct conditions induced by host-microbe interactions, as others have noted 30 , 31 . A more granular understanding of the relationships between diagnostic methods, as well as between individual clinical, molecular, and microbiome indicators, is needed to generate translational results that improve clinical care and inform the development of interventions for different types of BV. Human microbiome sequencing aims to comprehensively understand the true composition and function of our microbiota, but myriad molecular and computational challenges make perfectly capturing the vaginal microbiome difficult if not impossible. Those challenges notwithstanding, even imperfect estimates of the microbiome are incredibly valuable, and each type of sequencing recapitulates the microbiome in valuable ways. Metataxonomics provides robust and economical identification of taxa but lacks comprehensive and robust species resolution and is insensitive to bacterial activity and viability. Metagenomics improves species resolution over metataxonomics and is similarly based recovery of genomic DNA. Metatranscriptomics has similar species resolution to metagenomics while also being sensitive to microbial activity and viability since it is based on mRNA. By relying on DNA, both metataxonomics and metagenomics provide an estimate of community composition in terms of genome number, regardless of bacterial activity or viability, whereas metatranscriptomics estimates microbial contribution to the active microbiome. Each method is reflective of the true microbiome, but reflective in different ways 20 . Continued research is needed to determine how different sequencing methods affect microbiome interpretation in the context of BV and its clinical manifestations. Here, we help fill those knowledge gaps by comprehensively examining differences between clinical and microbiome-based diagnostic techniques and identifying highly granular associations between individual vaginal and cervical bacteria and vaginal health data.