Enlibio,Your biological research partner
Xinqidi Biotech Co.,Ltd,Wuhan,China 2008-2022
R&D 14th year

molBV reveals immune landscape of bacterial vaginosis and predicts human papillomavirus infection na

Issuing time:2022-01-18 10:07


Bacterial vaginosis (BV) is a highly prevalent condition that is associated with adverse health outcomes. It has been proposed that BV’s role as a pathogenic condition is mediated via bacteria-induced inflammation. However, the complex interplay between vaginal microbes and host immune factors has yet to be clearly elucidated. Here, we develop molBV, a 16 S rRNA gene amplicon-based classification pipeline that generates a molecular score and diagnoses BV with the same accuracy as the current gold standard method (i.e., Nugent score). Using 3 confirmatory cohorts we show that molBV is independent of the 16 S rRNA region and generalizable across populations. We use the score in a cohort without clinical BV states, but with measures of HPV infection history and immune markers, to reveal that BV-associated increases in the IL-1β/IP-10 cytokine ratio directly predicts clearance of incident high-risk HPV infection (HR = 1.86, 95% CI: 1.19-2.9). Furthermore, we identify an alternate inflammatory BV signature characterized by elevated TNF-α/MIP-1β ratio that is prospectively associated with progression of incident infections to CIN2 + (OR = 2.81, 95% CI: 1.62-5.42). Thus, BV is a heterogeneous condition that activates different arms of the immune response, which in turn are independent risk factors for HR-HPV clearance and progression. Clinical Trial registration number: The CVT trial has been registered under: NCT00128661.


Bacterial vaginosis (BV) is defined as vaginal dysbiosis with inflammation and accompanying symptoms including vaginal discharge1,2,3. According to Centers for Disease Control and Prevention (CDC) and NHANES studies, the prevalence of BV is 29.2% amongst reproductive-aged women living in the United States4,5. Globally this condition is estimated to have an economic burden of approximately $5-billion per year6,7.

In addition to its ubiquity, BV is a urogenital condition that has been associated with adverse reproductive health outcomes including infertility8, increased risk for pre-term birth9, and low birth weights10. Moreover, an active state of BV is associated with an elevated risk for transmission of a variety of sexually transmitted infections (STIs) ranging from bacterial pathogens such as Chlamydia11 and Mycoplasma12, to viral agents including HIV13,14 and human papillomavirus (HPV)15. There is increasing interest in understanding the relationship between the cervicovaginal microbiome (CVM) and HPV natural history and progression to cancer16,17,18,19,20. In fact, differences in the CVM might explain why some high-risk HPV (HR-HPV) infections resolve, while others persist and progress. Lastly, BV is also associated with non-reproductive health issues such as obesity21.

Clinical BV is primarily diagnosed using Amsel criteria22, which requires the presence of three out of four signs or symptoms: (1) homogeneous, thin, white discharge that smoothly coats the vaginal walls; (2) clue cells in a wet mount; (3) pH of vaginal fluid >4.5; and (4) a fishy odor from the vaginal discharge before or after addition of 10% KOH (i.e., whiff test). Although commonly used, this approach has been widely criticized for a considerable rate of misdiagnosis23. An alternative to Amsel’s criteria is the Nugent score that creates a composite score based on counts of key bacteria morphologically identified on a Gram stain (i.e., Lactobacillus, Gardnerella, and curved Gram-negative rods)24. Although this method is more sensitive than the Amsel criteria24,25,26, it has been shown to suffer from interobserver variability27 and its use has primarily been limited to research settings due to the amount of time, expertise, and costs required to perform the test28. The term, molecular BV has been introduced recently14 and there are various meanings depending on the system used for molecular detection and the correlation with clinical, bacteriologic, and/or microscopic BV2,28. It specifically refers to suboptimal states of the CVM that are usually associated with reduced levels of Lactobacillus as measured by molecular techniques.

Bacterial vaginosis has features of an inflammatory state and is associated with alterations of cervicovaginal cytokines29,30,31,32. A number of studies have reported the association of elevated IL-1β and BV31,32, whereas most immune markers associated with BV appear to differ across studies. It has been proposed that this variability may be due to small sample sizes, heterogeneity of study populations33, and/or differing microbial taxa within the CVM. It is important to identify the source of this variability since the pathogenic effects of BV appear to be associated with local inflammation14,32,34,35,36.

In this study, we describe a 16S rRNA gene amplicon sequencing-based algorithm, called molBV, that can reproducibly categorize BV using a Nugent-like 0-10 score across a variety of populations including those from the US and Africa. Using this molecular approach to identify BV, we report the association of a set of cervicovaginal cytokines with molBV-BV. In particular, we demonstrate that elevated levels of Lactobacillus iners may in part explain the detection of a BV-like inflammatory signature amongst molBV-BV negative women. Although there appears to be a predominant host immune response to molBV-BV, the CVM’s positive associations with alternative forms of inflammation are associated with specific microbial agents. We utilize these observations to explore risk factors for the rate of clearance and progression of oncogenic HPV. We provide evidence that an inflammatory cervical profile underlies the association of HR-HPV natural history with molBV-BV. Surprisingly, the alternative inflammatory pathway is associated with the progression of HR-HPV infections to neoplastic lesions. This study provides evidence of multiple host inflammatory pathways associated with the cervicovaginal microbiome that influence the outcome of cervicovaginal HPV infection and possibly other pathologic outcomes of bacterial vaginosis.


Developing a molecular bacterial vaginosis scoring system

Initially, 30 young women with and without symptoms of BV were recruited for evaluation of Nugent, Amsel and 16S amplicon sequencing. Three samples were inadequate for study leaving a training set of 57 participants (mean age = 21, range 15–25 years). Based on Amsel’s criteria, 22 were classified as BV-positive; whereas, Nugent score evaluations categorized 26 with BV, 8 as intermediate for BV, and 23 as inconsistent with BV (Supplementary Table 1).

The 16S rRNA gene V4 region was amplified from cervicovaginal samples of all 57 participants, as it has been shown to robustly detect bacterial species from the cervicovaginal region37. There was an average (SD) of 16,580 (487) 16S reads per sample. Fungal sequencing of an ITS1 region amplicon using recently validated primers38 resulted in an average (SD) of 16,290 (3486) ITS reads per sample. Following taxonomic assignments and clustering by Euclidean distances (Fig. 1A), samples formed two primary clusters that were either defined by a dominance of two major species from the genus Lactobacillus (n = 27) or a state of polymicrobialism (n = 30). There was a highly significant tendency of the BV-positive samples to sort to the polymicrobial clade and BV negative samples to the Lactobacillus clade based on either the Amsel (p < 0.001) or Nugent BV diagnosis (p < 0.001). Hierarchical clustering using fungal communities revealed two primary clades: one dominated by Candida albicans, the other with a dominance of Malassezia restricta(Supplementary Fig. 1). The fungal community clustering showed no significant association with binary BV diagnosis, although some clustering was observed for both the Nugent (p-value = 0.18) and Amsel BV-positive samples (p-value = 0.22).

Fig. 1: Microbial community features of bacterial vaginosis.
figure 1

Panel A shows a heatmap of the 15 most prevalent bacterial species that are indicated to the right of the heatmap. Each column represents a participant. Hierarchical clustering separates samples into two primary clades: one dominated by Lactobacillus and one with polymicrobialism. There is a significant tendency of the BV-positive cases to be found in the polymicrobial clade and BV-negative in the one dominated by Lactobacillus based on either the Amsel or Nugent diagnosis (see the “Methods” section) (p < 0.001 for both). Panel B shows the alpha diversity differences between BV diagnosed by either Amsel or Nugent criteria and the microbial communities based on the Chao1, Fisher, and Shannon diversity indices (all p < 0.001), as indicated at the right of the panels. Panels Cand D show beta diversity analyses using PCoA and the Jensen–Shannon diversity index for the Amsel BV (panel C) (R2 = 0.25, p < 0.001) and Nugent BV diagnosis (panel D) (R2 = 0.59, p < 0.001). Panel E shows the top 20 microbial markers (based on W-stat) for detecting BV using the “clean” BV status sample set (Amsel+/Nugent+ vs. Amsel−/Nugent−). The y-axis represents the ANCOM W-stat, while the x-axis represents the mean relative abundance difference between BV+ and BV− cases for each bacterial taxon. The size of the circles represents relative abundance. Source data are provided as a Source Data file.

All alpha diversity measures (Chao1, Fisher, and Shannon, all p < 0.001) were highly associated with both Amsel and Nugent outcomes of BV (Fig. 1B). Beta diversity analyses using Jensen–Shannon Divergence (JSD) distances subsampled to 10,000 reads revealed that both the Amsel and Nugent criteria for BV were also significantly associated with the vaginal microbiome R2 = 0.25, p < 0.001 (Fig. 1C) and R2 = 0.59, p < 0.001 (Fig. 1D), respectively. To identify specific taxa associated with BV using ANCOM, we focused on the samples with concordant results for BV by Amsel and Nugent criteria (Supplementary Table 1); 52 differentially abundant genera were identified (FDR < 0.05, Fig. 1E), with Lactobacillusbeing the dominant genus elevated in BV-negative women and a mixture of anaerobic Gram-negative bacteria such as Gardnerella elevated in BV-positive, as expected39,40,41. There were no significant associations of fungal alpha or beta measures or specific fungal taxa identified with BV states (Supplementary Fig. 2A–C).

We sought to define a single molecular score from the 16S rRNA gene amplicon next-generation sequencing (NGS) data that would maximize generalizability of such a metric. Thus, we limited the markers to those taxa present in > 80% of all samples at a relative abundance of ≥ 0.001% after subsampling to 10,000 reads. We identified 11 genera meeting these criteria including Lactobacillus, Prevotella, Gardnerella, Megasphaera, Parvimonas, Clostridium, Porphyromonas, Adlercreutzia, Dialister, Atopobium, and Sneathia. We then derived a molBV score using robust regression modeling as described in the “Methods” section and created an averaged score of microbial taxa ratios for each sample providing a score of 0–10 similar to a Nugent score output.

molBV prediction of Nugent BV using three independent clinical datasets

molBV was evaluated in publicly available datasets that included 16S rRNA gene amplicon NGS data and measures of bacterial vaginosis. One testing set contained 388 American women with available 16S data sequenced using the V1–V2 16S rRNA gene region (different from the V4 region used above) and clinical Nugent scores40. In addition, we identified two African populations, one collected in Cape Town (n = 90) and the other in Soweto (n = 78) that sequenced the V4 region of 16S and had Nugent measures of BV42. We ran the 16S amplicon NGS reads through the molBV pipeline to generate Nugent-like scores and observed a strong correlation between the clinical Nugent scores and the molBV scores in all three cohorts with ICC values between 0.71–0.81 (Supplementary Fig. 3A). We next assessed the molBV score as a discriminant tool for BV diagnostic categories similar to Nugent scores (BV-negative = 0–3 or BV-positive = 7–10)43. The molBV score showed high AUC values (0.88–0.98) in all three datasets and outperformed other measures of the microbiome such as alpha diversity measures Chao1 and Shannon and the relative abundance of Lactobacillus(Supplementary Fig. 3B-D). Thus, the molBV pipeline is a robust tool to convert 16S NGS data into BV categories independent of 16S amplicon region and population characteristics.

The inflammatory landscape of BV

Previous studies indicated that vaginal dysbiosis is associated with an innate immune response32,44. To further investigate the host immune landscape and bacterial vaginosis, we utilized the molBV tool to recapitulate categories of vaginal dysbiosis where other measures of BV were unavailable. We utilized 431 baseline samples from individual women participating in the placebo arm of the Costa Rica Vaccine Trial (CVT) that had 32 cytokine proteins (i.e., cytokines, chemokines and soluble receptors) quantitated from cervical secretions collected with a sponge (see Methods)45. Using three ordinal categories of BV derived from the molBV scores equivalent to Nugent BV negative (molBV 0–3, n = 179), intermediate (molBV 4–7, n = 70) and positive (molBV 7–10, n = 182), we identified 13 cytokines significantly associated with a trend across the three BV states (Fig. 2A, all markers q < 0.001). Cytokine levels were also tested with respect to age, smoking and HPV16 status and did not show any significant associations (Supplementary Table 2). In order to validate the use of ordinal BV categories, we performed additional sensitivity analyses using categorical BV states and found that the categorical models did not provide a better fit (Supplementary Fig. 4). The strongest positive association of ordinal molBV states was with IL-1β (unit increase OR = 1.73, 95% CI: 1.56–1.92), whereas IP-10 was inversely associated with BV (OR = 0.76, 95% CI: 0.68–0.85). Supplementary Table 3 shows the cytokines associated with molBV-BV in univariate analysis; 6 have been previously associated with BV and 7 additional cytokines are described in this report. Given that inflammation is a complex host response, we analyzed the interrelationships amongst the 32 cytokines to identify patterns of expression. A correlation network was constructed connecting those cytokines with a Pearson correlation >0.6 (see Fig. 2B). Since 18 of the markers show a strong correlation to at least one other marker, we sought to further improve the identification of differentially abundant inflammatory signals by using cytokine ratios that overcome some issues with compositional data and relative abundances46. Figure 2C shows a volcano plot indicating cytokine ratios (in red) with a q < 3.77*10−44 threshold. The ratios with the strongest affects were IL-1β/IP-10 and IP-10/TNF-α, shown at the far right and left, respectively. To further examine the relationships of cytokine ratios, we created a matrix showing the pairwise correlations associated with molBV-BV (Fig. 2D). Six highly correlated ratios share the IL-1β cytokine and 2 other highly correlated ratios share TNF-α (the 1st and 2nd highest ORs for a given cytokine, respectively, in the univariate analysis shown in Fig. 2A). Of the BV associated ratios, IL-1β/IP-10 had the strongest overall effect based on absolute odds ratio. Interestingly, IL-1β and IP-10 have been previously shown to be strongly associated with BV47, making this ratio very attractive for further consideration.

Fig. 2: Cytokines associated with molBV categorical states in the CVT dataset.
figure 2

Panel A forest plot showing the OR and 95% confidence intervals computed using a linear model between the cytokine quartiles and the three ordinal states of molBV (i.e., BV-negative, BV-intermediate, and BV-positive) derived using 431 biologically independent CVM samples. Only cytokines with an adjusted q-value < 0.001 are presented. ORs in panel A represent the odds of moving to the immediate next ordinal BV state. Panel B shows the correlation network between all cytokines with a Pearson correlation >0.6 using all baseline samples. Panel C is a volcano plot showing the ORs on the x-axis when comparing molecular BV-negative vs. BV-positive and the –log(q-value) on the y-axis. The ratios that had a q-value < 3.77*10-44 (i.e., −log(q-value)>100) are indicated and labeled in red. Panel D shows the pairwise Pearson correlation of the highly significant ratios presented in panel C (colored red). Ratios that had an OR < 1.0 were inverted for symmetry of data presentation. Clusters from the strongest molBV-associated cytokine ratios appear to fall into two primary groups; ones that include IL-1β and those with TNF-α. Panel Epresents a box and violin plot for the log(IL-1β/IP-10) ratio for BV-negative (colored in blue), BV-intermediate (colored in yellow), and BV-positive (colored in red) samples. Panel F shows the bacterial species identified by ANCOM predicting IL-1β/IP-10 inflammation (i.e., above the median) in women within the BV-negative group shown in the small figure below with a red border. Comparison is between samples above the median vs. below the median. Panel G shows the ORs and 95% confidence intervals for the top 10 ratio combinations of the 32-cytokines (based on adjusted q-value) when comparing BV-negative to BV-positive women that had IL-1β/IP-10 levels below the cohort median (n = 171 biologically independent samples). Source data are provided as a Source Data file.

Exploring the distribution of the IL-1β/IP-10 ratio across the molBV-BV ordinal states revealed that despite the clear and consistent association (linear trend p-value = 3 × 10−48) with BV (Fig. 2E), 24% (85/349) of women with molBV-BV did not have elevated IL-1β/IP-10 levels. Whereas, surprisingly 28% (88/309) of BV negative women, had elevated IL-1β/IP-10 levels (see left violin plot, Fig. 2E). We used ANCOM to explore potential causes of molBV-negative women having high IL-1β/IP-10 levels; this analysis identified elevated levels of L. iners and G. vaginalis. Whereas, those with elevated levels of L. crispatus had low levels of IL-1β/IP-10 (W-stat threshold>10, FDR < 0.05, Fig. 2F).

To determine whether women with molBV-BV who did not show an elevated IL-1β/IP-10 signature had an alternative form of BV-associated inflammation, we compared this group to BV-negative women with similarly low levels of IL-1β/IP-10 (i.e., below the cohort median, Fig. 2G). The analysis revealed that TNF-α/MIP-1β was significantly positively associated with molBV-BV (OR = 1.24, 95% CI: 1.19–1.29) in the absence of elevated IL-1β/IP-10.

Molecular BV, cervical inflammation, and the natural history of HPV

We previously reported that increased diversity of the cervicovaginal microbiome contributed to HPV natural history16. To evaluate the impact of bacterial vaginosis and directly test whether associated inflammation could be mediating the effect of BV on the natural history of HPV, we utilized the previously reported16 prospective cohort sub-study from the CVT trial. We utilized 16S NGS data from cervicovaginal DNA16 to calculate the molBV scores across two study visits (307/431 baseline participants had sequenced 16SV4 and cytokine data). Women who had sustained low levels of molBV vs. those that had sustained high molBV scores were more likely to clear HR-HPV over time (Fig. 3A, p = 0.02). Briefly, sustained levels refer to women that had molBV levels above or below the cohort median for both of the measured visits (203/307 were included). Similarly, sustained high-levels of BV-associated inflammation vs. low, as determined by IL-1β/IP-10, were associated with lower rates of HR-HPV clearance (Fig. 3B, p = 004). Sustained levels of this measure were also defined using stratification by the cohort median and agreement of the measure (above or below) at the two measured visits (183/307 were included). For detailed definitions of sustained-levels of molBV and IL-1β/IP-10 inflammation see the subsection “HPV natural history exposure/outcome definitions” in the “Methods” section of the manuscript.

Fig. 3: Molecular BV, inflammation, and HR-HPV clearance.
figure 3

Panel A shows the Kaplan–Meier curves for HR-HPV clearance colored by sustained BV status (i.e., having a molBV value above (red) or below (blue) the cohort median for both measured visits) with the unadjusted p-value presented in the bottom left corner of the plot and number of women at risk at each time point presented in the accompanying table directly below the plot. Panel B shows the Kaplan–Meier curves comparing women with either sustained high (red) or low (blue) IL-1β/IP-10 inflammation markers across the two analyzed visits. Sustained high and low IL-1β/IP-10 refers to women that had IL-1β/IP-10 levels above or below the cohort median for both measured visits, respectively, with the unadjusted p-value presented in the bottom left corner of the plot and number of women at risk at each time point presented in the accompanying table directly below the plot. Time is shown in days at the bottom of the figures and tables. Source data are provided as a Source Data file.

To determine whether molBV-BV and/or sustained BV-associated inflammation markers acted as independent risk factors for HR-HPV clearance, we used cox-proportional hazard models after covariate adjustment. Table 1 presents the effects of molBV and IL-1β/IP-10 levels adjusted for age, smoking status and HPV16. Model 1 considered the dichotomized molBV states and found that compared to having sustained low levels of molBV (reference), the transition from high (V1) to low (V2) was a significant protective factor against clearing a HR-HPV infection (HR = 0.55, 95% CI: 0.31–0.97). Model 2 considered the effect of both IL-1β/IP-10 and molBV states. In this cytokine-adjusted analysis, molBV levels were not associated with HR-HPV clearance, suggesting that IL-1β/IP-10 was an independent driver of HPV clearance (HR = 1.87, 95% CI: 1.08–3.20). In this context, molBV states were not significant with the exception of the marginal signal from the group that became low at visit 2 (HR = 0.38, 95% CI: 0.15–1.00). Given the strong correlation between molBV and IL-1β/IP-10 levels (Fig. 2E), we made an additional parsimonious model (Model 3) to more accurately measure the effect of sustained low IL-1β/IP-10 levels on increasing the likelihood of HR-HPV clearance (HR = 1.86, 95% CI: 1.19–2.90). This model did not change the hazard ratio of sustained low IL-1β/IP-10 with HR-HPV clearance, but it did significantly reduce the p-value supporting the analyses that it is the true driver of HR-HPV clearance.

Table 1 Risk factors for time to clearance of incident HR-HPV infections.

We next evaluated the association of BV and HR-HPV progression to CIN2+16. Briefly, the original study considered the binary outcome of persistent HR-HPV progressing to CIN2+ (diagnosed ~2 years after the second visit sample16) vs. HR-HPV infection clearance with CVM components serving as predictors. Using values from V2, we first tested whether molecular BV was associated with progression to CIN2+ using a generalized linear model (Table 2, Model 1). The results indicated that a continuous molBV score was prospectively associated with increased risk for CIN2+ progression with the odds of progression increasing by 1.24 (95% CI: 1.02–1.55) per unit increase of molBV. However, the dominant BV IL-1β/IP-10 signature was not significant when molBV was also in the model (Table 2, Model 2, OR = 1.15, 95% CI: 0.92–1.46). Remarkably, the alternative BV-associated inflammation signature represented by TNF-α/MIP-1β levels was a predictor of progression when either molBV (Table 2, Model 3, OR = 2.71, 95% CI: 1.46–5.61) or molBV and IL-1β/IP-10 were in the model (Table 2, Model 4, OR = 2.65, 95% CI: 1.39–5.68). Given the association of TNF-α/MIP-1β with molBV (Fig. 2C) and the modest correlation with IL-1β/IP-10 (Fig. 2D), an additional parsimonious model was constructed to more accurately measure the effect of TNF-α/MIP-1β with progression to CIN2+ without molBV or IL-1B/IP-10 (Table 2, Model 5, OR = 2.81, 95% CI: 1.65–5.42).

Table 2 Risk factors for progression to CIN2+.


In this study, we develop a relatively simple means to characterize cervicovaginal samples for BV using 16S rRNA gene amplicon next-generation sequencing (NGS) data. It is well known that features of the CVM are strongly associated with BV39,40,48,49,50,51,52, and our study takes this relationship further and provides a quantitative score of 0–10, equivalent to the Nugent score. This method is particularly useful for high-throughput analyses to determine a “molecular” Nugent-like score (molBV-BV) in women that might not have been evaluated for BV, but have an available cervicovaginal sample. Development of this method used a stringent diagnosis of BV including subjects concurrent for BV by both Nugent and Amsel criteria. The method was validated in three additional cohorts suggesting the generalizability of this particular molecular approach to generate a Nugent-like score. Although it is known that there is substantial variation between the cervicovaginal microbiome between African women and women with European ancestry53, use of the molBV algorithm in two African populations revealed high diagnostic AUCs for BV of 0.97 and 0.88 for Soweto and Cape Town sets, respectively. Moreover, the molBV diagnostic was also used to evaluate the local host inflammatory response in a population without available BV measures (i.e., Nugent or Amsel). We identified a total of 13 different cytokines associated with molBV-BV, of which 7, had not been previously reported (Supplementary Table 3 and ref. 32). In addition, we utilize the molBV score and cytokine data to demonstrate the contribution of each to the natural history of HR-HPV infections using a prospective study design. Another interesting feature of this study was to demonstrate a prospective association of TNF-α with CIN2+ (previously reported in a cross-sectional study by Łaniewski et al. 54).

Cytokine ratios were used in order to better address the interrelated data structure and we observed that increasing IL-1β/IP-10 ratios were strongly correlated with increasing molecular BV scores (linear trend p-value = 3.67 × 10−48). IL-1β/IP-10 ratio was previously postulated to be a relevant signature for BV and the identification of women at higher risk of STI transmission47.

Despite a strong correlation between IL-1β/IP-10 and molBV states, there were still 24–28% of women that had an elevated ratio despite a molBV score of 0–3. Upon further analyses, elevated L. iners and G. vaginalis were identified in this enigmatic group. This finding is of interest since dominance of Lactobacillus species in the CVM is typically associated with vaginal health40,49. In the CVT cohort, L. crispatus was inversely associated with IL-1β/IP-10 inflammation consistent with previous reports40. It is of interest that women in the discovery set within a Lactobacillus-dominated clade and having clinical features of BV also had elevated levels of L. iners (Fig. 1A). It is possible that the association of L. iners with a BV-associated inflammatory state is due to strain-level variation since the biology of L. iners is perplexing55,56,57. This particular species of Lactobacillus differs from other members of the genus in many respects including genome characteristics that give it a more perplexing character as compared to other lactobacilli56. Furthermore, certain strains of L. iners appear to carry unique genes, such as those that encode for inerolysin58. France et al.51, postulated that these genes appear to have been horizontally transferred to L. iners from G. vaginalisand this allows certain members of this species to directly extract nutrients from host cells, which may explain how it can persist in a sub-optimal CVM and possibly induce vaginal inflammation. Deeper sequencing of the CVM will be required to validate this hypothesis. Alternatively, the community context of L. iners might influence its behavior and association with inflammation. It appears that specific bacteria could have context specific functionality depending on the observed molecular BV state at the time of sampling. Although the molBValgorithm was set to clinical parameters of BV, the current work suggests that additional stratification to identify women at risk for sub-clinical inflammation is necessary. This is especially important given the presented data showing the significance of these inflammatory shifts in the context of HR-HPV infections and possible implications for other diseases in which BV acts as a risk factor. Lastly, L. iners may reflect associations with other unmeasured determinants leading to BV56.

An alternative consideration that may play a role in the pathogenesis of BV, in addition to bacterial biomarkers is the importance of microbial biomass. Amplicon sequencing of the microbiome is compositional in nature and does not provide a direct means to establish exact microbial biomass46. This variable was previously shown to be relevant with the qPCR technique that allowed intermediate BV states to be resolved by quantifying bacterial loads of G. vaginalis and A. vaginae59. It would be interest to incorporate a quantitative technique such as the one developed by Morton et al. 46, to further expand molecular characterization of BV in future studies.

Another remarkable aspect of the cytokine analysis of this study is the heterogeneity of immune markers that were elevated across the strata of molecular BV states. It would be reasonable to predict that having an elevation of the predominant cytokine ration (i.e., IL-1β/IP-10) would yield similar distributions of other cytokines given the correlation levels observed in Fig. 2B and D. However, when comparing the molecular BV-positive to BV-negative women with low levels of IL-1β/IP-10, we observed that the BV-positive women had elevated levels of TNF-α/MIP-1β inflammatory cytokines. In fact, TNF-α/MIP-1β was the dominant signature elevated in molecular BV-positive women with low levels of IL-1β/IP-10. These results indicate that BV seems to be heterogeneous as to the exact type of sub-clinical host inflammatory response. Given the identifiable variance of inflammatory levels within individual strata of BV due to specific organisms such as L. iners, it is likely that specific taxa, or bacterial networks, might be associated with this observation. Further studies with deeper sequencing to identify the possible microbial genetic basis for these observations are needed.

There is increased interest in uncovering the relationship between the cervicovaginal microbiome and HPV natural history, since a number of studies show a correlation between vaginal microbial diversity, BV, and HPV clearance16,17,18,19,20. It is not currently known why certain high-risk infections clear while a small minority persists for years and eventually progresses to pre-cancer60. We utilized HR-HPV detection data within the Costa Rica vaccine trial cohort45 to evaluate cervical cytokine profiles, 16 S rRNA gene amplicon NGS data and the newly developed molecular BV states to interrogate possible mechanisms of HR-HPV infection clearance. Kaplan-Meier analyses revealed similar clearance patterns amongst women with sustained low vs. high IL-1β/IP-10 levels and low vs. high molBV scores (see Fig. 3), although the cytokine measure showed a stronger association with HR-HPV clearance. To better understand these relationships, we utilized Cox-proportional hazard modeling with adjustment for age, smoking status, and HPV16 infection. The association with molBV-BV was eliminated once we adjusted for an elevated IL-1β/IP-10 ratio, possibly indicating that a specific type of inflammation associated with BV was driving the relationship between HR-HPV persistence or clearance.

Another relevant HPV outcome is a progression of a persistent HR-HPV infection to pre-cancer (CIN2+)58. The CVM was previously reported to be predictive of this outcome16,20,36,61. However, when we tested the continuous IL-1β/IP-10 levels in a model, this cytokine signature did not add any additional information beyond molBV for HR-HPV progression. Surprisingly, when we tested the TNF-α/MIP-1β signature that is also associated with certain characteristics of molBV, we found that it was associated with CIN2+ progression even after adjustment for molBV as well as molBV and IL-1β/IP-10 (Table 2). The final analysis revealed that a single unit increase in the TNF-α/MIP-1β ratio was positively associated with an odds of 2.81 (95% CI: 1.65–5.42) of developing CIN2+ within 2-years of V2 sample collection. The data presented in the current report suggest that BV and the host response is a highly heterogeneous relationship and although BV is consistently associated with certain microbial shifts and overall community structure (e.g., higher alpha diversity), the host response can also be modified by the presence of specific taxa. These results may explain why certain studies do not see an association with CVM diversity (a surrogate for BV), but do see signals when analyzing specific bacteria62. Based on the observations reported in this study, these variations appear to have a substantial effect on the immune response, which in turn has an effect on HR-HPV clearance and progression to CIN2+.

The currently reported analyses have several weaknesses that should be taken into account when interpreting the data. All of the analyses of the CVM utilized 16S rRNA gene amplicon sequencing. This method limits the taxonomic resolution of bacteria and other organisms constituting the microbiome. A deeper exploration of BV using techniques such as shotgun metagenomics may provide a more thorough explanation as to why there is significant heterogeneity in the local host inflammatory response to BV amongst different women. Additionally, although our core analysis utilized compositionally aware approaches, rarefaction was used when calculating beta diversity, which may bias the magnitude of the PERMANOVA result. The relationships between molBV and the analyzed cytokines were considered in a linear context, other non-ordinal relationships might exist and are worthy of future investigation. Moreover, the developed molBV score was dependent on a relatively small set of bacteria and was based on “clean” BV diagnoses in which Amsel and Nugent’s tests agreed. This choice was made in order to facilitate the robustness of the measure. This appears to have been effective based on the associations seen within African populations but may present a limitation in populations with different structures of the CVM, especially ones that may occur in women where the two clinical scoring systems are discordant. Most of the women in the reported study are adolescents or young adults and we do not know if the analyses extend to women of all ages and geographic locations. Moreover, it is possible that other inflammatory signatures might exist in the cervicovaginal region that were not measured in the current report. Finally, in using the clinical Nugent score to guide our analyses we may have inadvertently missed important physiological phenomena of BV that are inherently diluted by this clinical score; future studies should utilize additional CVM reduction techniques such as the recently developed VALENCIA51, which produces community state types that are clinically agnostic.

Here we present a comprehensive molecular characterization of BV using 16S rRNA gene amplicon sequencing and a curated panel of cytokines. We demonstrate using multiple cohorts that 16S amplicon sequencing can be reliably used to diagnose BV employing the newly developed molBV score. We further demonstrate that this score was strongly correlated with a heterogeneous inflammatory landscape within the cervicovaginal region. Exploring these inflammatory markers further revealed a complex system of interactions between individual taxa, specific cytokines, and molecular BV states. In addition, we demonstrated that there is the potential for clinical relevance of the findings through the use of HR-HPV outcomes. We specifically show that different possible inflammatory states in BV are either associated with persistence of HR-HPV infection (i.e., IL-1β/IP-10) or the progression of infection to precancer (i.e., TNF-α/MIP-1β). In support of a role for TNF-α in the progression of HR-HPV infections, a recent report indicated that TNF-α was the main discriminatory biomarker associated with invasive cervical cancer63.

Whether the adverse health outcomes from BV are all based on the host inflammatory response remains to be rigorously evaluated. Deeper exploration of these associations is warranted using more robust techniques such as machine learning in order to further understand why certain women experience inflammation with BV, while others do not and how the host-microbiome relationships impact health. Lastly, the implication of BV inducing a local inflammatory response might imply signaling systemic inflammation, which was not assessed in this study. The role BV-induced inflammation might have in immune conditions more prevalent in women remains an interesting hypothesis that could have profound diagnostic and therapeutic ramifications.

Article classification: Biological abstract
Share to:
Add:Room A11-329, 1st Floor, No.1, SBI Venture Street, Optics Valley, East Lake
New Technology Development Zone, Wuhan, China.
Certificate NO.:U18Q28010569R0S