Effect of metronidazole on vaginal microbiota associated with asymptomatic bacterial vaginosis.

Vaginal dysbiosis-induced by an overgrowth of anaerobic bacteria is referred to as bacterial vaginosis (BV). The dysbiosis is associated with an increased risk for acquisition of sexually transmitted infections. Women with symptomatic BV are treated with oral metronidazole (MET), but its effectiveness remains to be elucidated. This study used whole-genome sequencing (WGS) to determine the changes in the microbiota among women treated with MET. WGS was conducted on DNA obtained from 20 vaginal swabs collected at four time points over 12 months from five randomly selected African American (AA) women. The baseline visit included all women who were diagnosed with asymptomatic BV and were untreated. All subjects were tested subsequently once every 2 months and received a course of MET for each BV episode during the 12 months. The BV status was classified according to Nugent scores (NSs) of vaginal smears. The microbial and resistome profiles were analysed along with the sociodemographic metadata. Despite treatment, none of the five participants reverted to normal vaginal flora - two were consistently positive for BV, and the rest experienced episodic cases of BV. WGS analyses showed Gardnerella spp. as the most abundant organism. After treatment with MET, there was an observed decline of Lactobacillus and Prevotella species. One participant had a healthy vaginal microbiota based on NS at one follow-up time point. Resistance genes including tetM and lscA were detected. Though limited in subjects, this study shows specific microbiota changes with treatment, presence of many resistant genes in their microbiota, and recurrence and persistence of BV despite MET treatment. Thus, MET may not be an effective treatment option for asymptomatic BV, and whole metagenome sequence would better inform the choice of antibiotics.

Introduction

Bacterial vaginosis (BV) is a polymicrobial condition with a dysbiosis of the vaginal microbiota where the levels of Lactobacilli decrease, while that of opportunistic, anaerobic bacteria increases [1]. BV is most common among women at reproductive age [2] and affects approximately 29 % of women in the United States [3]. It is associated with vaginal discharge, malodour, itching, and increased vaginal pH [3, 4]. However, about 50 % of cases of BV are asymptomatic [5]. BV has also been found to increase the risk of sexually transmitted infections (STIs) such as human immunodeficiency virus (HIV) [4, 6], herpes simplex virus type-2 [7], human papillomavirus (HPV) [8], [9], [9], and Trichomonas vaginalis [10]. Adverse reproductive sequelae, such as spontaneous abortion, preterm labour, and pelvic inflammatory disease, are associated with BV [2, 4, 6, 11]. Additionally, the aetiology of BV remains unknown, and the long-term effective treatment methods are yet to be determined [10, 12]. The longitudinal changes in the composition of the vaginal microbiota can allude to bacteria associated with healthy and unhealthy microbiota as well as the effect of treatment on their composition.

species dominate a healthy vaginal microbiota [2, 13–15]. In general, the genus has been found to lower the vaginal pH through the production of lactic acid, which promotes an ideal vaginal ecosystem [14, 16]. Among the species, and promote a healthy vaginal microbiota, whereas and do not [6, 14, 15]. In cases of BV, the vaginal microbiota is frequently dominated by spp., spp., spp., Atopobium vaginae, and [2, 6, 15]. The most typical antibiotics used to treat BV include metronidazole (MET) and clindamycin to inhibit bacterial growth [2, 10]. Though treatment often appears to be effective, there is an increased long-term recurrence (within 1 year) [2, 15]. Recurrent BV may probably be due to reinfection or increasing microbial resistance [15]. Demographic and behavioural characteristics further exacerbate the prevalence. For example, hormonal contraception has been shown to reduce the risk of BV [14]. An increased risk of BV can be attributed to race (African American), history of STIs, douching, and sexual behaviour [11].

BV diagnosis is based on either Amsel’s or Nugent tests [3]. Amsel’s criteria are widely used in clinical settings, whereas Nugent scores (NSs) are used extensively in research [10]. Amsel’s criteria are based on the presence of three of the four following signs: milky white discharge, presence of clue cells in wet mount microscopy, elevated vaginal pH, and whiff test [10, 17]. NSs range from 0 to 10 and are based on the relative abundance of bacteria in the Gram-staining of vaginal smears. A healthy vaginal microbiota is expected to have a NS of 0–3. Scores of 4–6 and 6–10 suggest intermediate flora and BV, respectively [18]. The vaginal microbiota can be further characterized by community state types (CSTs) [16]. There are five CSTs that describe the vaginal microbiota in terms of dominant bacteria; CSTs I, II, III and V are dominated by species, especially L. crispatus, L. gasseri, L. iners, and , respectively; and the CST IV is highly heterogenous and dominated by strictly anaerobic bacteria [16]. Most early vaginal microbiota studies commonly used 16S rRNA sequencing to describe the composition [6, 13, 19]. In this study, we examined the change in the composition of the vaginal microbiota of women with asymptomatic BV using whole-genome sequencing (WGS) and investigated the longitudinal change in the vaginal microbiota among women treated with MET. The changes in the relative abundance of bacterial species, alpha diversity, CST, and BV status across time points were examined. Furthermore, the correlation between NS and identified bacterial species were assessed. We report the exploratory findings of a pilot study with five women (20 samples) who were followed for 1 year at four time points.

Methods

Clinical samples

Vaginal swabs in this study came from a previously completed clinical trial funded by the National Institutes for Health, known as the BRAVO study [20]. The BRAVO study sought to investigate the effect of home screening and treatment of asymptomatic BV on incident and (NG/CT). This trial gathered data from women aged 15–25 years from five states in the US, namely Alabama, California, Maryland, North Carolina, and Pennsylvania. All women were diagnosed with asymptomatic BV at baseline using modified Amsel’s criteria. Women self-collected vaginal swabs at baseline and every 2 months over the next 12 months. Self-collected vaginal samples were shipped to a central laboratory in Pittsburg for testing. The vaginal swabs were tested for BV and NG/CT using NSs and the BD ProbeTec Amplified DNA Assay (Becton-Dickson, Inc. Sparks, MD), respectively.

Due to the pilot nature of the current study and limited funds, five women were randomly selected for WGS, which is known to provide a more detailed microbial profile than 16S sequencing. These five women were chosen from 124 subjects because they had the least possible confounding factors among the cohort. All five women chosen were African American, non-Hispanic, between the ages of 19–22, had never been pregnant before, did not have chlamydia or gonorrhoea, did not consume antibiotics during the survey, were sexually active, had prior BV exactly once, and stayed through the entire programme. The broader cohort of 124 subjects was not all African American; and some were Hispanic, many were pregnant during or prior to the study, and started with more severe BV at the start. Treatment was administered after baseline. Therefore, all baseline samples were untreated. Treatment was administered at each time point for which a woman tested positive for BV.

Epidemiological data

Baseline characteristics were assessed using a self-administered questionnaire in the parent study. Information on age, race (Black or African American/White or Caucasian/Hawaiian or Pacific Islander/Asian/Other), ethnicity (Hispanic/Latino, not Hispanic/Latino), and the highest level of education (eighth grade or less/Incomplete high school/High School graduate or GED/Incomplete college/Associate degree/College degree/Masters or advanced degree/Vocational certificate) were collected. Clinical and behavioural questions related to medical history explored prior antibiotic use (past 30 days), previous lifetime episodes of BV (never, once, 2–4 times, ≥5 times), prior lifetime treatment for BV (≤1 month ago, 1–6 months ago, 7–12 months, ≥12 months ago), prior pregnancy (yes/no), frequency of vaginal douching (daily/weekly/monthly/yearly/not in the past year/never), birth control methods used in the past year (pills, patch, Nuva-ring, condoms, a spermicide cream, Depo-Provera shot, intra-uterine device (IUD) and other), and the number of different sex partners and the types of sexual behaviour in the past year, (receptive oral sex, unprotected anal sex, unprotected vaginal sex, women who have sex with women, and new sex partners). The treatment arm assignment was also recorded.

BV assessment

BV was diagnosed using the NS [3]. A score of 0–3 is considered healthy vaginal microbiota, 4–6 intermediate flora, and 7–10 indicates BV. The vaginal microbiota can be further classified into five community state types (CSTs) [16]. The CSTs I, II, III, and V are dominated by , L. gasseri, L. iners, and L. jensenii, respectively, and are typically associated with a healthy vaginal ecosystem [16]. On the other hand, CST IV represents a diverse group with no individual dominant species but is associated with BV [16].

DNA preparation and whole genome sequencing

Vaginal swabs were processed and analysed with the ZymoBIOMICS Service (Zymo Research, Irvine, CA). Genomic DNA samples were profiled for shotgun metagenomic sequencing. Following the manufacturer’s protocol, sequencing libraries were prepared with KAPA HyperPlus (MA, US) with 100 ng DNA input. The libraries used internal 8 bp barcodes and TruSeq adapters. All libraries were quantified with TapeStation and pooled evenly. The final pool was quantified with quantitative PCR and sequenced with Illumina HiSeq. The operational taxonomic units were then counted, and the abundance of each taxon was normalized to a constant sum of one for each sample.

Bioinformatics and statistical analysis

Statistical tools

Data were analysed using the R statistical software package version 3.5.1 [21]. The plyr [22] and Hmisc [23] were used for matrix operations and data manipulation.

Microbial abundance analysis

The abundance profile for each subject was generated using in-house scripts. The resulting microbial profile of each sample was combined with clinical attributes of interest and labelled as treated/untreated correspondingly. Names of microbes were trimmed to focus only on the species level. Used the vegan [24] tool was to calculate the community diversity based on the Inverse Simpson Index [25].

Correlational analyses

Pearson and Spearman correlation coefficients [26] were calculated to find the taxa whose abundance values in the samples co-varied the most and the taxa that were most correlated with NSs, as described earlier [27].

Variance analysis

Variance analysis was performed on the sample microbial profiles. Principal component analysis (PCA), a dimensionality reduction tool, was used to assess the variance in the vaginal microbiota and to determine the taxa that contributed the most to the variance. It identifies the direction along which the samples display the most variance and also extracts the principal components that reduce the dimensionality of the data while preserving as much variance as possible. The balanced error rate of the Mahalanobis distance [28] was calculated for an increasing number of components. The mixomics R-package [29] was used for the PCA analyses.

Discriminant analysis

Discriminant analysis was performed on the sample microbial profiles to determine the taxa that contributed the most to the discrimination between the samples. The abundance of these taxa was then compared using Welch’s t-test [30], which was selected as there were unequal sample sizes, five samples untreated (baseline) versus 15 treated (follow-up). The Benjamini-Hochberg correction [31] was also used to decrease the false discovery rate. The Sparse Partial Least Squares Discriminant Analysis (sPLS-DA) [32] is similar to PCA, but it finds the direction along which the separation between the classes is maximized (discriminatory). The four classes for this analysis correspond to the sampling time points. Three-fold cross validation was repeated 100 times. The mixomics [29] was also used for the sPLS-DA analyses.

Analyses of ABR genes

The antimicrobial resistance profile for each subject was generated using scripts from PeTRi metagenomic pipeline [33]. First, accessed the list of 2786 antibiotic-resistance reference genes from the Comprehensive Antibiotic Resistance Database (CARD) version v3.1.1 [34]. Second, aligned metagenomic reads from each subject against the reference sequences using Bowtie2 [35]. Third, using the alignments computed ‘counts per million reads’ for each ABR gene, which served as an approximation of the metagenomic sample’s resistance ‘potential’. Since the CARD database does not contain nim genes associated with MET resistance, their presence was manually curated.

Visualization tools

The R libraries gplots [36], ggplot2 [37], and plotly [38] were used for the creation of all figures presented in this manuscript.

Results

Participant characteristics

The study sample included five AA women, who were less than 22 years of age (80%), did not have a college degree (80%) and lacked a history of pregnancy (100%). None of the women reported prior antibiotic use 30 days at the start of the study or a BV treatment history. Four of the five women reported douching in the past but did not douche regularly. These women actively used a form of birth control, most commonly condoms (75%), and had at least one sex partner in the past year. Most women had two or more oral sex partners in the past year (60%), engaged in vaginal sex without a condom (80%), and reported a new sex partner in the past year (80%). For this study, all five women received 1 week oral MET treatment after baseline assessment and every subsequent time point for which they tested positive for BV. Thus, at baseline, all vaginal samples were positive for BV and untreated. None of the participants had NG/CT.

Vaginal microbial profile

All the participants of this study had asymptomatic BV. Two participants (1 and 2) had BV (as per NSs) despite treatment for the entire year. Both participants reported having new and multiple sex partners as well as having received oral sex. In addition, Participant 2 also reported engaging in anal sex.

Fig. 1 describes the changes in the composition of vaginal microbiota over time for the five participants. Participants 1, 2, 3, and 4 were classified as having community state type (CST) V throughout the sampling period, even though they presented small variations in their microbial compositions. On the other hand, Participant 5 presented states III, I, V, and III in that order and was the only participant who transitioned to health-associated CSTs.

Fig. 1.

Vaginal microbial and resistome profile of five African American participants at four time points. Each column represents a specific time point for each participant beginning at baseline and followed by 4 month intervals, a total of four columns per participant ordered sequentially. Each column represents a specific time point for each participant beginning at baseline and followed by 4 month intervals, a total of four columns per participant ordered sequentially. The labels M00, M02, M08, and M12 correspond to baseline, and months four, eight, and twelve, respectively. The vaginal microbiota characterized by community state types I to IV are indicated in the x-axis. (a) Relative abundance of microbes. (b) The Nugent scores (NS) associated with the participant. The BV is indicated by red, intermediate vaginal flora by orange, and healthy vaginal microbiota by green. (c) Alpha diversity calculated using the Inverse Simpson Index. (d) Resistome profile is given as counts per million reads (CPM).

Participant 4 reportedly engaged in unprotected vaginal, oral, and anal sex throughout the study but also reported the use of condoms twice during the entire study. Participants 3 and 4 observed transitions between a BV state and intermediate state; the relative abundance of spp. increased as the study progressed. spp. consistently dominated the vaginal microbiota despite fluctuations in the relative abundance in the cases of consistent BV.

Participant 5 responded to treatment and had a healthy vaginal microbiota at month two but transitioned to BV at month eight, followed by intermediate BV at month 12. This participant consistently used condoms and had received Depo-Provera shots as forms of birth control throughout the study. She did not report oral, anal, or vaginal sex without a condom. Additionally, she had no new sex partners and had only one sex partner throughout the study period. This participant experienced the only instance of a healthy vaginal microbiota, which was dominated by . Markedly decreased levels of spp. were present in the healthy vaginal microbiota compared to intermediate and BV states. Once this participant transitioned from a healthy state to BV, the relative abundance of spp. increased and subsequently decreased after a transition from BV to the intermediate stage.

Fig. 2 shows the evolution of the most abundant bacterial species. We found the most abundant species associated with BV were spp., , Prevotella amnii, and Prevotella bivia. Gardnerella spp. increased from baseline over time. The relative abundance of and sp. decreased from baseline to the end of the study period. The relative abundance of the remaining bacterial species fluctuated throughout the 12 months. Contrary to expectations, there was an increase in the abundance of spp. over time with MET treatment.

Fig. 2.

Evolution of the most abundant taxa throughout the study. The boxplots show the 16 taxa whose abundance was the highest over all subjects. The different colours represent different time points. None of the changes were statistically significant due to the size of the sample.

Variance in vaginal microbiota

Fig. S1 (available in the online version of this article) illustrates the results of PCA. Based on these results, the attributes (taxa) with maximum variance (Table 1) were determined and investigated further. Consistent with results presented earlier in this section, the vaginal microbiota of Participant 5 (only participant with a healthy vaginal microbiota at any time point) was sufficiently distinct from the rest. Based on the PCA analysis results, spp. and species are strongly associated with BV, while and were associated with a healthy vaginal microbiota (Table 1). This was also supported by the fact that all species had a negative contribution in the last column of Table 1 (Principal component 1), while the BV-associated bacteria contributed to BV status.

Table 1.

Species determined by PCA to have maximum variance among five AA women

Taxa	Principal component 1
Gardnerella spp.	0.796033
Lactobacillus iners	−0.529284
Lactobacillus crispatus	−0.231106
Lactobacillus sp 7 1 47FAA	−0.172205
Megasphaera sp UPII 199 6	0.032557
Prevotella amnii	0.023059

spp.

sp 7 1 47FAA

sp UPII 199 6

Effect of metronidazole treatment

After the sPLS-DA analysis, the subjects’ data were plotted using their projections along with the first and second principal components (Fig. 3). Unlike PCA, sPLS-DA finds the direction that best separates the labelled classes. The samples were classified as untreated at time point 0 (baseline) and treated after receiving the treatment (months 2, 8 and 12). The treated group clusters together tightly, suggesting that treatment somehow ‘homogenizes’ the surviving microbiome. Subject 5 was a clear outlier. sPLS-DA determined that the genera, and Anaerococcus, were the most differentially abundant taxa between untreated and treated vaginal samples. However, this finding was likely due to overfitting the model as the Balanced Error Rate was 43%, close to the chance for a two-class classifier. Also, there was no statistically significant difference between treated and untreated samples for any of the taxa because of the small sample size. This argues that sPLS-DA findings should be interpreted with caution.

Fig. 3.

sPLS-DA projection of the samples into the first and second principal components (PC). Treated participants are coloured in blue while untreated ones are coloured in orange. Each sample is represented by X.MYY, where X means the participant number, and YY represents the month of data collection. The treated group clusters together in a tighter way than the heterogeneous untreated group (i.e. participants at baseline).

Resistome profile

The prevalence of resistance genes was analysed against the CARD database [34]. The metagenome of all the participants contained antibiotic resistance genes (Fig. 1d). The 15 most abundant genes analysis shows that tet genes that would confer resistance to tetracycline are most prevalent, with tetM being the most abundant, followed by tetQ. The next most abundant gene, lscA, is associated with clindamycin, lincomycin, dalfopristin, and tiamulin resistance. Analysis for the presence of nim genes related to the MET resistance showed the presence of nimJ in a few of the samples at extremely low levels.

Correlational studies

Fig. 4 shows the pairwise correlations between the taxa with the correlations displayed as a heatmap. Each row represents a different taxon. Rows and columns were ordered to depict the hierarchical clustering (obtained by iteratively grouping the most similar groups) shown as a dendrogram on the side of the heatmap. For readability, included only the strongest Spearman correlations (excluding self-correlations). A threshold of 0.8 (-0.8, resp.) was used for filtering positive (negative, resp.) correlation values. There were no correlations lower than −0.8. While these are not displayed in the figure, we note that the five strongest negative correlations were for the following taxa pairs: and spp. (−0.74), and (−0.69); N. unclassified and sp. Type 1 (-0.70); C. dubliniensis and sp. Type 1 (-0.64); sp. 7 47FAA and spp. (−0.62).

Fig. 4.

Top Spearman correlations for taxa pairs. Each interaction between a row and a column in this heatmap represents the Spearman correlation between the two taxa involved. Rows and columns are the different taxa ordered according to hierarchical clustering. For readability, only the top correlations were selected. Yellow colours indicate correlations closer to 1, while red colours imply lower correlations.

spp. and Staphylococcus spp. were highly correlated with many other species throughout the sample (Fig. 4). spp. was highly positively correlated with spp. and spp. Furthermore, was highly positively correlated with , , Actinomyces neuii, and unclassified. Of note is the correlation between and spp. as well as and spp.

Table 2 shows the top Spearman correlations between taxa abundance and the NS. There was just one episode of a low NS (healthy vaginal microbiota), four episodes with intermediate scores, and 14 episodes with high NSs (BV). Even though the data were small and skewed, the strongest correlations highlight the NSs underlying logic. All Lactobacilli spp. were associated with a low NS. In other words, as the abundance of Lactobacilli spp. increased, the NSs tend to decrease (negative association). The abundance of BV-associated bacteria such as , , Megasphaera, and were correlated with high NSs. Surprisingly, and spp. were absent from this top list, with Spearman correlation values in the range −0.35 through 0.34.

Table 2.

Spearman correlation between taxa and the Nugent score

Taxa name	Spearman correlation
Lactobacillus crispatus	−0.6347917188
Lactobacillus gasseri	−0.5947203609
Lactobacillus jensenii	−0.5345614474
Candida lusitaniae	−0.5345614474
Mycoplasma hominis	0.5236533323
Prevotella buccalis	0.5041102652 s
Lactobacillus amylovorus	−0.5000000000
Megasphaera sp. type 1	0.4949254787
Porphyromonas uenonis	0.4084817754
Lactobacillus vaginalis	−0.3983923509

Discussion

Due to the study’s pilot nature, the limited sample size, findings should be interpreted with caution and should not be generalized. Since this study depended on the participant self-reporting, one cannot rule out the participants' failure to adhere to the antibiotic regimen. Another confounding factor might be when the samples were collected; even though specific guidelines were provided for the sample collection, we cannot guarantee it was the same for each participant. With these caveats, three important observations are noted in asymptomatic BV patients treated with MET – (i) specific changes in the vaginal microbiota with treatment, (ii) recurrence and persistence of BV despite treatment, (iii) presence of resistance genes.

According to previous studies, effective treatment of BV resulted in a decrease of the relative abundance of facultative and strict anaerobes in the vaginal microbiota accompanied by an increase in the relative abundance of spp. [2, 39]. However, the current study showed a decrease over time in the relative abundance of spp. Similar findings were observed in a study conducted by Hummelen et al. [40] with a Tanzanian population. The MET treatment decreased species diversity in their study, but did not shift it to a BV-associated vaginal microbial profile dominated by a spp. [40]. The characteristic absence of among African populations compared to non-African could account for this observation [41].

Throughout the study, there was an observed increase in spp. for participants that did not maintain healthy vaginal flora. This increase may be due to the ineffectiveness of MET as a treatment option or increased spp. resistance. Xioa and colleagues reported findings that indicated spp. might be resistant to MET since the treatment exerts no influence on the bacteria at the molecular level [2]. Deng and colleagues also presented findings highlighting MET resistance in spp. [42]. Our resistome analyses demonstrated the presence of nimJ associated with MET resistance in few participants, albeit at extremely low level [43]. Though the nim alleles are widespread in both Gram-positive and -negative genera of aerobic and anaerobic bacteria and archaea, their presence in spp. has not been reported level [43]. Thus, it is not possible to infer if the presence of nimJ actually contributed to any resistance. In addition, majority of the samples harboured tetracycline resistance genes. The most abundant of the tet gene tetM has been associated with spp. [44]; and the gene lsaC associated with clindamycin resistance. Clindamycin, a lincosomide, has been used as an alternate to MET for BV treatment [10]. The participants of this study are likely to respond poorly to clindamycin.

At baseline, the only participant (Participant 5) with a relatively healthy vaginal microbiota had the lowest level of spp. which may have influenced treatment efficacy. This is also supported by the presence lowest level of tet genes. Furthermore, at baseline, all participants were classified as CST IV (Fig. 1), except for Participant 5, who was classified as CST III. Brooks and colleagues have determined that the ability to transition to another CST is determined by the current CST [45]. The difference in baseline CST could have facilitated the observed transition of Participant 5 to a healthy vaginal microbiota. Also, CST III is dominated by [45]. There was an inverse relationship between the abundance of spp. and L. iners, which was only observed for Participant 5 with a healthy vaginal microbiota (Fig. 1), another factor that could have potentially facilitated the transition to a healthy vaginal microbiota for this participant. It has also been determined that is the dominant species posttreatment [15]. Participant 5 was also the only participant to report the use of hormonal birth control, which has been shown to reduce the risk of BV [46]. The results of PCA and sPLS-DA further highlighted the outlier status of Participant 5 when compared to the other participants.

spp., L. iners, and were determined by PCA to contribute to maximum variance among the samples. The relative abundance of taxa was not homogeneous for all participants. This heterogeneity could have resulted from the use of different birth control methods and differences in sexual practices. The results of the Spearman correlation between taxa abundance and NSs were consistent with current literature [47]. The species, L. crispatus, L. gasseri, and L. jensenii, were shown to be inversely associated with the NS, which suggests that they promote a healthy vaginal microbiota. Conversely, the relative abundance of Mobilincus hominis increased with the increase in NSs. Though spp. promotes a healthy vaginal microbiota; not all Lactobacilli appeared to be protective. has been shown to protect against BV whereas has been found to be associated with BV [14]. Our studies showed a correlation between and spp. as well as and spp. This observation is interesting because when BV was first identified, it was believed that the causative agent was [21]. These three species appear highly correlated with many other species in the vaginal microbiota.

This study suggests MET may not be an effective course of treatment for women with a low relative abundance of spp. and maybe more effective when the vaginal microbiota is dominated by instead of . Additionally, in cases of MET treatment, there is a high probability of long-term recurrence. The findings presented in this study warrant further research. For future research, the vaginal microbiota should be examined at shorter time intervals since the changes that occur are dynamic. Daily assessments may provide more useful information than using 4 month intervals. In addition to whole-genome sequencing, metabolomics should be conducted to determine whether BV may result from a specific pathogen, its metabolites, or a combination of the pathogen and its metabolites. Finally, research should also consider the contributing factors of a shift from a BV-associated vaginal microbial profile one that is dominated by spp. after treatment with MET. Moreover, the resistome profile would better inform the treatment options.

Data availability

All sequence data generated in this study were deposited in NCBI under BioProject number: PRJNA669294 with Accession numbers SRR12830910 through SRR12830929.

Click here for additional data file.