In vivo microbiome and associated immune markers: New insights into the pathogenesis of vaginal dysbiosis

The microbiota fulfils a key role in the training and function of the immune system, which contributes to the symbiosis between the host and complex microbial communities. In this study, we characterized the interplay between vaginal bacteria and local immune mediators during dysbiosis in selected women of reproductive age who were grouped according to Nugent’s criteria. The abundance of Gardnerella vaginalis and Bifidobacterium breve was increased in the intermediate dysbiotic status, while the presence of a plethora of non-resident bacteria characterized the group with overt vaginosis. In response to these increases, the anti-inflammatory IL1ra and pro-inflammatory IL2 increased, while the embryo trophic factors FGFβ and GMCSF decreased compared to the healthy milieu. A specific pattern, including IL1α, IL1β, IL8, MIG, MIP1α and RANTES, distinguished the intermediate group from the vaginosis group, while IL5 and IL13, which are secreted by Th2 cells, were significantly associated with the perturbation of the commensals Lactobacilli, Gardnerella and Ureaplasma. Summarizing, we postulate that although the dysbiotic condition triggers a pro-inflammatory process, the presence of a steady state level of Th2 may influence clinical manifestations. These results raise clinically relevant questions regarding the use of vaginal immunological markers as efficacious tools to monitor microbial alterations.

We profiled the vaginal microbiome and the local immune fluctuation in selected women of reproductive age grouped according to the Nugent's criteria with the goal of investigating the role of local immune mediators during the vaginal dysbiotic process.
In women with an intermediate Nugent score, we observed that the overgrowth of Gardnerella vaginalis was balanced by a high colonization of Bifidobacterium breve, which, overall, was able to guarantee a healthy vaginal equilibrium. Conversely, in women with BV, we observed a massive increase of non-resident vaginal species that paralleled a low presence of G. vaginalis.
The local immune response of the host was specifically stimulated by the vaginal dysbiosis, with a simultaneous increase of the anti-inflammatory IL1ra and the pro-inflammatory IL2 and a decrease of the embryo trophic factors FGFβ and GMCSF. Moreover, a specific pattern of mediators, including IL1α, IL1β, IL8, MIG, MIP1α and RANTES, distinguished the grade of vaginal dysbiosis observed in women with an intermediate Nugent score and in women with overt vaginosis.
Among the tested immune mediators, a significant association of the modulation of specific resident bacterial communities was highlighted for IL5 and IL13. The increase in the concentration of these two anti-inflammatory proteins was accompanied by a depletion of Lactobacilli ssp., G. vaginalis and Ureaplasma spp. The concomitant increase of the anti-inflammatory IL5 and IL13, secreted by Th2 cells, and the establishment of a vaginal dysmicrobism suggests a role for the Th2 activation in counteracting the Th1 response and in counteracting the presence of clinical symptoms.

Results
Characterization of women with vaginal dysbiosis. To characterize the vaginal dysbiotic process, we collected 62 cervical-vaginal samples from women who were diagnosed according to Nugent's criteria, including 30 women with a Nugent score of 0-3 (Healthy), 15 women with Nugent score of 4-6 (Intermediate) and 17 women with Nugent score of 7-10 (Vaginosis). The sequencing of the vaginal microbial community yielded 3,496,234 high-quality reads, excluding samples with less than 10,000 reads, with an average of 53,788 reads/ sample (range 10,162). We identified 15,769 operational taxonomic units (OTUs) across all the samples.
Bacterial diversity and comparisons of alpha diversity between cohorts. Out of the five distinct alpha diversity metrics (within-sample diversity), the Simpson, Shannon, observed species and PD_whole_tree metrics significantly characterized the dysbiotic vaginal microbiome ( Table 1) from that of the healthy microbiome, while Chao1 did not exhibit significant variations. Overall, across the five metrics used, the intermediate and vaginosis groups showed higher alpha diversities than the healthy group.
When accounting for the bacterial communities in the cohorts (Fig. 1 Among the other microorganisms detected, the intermediate group showed an increase of Gardnerella vaginalis (Actinobacteria) and Ureaplasma parvum (Tenericutes) compared to the vaginosis and healthy groups. Bifidobacterium breve was only detected in the intermediate group.
The vaginosis cohort exhibited a vast change in specific species, which are shown in Table 2.
We assessed which OTUs were discriminative between cohorts, and significant differences of bacterial species were not observed according to the FDR p value.
We next compared the overall microbial diversity between groups using the unweighted and weighted UniFrac distance matrices. The results of the β-diversity (between sample diversity comparison) were visualized by a Principal Coordinates Analysis (PCoA). As expected, the bacterial composition of these vaginal communities clustered according to the clinical grouping. Indeed, the analysis of the β-diversity highlights three distinct microbial clusters in the unweighted UniFrac PCoA (Fig. 2). The same clustering is partially confirmed by the weighted UniFrac PCoA (Fig. 3), although it is graphically less evident.
A one-way Analysis of Similarity (ANOSIM) statistical test was applied to the UniFrac distance matrices to test for significant differences according to the clinical grouping. ANOSIM attributed a significant difference to the grouping, shown by the weighted UniFrac (p = 0.001) and unweighted UniFrac (p = 0.001) measures, although the ability of the clinical grouping to explain the differences between cohorts is small (R = 0.31 and 0.21, respectively).

Correlations between microbiome and immune soluble factors.
Forty-eight soluble immune mediators, including cytokines, chemokines and growth factors, were measured in the vaginal samples to address the local immune response to microbial changes, the results of which are shown in Table 3. We observed that an increase of the anti-inflammatory IL1ra and the pro-inflammatory IL2, and a decrease of the embryo trophic factors FGFβ and GMCSF, distinguished the dysbiotic from the eubiotic microbiomes. A specific pattern of soluble immune mediators distinguished the grade of vaginal dysmicrobism. Specifically, a significant decrease of the proteins IL3, LIF, SCF, TNFβ, MCP1 and IL17 was associated with the intermediate group compared to the healthy group (Table 3). A significant increase of IL18, MIF and TNFα characterized the vaginosis group compared to the healthy group (Table 3). In addition, a panel of cytokines, including IL1α, IL1β, IL8, MIG, MIP1α and RANTES, exclusively differentiated the vaginal inflammatory process of the women who were diagnosed with an intermediate dysmicrobism from those with vaginosis (Table 3).
To evaluate the possible relationship between the local immune response and the bacterial composition, the concentrations of each immune factors were included as variables in the BIO-ENV rank-correlation procedure, where single variables or combinations of variables are selected based on how they best explain differences among samples.  Based on the weighted UniFrac distance matrix, several immune factors showed a high correlation with the different microbial patterns observed among samples. In the dysbiotic microbiota, many of the significantly upregulated cytokines, such as IL18, IL2, IL1ra, MIF, RANTES, TNFα, MIP1α and IL8, were confirmed to correlate with the alteration of the vaginal composition. Other immune factors, such as IL5, IL13, IL6, IL15, IL9, GROα, MIP1β and IFNγ, exhibited a correlation with the altered vaginal milieu, although their amount did not significantly change between cohorts.
Based on the unweighted UniFrac distance matrix, the soluble factors that correlated with the microbiome composition were MIF, IL18, MIP1β, TNFα and IL1ra (Table 4).  Next, we explored whether any individual bacterium correlated with an increase or decrease of specific immune factors. The associations were performed after rarefying the otu_table_biom (depth 10,000 sequences) and using the observation_metadata_correlation.py script with the Fisher_z_transform.

Discussion
The use of high-throughput sequencing techniques shed new light on the high variability and complexity of the vaginal microbiome. Our in vivo results agree with recent studies that have shown that Nugent score, based on   the quantification of only three bacterial morphotypes, is not sufficiently informative to clearly define a vaginal dysbiotic status 23 . Consistent with this issue is the observation of a high colonization by Bifidobacterium breve, neglected by Nugent Score, in the vaginal microbiome of women with an intermediate dysmicrobism. B. breve is able to counteract the suboptimal colonization by the lactate-producing Lactobacilli spp. 24 , guaranteeing a healthy vaginal equilibrium through the same mechanism exploited by Lactobacilli, that is, the production of lactic acid. It is noteworthy that, by neglecting the presence of all lactate producing bacteria and the identification of Lactobacilli species, the Nugent score criteria led to the overestimation of the dysbiotic status and to an unnecessary therapeutic intervention being suggested.
In addition, the Nugent score neglects the dynamic of interspecies communication between specific species of Lactobacilli and opportunistic pathogens. In the group of healthy women, L. crispatus 25  The novelty of our study consists of the in vivo demonstration that specific vaginal immune profiles significantly correlate with the microbial composition and clinical manifestation. This finding may help in the diagnosis of BV, predict the recurrence of vaginal dysbiosis and help follow the recovery after treatment. We observed that a subset of pro-inflammatory mediators, typically involved in the chronic inflammation process, synergized with a subset of cytokines that are involved in the switch toward the Th2 immune response. This network likely supports a dysbiotic condition, regardless of the host pro-inflammatory Th1 response, which is normally effective at restoring the eubiosis state.
More precisely, the women with an intermediate Nugent Score showed a significant decrease in several pro-inflammatory cytokines, including IL3 28 and IL17 13 . This event, together with the poor T cell-stimulating activity of G. vaginalis 29 , which was predominant in the intermediate group, may explain the absence or weakness of clinical symptoms 11 . Moreover, B. breve, which was uniquely identified in the intermediate group, can induce Th1 polarization of lymphocytes 30 , supporting their own growth by the production of IL2 31 and determining the resistance toward the infections. This mechanism may support a biotic balance assuring a healthy vaginal milieu.
In women with vaginosis, a massive pro-inflammatory response was observed, supported by a significant increase of IL18 32 , a constituent of the inflammasome complex 32 . The increase in concentration of other components of the inflammasome, such as IL1β, and the increase of several chemokines, such as IL8, MIG, MIP1α and RANTES, which are involved in the recruitment of leukocytes at the site of inflammation 33 , strongly differentiated women with vaginosis from women with an intermediate Nugent Score. This observation confirms the results of previous in vitro studies that showed an immunological shift towards a Th1-dominated profile during episodes of bacterial dysmicrobism 34 .
A new observation that emerged from our study is the significant association between specific commensal microorganisms and the modulation of local cytokines and chemokines, some of which have never been described before. Specifically, IL5 and IL13, which are secreted by Th2 cells 35 and are known to inhibit both the cell-mediated immune response and several macrophage functions, were statistically associated with Lactobacilli, Gardnerella and Ureaplasma species.
Based on the slight increase of these two cytokines (Table 5) and the depletion of some commensal bacteria (Fig. 1) in women with dysmicrobism, we suggest that the Th2 response is maintained at a steady state level of activation alongside the host Th1 response. The increase of the two anti-inflammatory cytokines is expected to significantly rise in asymptomatic women with an altered microbiome, explaining the absence or the weakness of symptoms. Furthermore, the steady-state Th2 response, together with a blunted Th1 response, could lead to immunologic tolerance causing chronic recurrent vaginal dysbiosis 36 .
Since IL5 and IL13 are strongly influenced by the alteration of eubiotic conditions, regardless of the specific pathogens causing the alteration of vaginal milieu, they could be further studied as indirect markers of vaginal disorder.
Taken together, our findings provide new in vivo insights into how different commensal bacterial species preserve the vaginal symbiotic equilibrium through their interplay with specific local immune mediators. These results raise clinically relevant questions regarding the role of vaginal immunological markers as crucial tool of surveillance of microbial alteration.

Methods
Patients and samples. Sixty-two immunocompetent women who fulfilled the inclusion eligibility criteria were included in this study. All women were Caucasian, of reproductive age (32-40 years old), were not pregnant, had no current use of hormonal or barrier contraceptive products, vaginal douching, tobacco or alcohol abuse, were not hospitalized or had systemic use of medication for chronic diseases or antibiotics/probiotics (oral or  Table 5. Soluble immune factors related to changes of the microbial composition. Data are shown as the mean value (pg/mL) ± standard error of mean. Comparisons between groups were performed using a nonparametric Kruskal-Wallis test. The p values are adjusted for multiple comparisons.
topical) within the 6 months prior to sample collection, and had no intercourse in the day prior to sampling. Microbiological criteria excluded concomitantly sexually transmitted viral infections, including HSV I-II and HPV. We identified patients among women that attended the Gynaecology Service of the Institute for Mother and Child Health IRCCS Burlo Garofolo of Trieste, Italy, as outpatients.
Vaginal samples were collected 7 days before the first day of the menstrual period. Under speculum examination, samples were collected using a 200 mm polyethylene Cervex brush device 37 (Rovers Medical Devices B.V., The Netherlands) by a single gentle 360° rotation of the cytobrush at the cervical os and were suspended in 1.5 ml of TE buffer. Each sample was divided into 3 (500 µl) aliquots and stored at −80 °C.
Sample processing and Ion Torrent Sequencing. DNA extraction was carried out using the NucliSENS ® easyMAG ® system (BioMèrieux, Gorman, North Carolina, USA) with an elution volume of 50 µl.
All DNA samples were stored at −80 °C prior to further processing.
A real time EvaGreen PCR (EvaGreen ® dye, Fisher Molecular Biology, Waltham, USA) was performed with the degenerate primer 27FYM (5′-AGR GTT YGA TYM TGG CTC AG-3′) and the primer U534R, targeting the V1-V3 region (500 bp) to allow for the construction of rich libraries. A nested PCR, targeting the V3 region 38 , was performed with the primers B338F_P1-adaptor (B338F 5′-ACTCCTACGGGAGGCAGC-3′) and U534R_A_barcode (U534R 5′-ATTACCGCGGCTGCTGG-3′) in conjunction with the IonXpress Barcode Adapter to obtain the 200 bp V3 region template for sequencing analysis 39 . Negative controls, including a no template control, were processed with the clinical samples. The PCR reactions were performed using the Kapa 2 G HiFi Hotstart ready mix 2 × (Kapa Biosystems, Massachusetts, USA), which has robust amplification that is necessary for NGS sequencing 40,41  Data Analysis. QIIME 1.8.01 42 was used to process the sequence data. High quality (Q > 20) sequences were demultiplexed and filtered by quality using split_libraries_fastq.py with the default parameters, except for the length parameter (at least 150 bp). Operational taxonomic units (OTUs) were defined at 97% similarity and clustered against the Vaginal 16S rRNA gene Reference Database, which was constructed by Fettweis et al. 43 , using open-reference OTU picking 44 with a uclust clustering tool 45 . Prior to further analysis, singleton OTUs and samples with low sequencing depth were removed (less than 10,000 reads). Chao1, PD whole tree, Shannon, Observed species and Simpson reciprocal metrics were used to assess alpha diversity (within-sample diversity), while beta diversity (between sample diversity comparison) was assessed with weighted and unweighted UniFrac distance matrices 46,47 and presented with principal coordinates analysis (PCoA). The robustness of the identified clusters was investigated using jackknifing (randomly resampling sequences without replacement). Differences in community composition between cohorts were investigated using analysis of similarity (ANOSIM, 999 permutations) and Kruskal-Wallis test, implemented in QIIME.
Correlations between the immune mediators, for which normalized amounts were written in the mapping file, and both weighted and unweighted UniFrac distance matrices (generated from the beta_diversity.py script) were assessed by the compare_categories.py script (metric BIOENV) of QIIME, exploiting the Spearman rank-order coefficient. A Spearman coefficient close to 1 signifies a highly positive correlation between immune mediators and the distance matrices; a value close to 0 signifies no association; and a value approaching -1 signifies a highly negative correlation.
To survey the association between microbial identities and the increase or decrease of specific immune factors, the observation_metadata_correlation.py script (with the Fisher_z_transform p-value assignment) of QIIME was used. The otu_table_biom was rarefied to a depth of 10,000 sequences/sample and the normalized amounts of the immune markers were written in the mapping file.
The dataset was deposited in the SRA database (PRJNA361297).
Immune soluble factor quantification. The quantification of soluble immune factors was performed using a recently described platform that is based on a magnetic bead multiplex immunoassay (Luminex, Bio-Plex, BIO-RAD Laboratories, Milano, Italy), which simultaneously detects 48 analytes, including cytokines, chemokines and growth factors 48 . Briefly, 50 μL of biological samples and standards were added in duplicate into a 96 multiwell plate containing the analyte beads. After incubation for 30 minutes at room temperature and washing, the antibody-biotin reporter was added and incubated for 10 minutes with streptavidin-phycoerythrin. The concentrations of the cytokines were determined using the Bio-Plex array reader (Luminex, Austin, TX). The Bio-Plex Manager software optimized the standard curves automatically and returned the data as Median Fluorescence Intensity (MFI) and concentration (pg/mL). To normalize the results, the total protein concentrations of samples were determined using a Bradford assay (Sigma-Aldrich, St. Louis, MO). Next, all cytokine and chemokine concentrations were normalized to total protein in the sample and were expressed as pg of immune marker⁄mL of total protein. Stata (v. 13.1) and GraphPad Prism (v. 5) were used for statistical data analysis. The Kruskal-Wallis one-way analysis of variance was used for comparisons between groups. When a significant p-value was observed (p < 0.05), a multiple comparison test was used to determine which groups were different.

Accession Codes. Sequences have been deposited at NCBI under the accessions PRJNA361297.
Ethics approval and consent to participate. The protocol used in this study was approved by the Ethics Committee of the IRCCS Burlo Garofolo Institute, Trieste (RC 26/13). All participants provided written informed consent for the experiments involving human participants and gave permission to access their medical records. All experiments were performed in accordance with the Declaration of Helsinki.