A multi-platform metabolomics approach identifies highly specific biomarkers of bacterial diversity in the vagina of pregnant and non-pregnant women

Bacterial vaginosis (BV) increases transmission of HIV, enhances the risk of preterm labour, and is associated with malodour. Clinical diagnosis often relies on microscopy, which may not reflect the microbiota composition accurately. We use an untargeted metabolomics approach, whereby we normalize the weight of samples prior to analysis, to obtained precise measurements of metabolites in vaginal fluid. We identify biomarkers for BV with high sensitivity and specificity (AUC = 0.99) in a cohort of 131 pregnant and non-pregnant Rwandan women, and demonstrate that the vaginal metabolome is strongly associated with bacterial diversity. Metabolites associated with high diversity and clinical BV include 2-hydroxyisovalerate and γ-hydroxybutyrate (GHB), but not succinate, which is produced by both Lactobacillus crispatus and BV-associated anaerobes in vitro. Biomarkers associated with high diversity and clinical BV are independent of pregnancy status, and were validated in a blinded replication cohort from Tanzania (n = 45), where we predicted clinical BV with 91% accuracy. Correlations between the metabolome and microbiota identified Gardnerella vaginalis as a putative producer of GHB, and we demonstrate production by this species in vitro. This work illustrates how changes in community structure alter the chemical composition of the vagina, and identifies highly specific biomarkers for a common condition.

In most instances, diagnosis is dependent upon microscopy of vaginal samples to identify BV-like bacteria by morphology alone (Nugent Scoring 11 ), or in combination with clinical signs (Amsel Criteria 12 ). The precision and accuracy of these methods are poor due to the diverse morphology of vaginal bacteria, the observation that many women with BV are asymptomatic, and subjectivity in microscopic examination [13][14][15] . Misdiagnosis creates stress for the patient, delays appropriate intervention and places a financial burden on the health care system. A rapid test based on stable, specific biomarkers for BV would improve diagnostic accuracy and speed, and reduce costs through improved patient management.
The metabolome, defined as the complete set of small molecules in a given environment, has been studied in a variety of systems to identify biomarkers of disease 16,17 , and advance our understanding of how the microbiota contributes to host metabolism 18 . Using an untargeted multiplatform metabolomics approach, combined with 16S rRNA gene sequencing, we demonstrate that the vaginal metabolome is driven by bacterial diversity, and identify biomarkers of clinical BV that can be reproduced in a blinded validation cohort. We further demonstrate that Gardnerella vaginalis, which has long been thought to be an important contributor to BV, is the likely source of one of the most specific compounds, GHB. This work provides a foundation for improved detection of disease and demonstrates how metabolomics can be utilized to identify validated sources of metabolites in microbial communities.

Results
The vaginal metabolome is most correlated with bacterial diversity. We completed a comprehensive untargeted metabolomic analysis of vaginal fluid in two cross-sectional cohorts of Rwandan women: pregnant (P, n = 67) and non-pregnant (NP, n = 64) (Supplementary Table S1). To normalize the amount of sample collected, vaginal swabs were weighed prior to and after collection and normalized to equivalent concentrations. This enabled us to collect precise measurements of metabolites in vaginal fluid. Metabolite profiling was carried out using both gas chromatography-mass spectrometry (GC-MS) and liquid chromatography-mass spectrometry (LC-MS), and microbiota composition by 16S rRNA gene sequencing.
The metabolome determined by GC-MS contained 128 metabolites (Supplementary Table S2). We conducted a series of partial least squares (PLS) regression analyses to determine the single variable that could best explain the variation in the metabolome. In both cohorts, the diversity of the microbiota, as measured using Shannon's Diversity 19 , was the factor that explained the largest percent variation in the metabolome (Supplementary Table S3), demonstrating that the vaginal metabolome is most correlated with bacterial diversity (Fig. 1A). Metabolites robustly associated with this diversity (95% CI) (Fig. 1B) were determined by jackknifing, and within this group, metabolites associated with extreme diversity tended to have less variation in the jackknife replicates, and were common to both pregnant and non-pregnant women. This identified a core set of metabolites associated with diversity.
The two cohorts overlapped by principal component analysis (PCA) ( Supplementary Fig. S1), and no metabolites were significantly different between pregnant and non-pregnant women (unpaired t-test, Benjamini-Hochberg p > 0.01). Thus, the cohorts were combined for all further analysis.
Metabolites and taxa associated with diversity. A single PLS regression was performed on all samples with Shannon's diversity as a continuous latent variable ( Supplementary Fig. S2). Samples were then ordered by their position on the 1 st component of this PLS. The diversity indices, microbiota and metabolites associated with diversity of PLS ordered samples are shown in Fig. 2. The vaginal microbiota of Rwandan women were similar to women from other parts of the world, with the most abundant species being L. iners followed by L. crispatus [1][2][3]20 (Fig. 2B, Supplementary Table S4). Women with high bacterial diversity were dominated by a mixture of anaerobes, including Gardnerella, Prevotella, Sneathia, Atopobium, Dialister and Megasphaera species. Figure 2D displays metabolites robustly associated with bacterial diversity in both cohorts based on the PLS loadings in Fig. 1B. Metabolites associated with high diversity include amines, which contribute to malodor [16][17][18] , and a number of organic acid derivatives such as 2-hydroxyisovalerate (2HV), γ -hydroxybutyrate (GHB), 2-hydroxyglutarate and 2-hydroxyisocaproate. Low diversity was characterized by elevated amino acids, including the amine precursors lysine, ornithine and tyrosine. Many of these metabolites were also detected by LC-MS, and trimethylamine (high diversity) and lactate (low diversity) were detected exclusively by this method (Supplementary Table S5). The identities of metabolites of interest were confirmed with authentic standards when available (Fig. 2, asterisks).

Succinate is not associated with diversity or clinical BV, and is produced by L. crispatus.
Succinate and lactate abundance are shown in panel E of Fig. 2. Succinate levels, and the succinate:lactate ratio have historically been associated with BV [21][22][23] , and succinate has been postulated to play an immunomodulatory role 23 . Here we show that succinate is not associated with bacterial diversity, nor is it significantly elevated in clinical BV as defined by Nugent scoring. This trend was independent of the detection method used. In addition, succinate was elevated in women dominated by L. crispatus compared to L. iners-dominated women (unpaired t-test, Benjamini-Hochberg p < 0.01) ( Supplementary Fig.  S3), indicating L. crispatus may produce succinate in vivo, a phenomenon that has been demonstrated in vitro 24 . We extracted metabolites from vaginal isolates grown on agar plates and confirmed that succinate Scientific RepoRts | 5:14174 | DOi: 10.1038/srep14174 is produced by L. crispatus in vitro, but not by L. iners ( Supplementary Fig. S4). Succinate was also produced by Prevotella bivia, and Mobiluncus curtisii, but not by G. vaginalis.
Metabolites associated with diversity are sensitive and specific for clinical BV. We defined clinical BV by the Nugent method, which is the current gold standard for BV diagnosis 11 . This microscopy-based technique defines BV as a score of 7-10 when low numbers of lactobacilli morphotypes are observed, and high numbers of short rods are present, which are presumed to represent BV associated bacteria. Nugent Normal (N) is defined as a score of 1-3, indicating almost exclusively Lactobacillus morphotypes. Intermediate samples are given a score of 4-6 and do not fit into either group. Although Nugent scores correlated well with bacterial diversity in our study, it was apparent from the microbiota and metabolome profiles that two samples (41 and 145) had clearly been misclassified by Nugent ( Fig. 2A, red dots). The Nugent status of these samples was therefore corrected prior to further analyses.
In total we identified 49 metabolites that were significantly different between clinical BV and N (unpaired t-test, Benjamini-Hochberg p < 0.01, Supplementary Table S2). We determined the odds ratio (OR) for BV based on conditional logistic regressions of all individual metabolites detected by GC-MS (Supplementary Table S2) to determine if the metabolites we associated with high bacterial diversity could accurately identify clinical BV as defined by Nugent scoring. Metabolites significantly elevated in Nugent BV (unpaired t-test, Benjamini-Hochberg p < 0.01) with OR > 1 are shown in Fig. 3A. Succinate was included as a comparator, although it did not reach significance. Both GHB and 2HV were significantly higher in women with BV, and had OR > 2.0, demonstrating they are indicators not only of high bacterial diversity, but also clinical BV. Receiver operating characteristics (ROC) curves built from LC-MS data determined that high 2HV, high GHB, low lactate and low tyrosine were the most sensitive and specific biomarkers for BV, with the largest area under the curve (AUC) achieved using the ratio of 2HV:tyrosine (AUC = 0.993) (Fig. 3B-D). ROC curves of GC-MS data identified similar trends, with the largest AUC achieved by the ratio of GHB:tyrosine (AUC = 0.968) (Supplementary Table S6).
We determined the optimal cut points for the GHB:tyrosine (0.621) and 2HV:tyrosine (0.882) ratios by selecting values which maximized the sensitivity and specificity for BV. These cut points were generated excluding Nugent intermediate samples, however when cut points were applied to intermediates, they grouped equally with N or BV, and samples with smaller proportions of lactobacilli tended to group  Table S7). The GHB:tyrosine ratio cut point was slightly less specific (88%), with an AUC of 0.948. We confirmed that succinate was not significantly different between Nugent N and BV in the validation set, nor was the succinate:lactate ratio.

Identification of G. vaginalis as a producer of GHB.
Correlations between metabolites and the OTU abundances were performed using a method that took into account both the compositional nature of 16S rRNA gene survey data and the technical variation [26][27][28] . Metabolites and taxa which contained any correlation below a Benjamini-Hochberg corrected p < 0.01 are displayed as a heatmap in Fig. 6. Tyramine, putrescine, and cadaverine were most correlated with Dialister (Spearman's R = 0.54, 0.51, 0.61, p < 0.001) (Supplementary Table S8), indicating this genus may contribute to malodor. We found that GHB was most correlated with G. vaginalis (Spearman's R = 0.56, p < 0.001), while 2HV was most correlated with Dialister, Prevotella, and Gardnerella (Spearman's R = 0.55, 0.48, 0.47, p < 0.001).   We chose to investigate the correlation between GHB and G. vaginalis, since this was an unexpected metabolite that was predictive for both Shannon's diversity and Nugent BV. Examination of available genomes showed that many strains of G. vaginalis possess a putative GHB dehydrogenase (annotated as 4-hydroxybutyrate dehydrogenase). We extracted metabolites from bacterial colonies grown on agar plates and reproducibly detected GHB in G. vaginalis extracts well above control levels (unpaired t-test, p < 0.05), but did not detect GHB from other species commonly associated with BV (Fig. 7, Supplementary Table S9). These data suggest that G. vaginalis is the primary source of GHB detected in vivo.

Discussion
We have demonstrated that alterations in the vaginal metabolome are driven by bacterial diversity in both pregnant and non-pregnant Rwandan women, and identified 2HV and GHB as highly specific biomarkers of clinical BV, the latter of which we attribute to production by G. vaginalis. We obtained  extremely accurate results by controlling for the mass of vaginal fluid collected, however we recognize this may not be logistically possible in a clinical setting. To circumvent this need we expressed biomarkers as ratios to the amino acid tyrosine, which we identified as the most differential amino acid in health (Supplementary Table S2). Using optimal cut points of these ratios we predicted 91% of Nugent BVs in a blinded replication cohort, demonstrating the reproducibility of our findings. These cut points also accurately classified Nugent intermediate samples into groups with similar microbiota profiles dominated by either lactobacilli (N) or anaerobes (BV).
Although we demonstrate production of GHB by G. vaginalis, it is important to note that no single organism has been identified as the cause of BV, and G. vaginalis is present in many women with a lactobacilli-dominated microbiota. However, GHB is metabolized from succinate in other bacteria 29,30 , suggesting a similar pathway could exist in G. vaginalis. Succinate-producing genera may therefore be required, making G. vaginalis essential, but not sufficient for GHB production in the vagina. This remains to be tested.
The predominant pathway for succinate production in bacteria is from pyruvate via anaerobic respiration. The genes for this pathway are expressed in vivo and differentiate BV from N 31 . Srinivasan et al. 32 recently proposed an alternate pathway whereby succinate is produced from putrescine via gamma-aminobutyric acid (GABA). Although this pathway is plausible, we find it unlikely given many of the enzymes are either not expressed in vivo or are absent from the genomes of common vaginal organisms 31 .
Despite previous findings that succinate is elevated in BV 21-23,32 it was not differential in our study. This unexpected outcome is likely a result of normalizing sample weights prior to analysis, which we used to ensure the most consistent measurements of metabolites. Succinate was one of the most abundant metabolites detected in vaginal fluid in our study ( Supplementary Fig. S2), and was present in nearly all samples regardless of BV status. The universal presence of succinate make it more susceptible to dilution effects compared to GHB and 2HV, which were less abundant and below our detection limit in many non-BV samples. Other groups have reported large ranges in succinate abundance in women with BV 21,22 , or used pooled samples 22 , which could account for additional disparities in results. Differences in succinate abundance may have been more pronounced in previous studies if there were a lack of L. crispatus-dominated women, which our data demonstrate is a succinate producer (Fig. S3,  S4). L. crispatus contains all the enzymes necessary to produce succinate from malate with the exception of malate dehydrogenase (MDH). However, the pathway is annotated as complete at http://biocyc.org/ LCRI491076-HMP/missing-rxns.html, with the closely related enzyme lactate dehydrogenase (LDH). As there is increased expression of succinate-producing pathways during BV 30 , it is probable that large amounts are produced initially, but then rapidly converted to other compounds, such as GHB, by the microbiota and/or host.
In addition to GHB, 2HV was identified as a highly specific biomarker for BV. 2HV is produced from breakdown of branched chain amino acids in humans 33 and some bacteria [34][35][36] . When the trend for amino acid depletion in BV is considered, these findings suggest increased amino acid catabolism in this condition. Some of these amino acids are converted to the amines cadaverine, tyramine, and putrescine, which are also associated with BV. These odor-causing compounds were most correlated with Dialister. Yeoman et al. 37 also linked amines to Dialister species, and the decarboxylating genes required for amine production are expressed by this genus in vivo 30 . These data strongly suggest that Dialister is one of the genera responsible for malodor in the vagina. Given the small proportion of this genus in women with BV (0.2-8% in our study), this emphasizes the need for functional characterizations of the microbiome using metabolomic and transcriptomic approaches.
The taxa that constitute the vaginal microbiota are highly conserved across different populations 1-3,20,38 , although prevalence of certain taxa differs between ethnicities 1,38 . These observations lead us to believe that GHB, 2HV and their tyrosine ratios will be globally applicable for the diagnosis of BV. Our ability to replicate findings in a distinct population strongly supports this theory. Srinivasan et al. 32 concurrently identified elevated GHB in the vaginal fluid of American women with BV 39 , however they were not able to replicate this result in a second cohort. This could be due in part to the use of cervicovaginal lavages for sample collection or the use of different detection methods between cohorts. 2HV (annotated as alpha-hydroxyisovalerate) was also identified as differential in their study, but was not tested in the replication cohort. These observations further validate our findings and demonstrate these biomarkers are robust to the effects of dilution.
The exact role, if any, of GHB and 2HV in the etiology of BV is unknown. Systemically GHB has both inhibitory and excitatory effects through activation of the GABA(B) and perhaps GABA(A) receptors in the brain, resulting in stimulatory and sedative effects if taken at high doses [40][41][42] . The effects of GHB at other sites remain elusive. Future work should attempt to elucidate biological function of GHB and other novel metabolites to determine what effect (if any) they have on lactobacilli and the vaginal environment.
Although we did not identify any metabolites that differed significantly between pregnant and non-pregnant women, it should be noted that patient selection was biased to include an even number of women with and without Nugent BV. Our study was not designed to test if the metabolome differed during pregnancy, but rather if the metabolic signatures of BV were similar between pregnant and non-pregnant women. Other groups have noted decreased bacterial diversity during pregnancy and across gestational age [43][44][45] . These observations suggest differences in the metabolome of pregnant women Scientific RepoRts | 5:14174 | DOi: 10.1038/srep14174 would be observable in a larger randomly sampled population, and may include elevated levels of metabolites associated with low diversity such as amino acids.
In summary, we have demonstrated using an untargeted, multiplatform approach that differences in the vaginal metabolome are driven by bacterial diversity. Other metabolomic studies have focused on symptom-associated metabolites 32,37 , changes after treatment 46 , or longitudinal changes in a few subjects 47 , and included exclusively non-pregnant women. We identified several highly specific biomarkers for clinical BV that are independent of pregnancy status, and replicated this result in a blinded cohort. By combining high-throughput sequencing with advanced mass spectrometry techniques we have shown how in vivo metabolite information can be used to identify validated sources of metabolic end products in bacterial communities. These techniques can be applied to many systems where organisms may be fastidious or difficult to culture, and provide a much-needed link between microbial composition and function.

Methods
Clinical samples. Premenopausal women between the ages of 18 and 55 were recruited at the University of Kigali Teaching Hospital (CHUK) and the Nyamata District Hospital in Rwanda. The Health Sciences Research Ethics Board at The University of Western Ontario, Canada, and the CHUK Ethics Committee, Rwanda granted ethical approval for all experiments involved in the study. The methods were carried out in accordance with the approved guidelines and all women provided written informed consent. Participants were excluded if they had reached menopause, had a current infection of gonorrhoea, Chlamydia, genital warts, active genital herpes lesions, active syphilis, urinary tract infections, received drug therapy that may affect the vaginal microbiome, had unprotected sexual intercourse within the past 48 hours, used a vaginal douche, genital deodorant or genital wipe in past 48 hours, had taken any probiotic supplement in past 48 hours, or were menstruating at time of clinical visit. As materials for sample collection were limited, we set out to obtain an equal number of women with and without Nugent BV to ensure the study would be powered to test for BV biomarkers. To accomplish this, only women with suspected BV were recruited after the quota of Nugent N women was met. After reviewing details of the study, participants gave their signed consent before the start of the study. For metabolome analysis, sterile Dacron polyester-tipped swabs (BD) were pre-cut with sterilized scissors and weighed in 1.5 ml microcentrifuge tubes prior to sample collection. Using sterile forceps to clasp the pre-cut swabs, a nurse obtained vaginal samples for metabolomic analysis by rolling the swab against the mid-vaginal wall. A second full-length swab was obtained for Nugent Scoring and 16S rRNA gene sequencing using the same method. Nugent Scoring was performed at CHUK by Amy McMillan. Vaginal pH was measured using pH strips. Samples were frozen within 2 hours of collection and stored at − 20 °C or below until analysis.

Microbiome profiling.
Vaginal swabs for microbiome analysis were extracted using the QIAamp DNA stool mini kit (Qiagen) with the following modifications: swabs were vortexed in 1 mL buffer ASL before removal of the swab and addition of 200 mg of 0.1 mm zirconia/silica beads (Biospec Products). Samples were mixed vigorously for 2 × 30 seconds at full speed with cooling at room temperature between (Mini-BeadBeater; Biospec Products). After heating to 95 °C for 5 minutes, 1.2 ml of supernatant was aliquoted into a 2ml tube and one-half an inhibitEx tablet (Qiagen) was added to each sample. All other steps were performed as per the manufacturers instructions. Sample amplification for sequencing was carried out using the forward primer (ACACTCTTTCCCTACACGACGCTCTTCCGATCTnnnn(8)CWACGCGARGAACCTTACC) and the reverse primer (CGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCTn(12)ACRACACGAGCTG ACGAC) where nnnn indicates four randomly incorporated nucleotides, and (8) was a sample nucleotide specific barcode. The 5′ end is the adapter sequence for the Illumina MiSeq sequencer and the sequences following the barcode are complementary to the V6 rRNA gene region. Amplification was carried out in 42 μ L with each primer present at 0.8 pMol/mL, 20 μ L GoTaq hot start colorless master mix (Promega) and 2 μ L extracted DNA. The PCR protocol was as follows: initial activation step at 95 °C for 2 minutes and 25 cycles of 1 minute 95 °C, 1 minute 55 °C and 1 minute 72 °C.
All subsequent work was carried out at the London Regional Genomics Centre (LRGC, lrgc.ca, London, Ontario, Canada). Briefly, PCR products were quantified with a Qubit 2.0 Flourometer and the high sensitivity dsDNA specific fluorescent probes (Life Technologies). Samples were mixed at equimolar concentrations and purified with the QIAquick PCR Purification kit (QIAGEN). Samples were paired-end sequenced on an Illumina Mi-Seq with the 600 cycle version 3 reagents with 2 × 220 cycles. Data was extracted from only the first read, since it spanned the entirety of the V6 region including the reverse primer and barcode.
Resulting reads were extracted and de-multiplexed using modifications of in-house Perl and UNIX-shell scripts with operational taxonomic units (OTUs) clustered at 97% identity, similar to our reported protocol 48 . Automated taxonomic assignments were carried out by examining best hits from comparison the Ribosomal Database Project (rdp.cme.msu.edu) and manually curated by comparison to the Greengenes database (greengenes.lbl.gov) and an in house database of vaginal sequences (Macklaim unpublished). Taxa with matches at least 95% similarity to query sequences were annotated as such. OTUs were summed to the genus level except for lactobacilli, and rare OTUs found at less than 0.5% abundance in any sample removed. Supplementary Table S1 displays the nucleotide barcodes and their corresponding samples. Reads were deposited to the Short Read Archive (BioProject ID: PRJNA289672). To control for background contaminating sequences, a no-template control was also sequenced. Barplots were constructed with R {r-project.org } using proportional values.
To avoid inappropriate statistical inferences made from compositional data, centred log-ratios (clr), a method previously described by Aitchison 49 and adapted to microbiome data was used with paired t-tests for comparisons of genus and species level data 27,28 . The Benjamini Hochberg (False Discovery rate) method was used to control for multiple testing with a significance threshold of 0.1. All statistical analysis, unless otherwise indicated, was carried out using R (r-project.org).

Sample Preparation GC-MS.
Vaginal swabs were pre-cut into 1.5 mL tubes and weighed prior to and after sample collection to determine the mass of vaginal fluid collected. After thawing, swabs were eluted in methanol-water (1:1) in 1.5 mL microcentrifuge tubes to a final concentration of 50 mg vaginal fluid/ mL, which corresponded to a volume ranging from 200-2696 μ L, depending on the mass of vaginal fluid collected. A blank swab eluted in 800 μ L methanol-water was included as a negative control. All samples were vortexed for 10 s to extract metabolites, centrifuged for 5 min at 10 621 g, vortexed again for 10 s after which time the brushes were removed from tubes. Samples were centrifuged a final time for 10 min at 10 621 g to pellet cells and 200 μ L of the supernatant was transferred to a GC-MS vial. The remaining supernatant was stored at − 80 °C for LC-MS analysis. Next, 2 μ L of 1 mg/mL ribitol was added to each vial as an internal standard. Samples were then dried to completeness using a SpeedVac. After drying, 100 μ L of 2% methoxyamine-HCl in pyridine (MOX) was added to each vial for derivatization and incubated at 50 °C for 90 min. 100 μ L N-Methyl-N-(trimethylsilyl) trifluoroacetamide (MSTFA) was then added and incubated at 50 °C for 30 min. Samples were then transferred to micro inserts before analysis by GC-MS (Agilent 7890A GC, 5975 inert MSD with triple axis detector). 1 μ L of sample was injected using pulsed splitless mode into a 30 m DB5-MS column with 10 m duraguard, diameter 0.35mm, thickness 0.25 μ m (JNW Scientific). Helium was used as the carrier gas at a constant flow rate of 1 ml/min. Oven temperature was held at 70 °C for 5 min then increased at a rate of 5 °C/min to 300 °C and held for 10 min. Solvent delay was set to 13 min to avoid solvent and a large lactate peak, and total run time was 61 min. Masses between 25 m/z and 600 m/z were selected by the detector. All samples were run in random order and a standard mix containing metabolites expected in samples was run multiple times throughout to ensure machine consistency.
Data Processing GC-MS. Chromatogram files were deconvoluted and converted to ELU format using the AMDIS Mass Spectrometry software 50 , with the resolution set to high and sensitivity to medium. Chromatograms were then aligned and integrated using Spectconnect 51 (http://spectconnect.mit.edu), with the support threshold set to low. All metabolites found in the blank swab, or believed to have originated from derivatization reagents were removed from analysis at this time. After removal of swab metabolites, the IS matrix from Spectconnect was transformed using the additive log ratio transformation (alr) 49 and ribitol as a normalizing agent (log2(x)/log2(ribitol)). Zeros were replaced with two thirds the minimum detected value on a per metabolite basis prior to transformation. All further metabolite analysis was performed using these alr transformed values.
Metabolites were initially identified by comparison to the NIST 11 standard reference database (http:// www.nist.gov/srd/nist1a.cfm). Identities of metabolites of interest were then confirmed by authentic standards if available.
Global metabolomic analysis. In order to visualize trends in the metabolome as detected by GC-MS, principal component analysis (PCA) was performed using pareto scaling. To determine the percentage of variation in the metabolome that could be explained by a single variable we performed a series of partial least squares (PLS) regressions where each variable was used as a continuous latent variable. We tested every taxa, pH, Nugent score, pregnancy status, Shannon's diversity index and sample ID and compared the percent variation explained by the first component of each PLS. The variable with the highest value was determined to be most closely associated with the metabolome (Shannon's Diversity). Analysis was conducted in R using the PLS package and unit variance scaling. Jackknifing with 20% sample removal and 10 000 repetitions was then applied to determine 95% confidence intervals for each metabolite. Metabolites with confidence intervals that did not cross zero in both cohorts (pregnant and non-pregnant) were considered significantly associated with diversity. Heatmaps of significant metabolites were constructed using the heatmap.2 function in R with average linkage hierarchical clustering and manhattan distances. Unless specified otherwise, all tests for differential metabolites between groups were performed using unpaired t-test with a Benjamini-Hochberg (False Discovery Rate) significance threshold of p < 0.01 to account for multiple testing and multiple group comparisons.
16S rRNA microbial gene profiles generate compositional data that interferes with many standard statistical analyses, including deter determining correlations [26][27][28] . We used the aldex.corr function from the ALDEx2 package to calculate the Spearman's rank correlation between each OTU abundance in 128 inferred technical replicates and that were transformed by center log-ratio transform 27 cor.test function in R. This approach is conceptually similar to that adopted by SPARCC 26 , but calculates the correlation between the OTU abundances and continuous metadata variables. Heatmaps of correlation p values were constructed using the heatmap.2 function in R with complete linkage hierarchical clustering and Euclidean distances.
Odds ratios of metabolites to identify Nugent BV from Normal were calculated from conditional logistic regressions performed on all metabolites using the glm function in R with 10 000 iterations and a binomial distribution. Metabolites with 95% CI > 1 and p < 0.01 (unpaired t-test, Benjamini-Hochberg corrected) were determined to be significantly elevated in Nugent BV. "Nugent BV" was defined by the clinical definition of a score of 7-10, with a score of 0-3 being "Nugent Normal". ROC curves and forest plots were built in R using the pROC and Gmisc packages respectively.

Sample Preparation LC-MS.
To confirm GC-MS findings, samples which had at least 100 μ L remaining after GC-MS were also analyzed by LC-MS. 100 μ L of supernatant was transferred to vials with microinserts and directly injected into an Agilent 1290 Infinity HPLC coupled to a Q-Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific) with a HESI source. For HPLC, 2 μ L of each sample was injected into a ZORBAX Eclipse plus C18 2.1 × 50mm × 1.6 micron column. Mobile phase (A) consisted of 0.1% formic acid in water and mobile phase (B) consisted of 0.1% formic acid in acetonitrile. The initial composition of 100% (A) was held constant for 30 s and decreased to 0% over 3.0 min. Mobile phase A was then held at 0% for 1.5 minutes and returned to 100% over 30s for a total run time of 5 min.
Full MS scanning between the ranges of m/z 50-750 was performed on all samples in both positive in negative mode at 140 000 resolution. The HESI source was operated under the following conditions: nitrogen flow of 25 and 15 arbitrary units for the sheath and auxiliary gas respectively, probe temperature and capillary temperature of 425 °C and 260 °C respectively and spray voltage of 4.8 kV and 3.9 kV in positive and negative mode respectively. The AGC target and maximum injection time were 3e6 and 500 ms respectively. For molecular characterization, every tenth sample was also analyzed with a data dependent MS2 method where a 35 000 resolution full MS scan identified the top 10 signals above a 8.3e4 threshold which were subsequently selected at a 1.2 m/z isolation window for MS2. Collision energy for MS2 was 24, resolution 17 500, AGC target 1e5 and maximum injection time was 60ms. Blanks of pure methanol were run between every sample to limit carryover, and a single sample was run multiple times with every batch to account for any machine inconsistency. A blank swab extract was also run as a negative control.
For increased sensitivity, a separate LC-MS method was used for relative quantification of GHB in human samples. This was accomplished by selected ion monitoring in the mass range of 103.1-107.1 m/z in positive mode, and integrating the LC peak area of the [M + H + ] ion (± 5 ppm).

Data Processing LC-MS.
After data acquisition Thermo RAW files were converted to MZML format using ProteoWizard 53 and imported into MZmine 2.11 54 (http://mzmine.sourceforge.net) for chromatogram alignment and deconvolution. Masses were detected using the Exact Mass setting and a threshold of 1E5. For Chromatogram Builder, minimum time was 0.05 min, minimum height 3E3, and m/z threshold set to 0.025 m/z or 8 ppm. Chromatogram Deconvolution was achieved using the Noise Amplitude setting with the noise set to 5E4 and signal to 1E5 for negative mode. Due to an overall greater signal and noise in positive mode, the noise was adjusted to 6E5 and signal to 6.5E5 for positive mode. Join aligner was used to combine deconvoluted chromatograms into a single file with the m/z threshold set to 0.05 m/z or 10 ppm, weight for m/z and RT set to 20 and 10 respectively, and a RT tolerance of 0.4 min. After chromatograms were aligned, a single .CSV file was exported and all further analysis was carried out in R.
To confirm metabolites identified as significant by GC-MS in the LC-MS data set, the masses of metabolites of interest were searched in the LC-MS data set, and identities confirmed by MS 2 using METLIN 55 and the Human Metabolome Database 56 online resources. Standards of metabolites of interest were also run to confirm identities when available. An unpaired t-test with Benjamini-Hochberg correction was used to determine metabolites significantly different between Nugent BV and Normal in the LC-MS data set. Metabolites with corrected p < 0.05 were considered statistically significant. Metabolites detected exclusively by LC-MS that have previously been associated with BV or health (lactate, trimethylamine) were also included in this analysis. Data was log base 10 transformed prior to data analysis and zeros replaced by two thirds the minimum detected value on a per metabolite basis. To determine optimal cut points of biomarkers for diagnostic purposes, cut points were computed from LC-MS data using the OptimalCutpoints package in R 57 and the Youden Index method 58 . Validation in blinded replication cohort. Women between the ages of 18 and 40 were recruited from an antenatal clinic at the Nyerere Dispensary in Mwanza, Tanzania as part of a larger study on the effect of micronutrient supplemented probiotic yogurt on pregnancy. The Medical research Coordinating Committee of the National Institute for Medical Research (NIMR), as well as the Health Sciences Research Ethics Board at The University of Western Ontario granted approval for all experiments involved in the study. The methods were carried out in accordance with the approved guidelines and all women provided written informed consent. The study was registered with clinicaltrials.gov (NCT02021799). Samples were collected using the methods mentioned above, and Nugent scores performed by research technicians at NIMR in Mwanza, Tanzania. A subset of samples was selected based on these Nugent scores by a third party, who ensured there was not repeated sampling of any women. Amy McMillan, who performed metabolite analysis, was blinded to the Nugent scores for the duration of sample processing and data analysis. Biomarkers were quantified in samples by LC-MS using the protocols mentioned above. The study was unblinded after the submission of BV status based on the ratio cut points established in the Rwandan data set.
Identification of putative GHB dehydrogenases in G. vaginalis strains. The protein sequence of a bona fide 4-hydroxybutyrate (GHB) dehydrogenase isolated from Clostridium kluyveri 29 (GI:347073) was blasted against all strains of G. vaginalis in the NCBI protein database. Blast results identified multiple isolates containing a putative protein with 44-46% identity to the GHB dehydrogenase from C. kluyveri. The strain used for in vitro experiments (G.vaginalis ATCC 14018) was not present in the NCBI protein database, however a nucleotide sequence in 14018 with 100% nucleotide identity to a putative 4-hydroxybutyrate dehydrogenases in strain ATCC 14019 (GI:311114893) was identified, indicating potential for GHB production by strain 14018.
In vitro extraction of GHB from vaginal isolates. Due to their fastidious nature, we found it difficult to obtain consistent growth of all vaginal strains in liquid media. To circumvent this, a lawn of bacteria was plated and metabolites were extracted from agar punches. All strains were grown on Columbia Blood Agar (CBA) plates using 5% sheep's blood for 96h under strict anaerobic conditions, with the exception of L. crispatus, which was grown on de Man Rogosa Sharp (MRS) agar for 48 h. To extract metabolites, 16 agar punches 5 mm in diameter were taken from each plate and suspended in 3 mL 1:1 Me:H 2 0. Samples were then sonicated in a water bath sonicater for 1h, transferred to 1.5 ml tubes after vortexing and spun in a desktop microcentrifuge for 10 min at 10 621 g to pellet cells. 200 μ l of supernatant was then aliquoted for GC-MS described above. The area of each peak was integrated using ChemStation (Agilent) by selecting m/z 233 in the range of 14-16 min. Initial peak width was set to 0.042 and initial threshold at 10. An authentic standard of GHB was run with samples to confirm identification. Succinate production by vaginal isolates was measured from the same GC-MS run, and quantified using Spectconnect as described above. Un-inoculated media was used as a control and experiments were repeated three times with technical duplicates.