Introduction
The vaginal microbiome appears to play an important role in many aspects of reproduction, insofar as its composition is associated with tubal factor infertility, with early miscarriage, the gestational age at delivery, and with neonatal health outcomes of term infants.1 Vaginal lactobacilli are typically associated with positive pregnancy outcomes, whereas bacterial vaginosis (BV), a poorly-defined polymicrobial syndrome, is associated with miscarriage, increased risk of preterm birth (PTB), and other adverse outcomes of pregnancy.2–5 BV is generally characterized by a shift in the vaginal microbiota from Lactobacillus dominance and a low pH environment to taxonomic diversity and a higher pH environment. Despite the association between BV and PTB, treatment of BV for the prevention of PTB has not been consistently efficacious,6–8 suggesting either that additional factors play a role, or that the microbiota sets into motion a host response that fails to respond to antimicrobial intervention.
Once thought to be an immune-suppressed state, necessary for tolerating the semi-allogeneic conceptus, successful pregnancy is now understood to be immunologically complex. Critical transitions between distinct immunological phases shift the Th1/Th2 T helper cell balance either towards Th1, when pro-inflammatory signals dominate, or towards Th2, termed the “Th2 phenomenon” when anti-inflammatory signals are dominant.9 A modest pro-inflammatory environment contributes to successful implantation followed by a gradual increase in Th2 cytokines including interleukin (IL)-4 and IL-10 that suppress pro-inflammatory cytokines in the 2nd trimester.10 As term approaches, anti-inflammatory IL-10 levels decrease resulting in an increase in pro-inflammatory cytokine production, prostaglandin (PG) synthesis, and the induction of labor.11,12
IL-1β is thought to be one of the main drivers in activating the inflammatory processes of labor.13–15 Together with tumor necrosis factor-α (TNFα), IL-1β modulates matrix metalloproteinase (eg, MMP-8 and -9) activity,16 which mediates the collagen degradation required for cervical ripening and weakening of the fetal membranes.17,18 Pro-inflammatory cytokines IL-1β, IL-6, and IL-8 also contribute to PG synthesis by activating cyclooxygenase-2 expression via nuclear factor kappaB (NFκB) signaling.11 PGs prostaglandin E2 (PGE2) and PGF2α induce spontaneous uterine contractility and contribute to cervical ripening and decidua/membrane activation.19 This interplay between cytokine and eicosanoid mediators is regulated reciprocally at both the synthetic and catabolic level such that while pro- and anti-inflammatory cytokines modulate expression of cyclooxygenase-2and production of 15-hydroxyprostaglandin dehydrogenase to stimulate synthesis and reduce catabolism of PGs,11 the pro-inflammatory prostanoid, PGE2, and anti-inflammatory PG, PGD2 can stimulate or repress the release of inflammatory cytokines in gestational tissues.11,19,20 Sphingolipids have also been implicated in the inflammatory processes of normal parturition. During pregnancy sphingolipid metabolism is highly activated in the decidua21 and the ceramide, sphingosine and sphingosine-1-phosphate, can induce synthesis of inflammatory prostanoids PGE2 and PGI2, respectively.22–24
Although a normal component of parturition, inflammation is also implicated in preterm delivery.25,26 While the number of leukocytes in the vagina during most BV cases is not increased,28 pro-inflammatory cytokines have been reported to be elevated in pregnant and non-pregnant women with BV.27,29,30 Likewise, eicosanoids and sphingolipids may be altered in vaginal dysbiotic conditions. In a study of the metabolic signature of BV, the non-prostanoid eicosanoid 12-hydroxyeicosatetraenoic (HETE) was identified as a key metabolite positively associated with BV-associated bacteria and sphingosine was found to be differentially associated with BV and BV-associated bacteria.31 The potential for vaginal bacteria to cause intrauterine infection and/or modulate the highly interrelated signaling pathways critical to successful pregnancy, parturition, and labor is thus suggestive of a hypothesis wherein the host inflammatory responses initiated by vaginal dysbiotic conditions such as BV and/or microbial invasion of the amniotic cavity can contribute to preterm labor and birth. However, it is essential to go beyond defining the vaginal microbiome and to probe host-microbiome signaling to gain a better understanding of the microbial players, temporal aspects of pathology, and the nature of the induced host response associated with pathological outcomes.
Herein, we describe a cross-sectional early pregnancy study under the umbrella of the National Institutes of Health Human Microbiome Project (HMP) to explore the potential of multi-omic data to describe the vaginal microenvironment during pregnancy. Rationalizing the clinical impact of early detection of conditions that impact risk of adverse outcomes, vaginal samples collected from women in their first trimester of pregnancy were subjected to 16S rRNA gene survey, targeted lipidome and multiplexed human cytokine analyses. These data were analyzed in tandem to present an integrated snapshot of the microbial profiles and the associated host response in early pregnancy. Taxonomic surveys revealed, as prior studies have shown, that microbial profiles cluster into multiple distinct community types. The majority of women exhibited community types dominated by Lactobacillus species; however, Lactobacillus-dominated community types were more prevalent among pregnant women relative to non-pregnant women. Using discriminant analysis methods, integration of multi-omic datasets revealed significant associations between microbial taxa and host lipid and cytokine signaling. The results support multifactor analyses as a means to garner valuable insight into the vaginal environment during pregnancy.
Results
Participant demographics
Sixty-six women who were enrolled in the Multi Omic Microbiome Study-Pregnancy Initiative (MOMS-PI study),32 a project supported by the second phase of National Institutes of Health Human Microbiome Project, the Integrative HMP,33 during the first trimester of pregnancy, and 81 non-pregnant women who were enrolled in the Vaginal Human Microbiome Project (VaHMP) study,34 a project supported by the first phase of the NIH HMP, and who were matched for demographic characteristics, were selected for analysis of the vaginal microbiome and comparison of community profiles during early pregnancy and during non-pregnancy. A minimum of 1000 reads were available from the 16S rRNA gene survey analysis of each of the included subjects. Demographic and relevant health history information is summarized in Table 1.
;)
Table 1: Demographic summary and selected health history of women selected for this study.
The vaginal microbiome of early pregnancy is dominated by Lactobacillus
Comprehensive 16S rRNA gene-based taxonomic survey of the vaginal microbiota of the MOMS-PI participants yielded a mean count of 55,474 reads/sample with minimum and maximum read counts of 1792 and 178,837, respectively. Vaginal microbiome profiles of non-pregnant participants and participants in their first trimester are depicted in Figure 1 A and B, respectively. As previously described,35 the most abundant taxonomic group with at least 30% prevalence was used to define vaginal community types. Profiles for which no taxonomic group was present at >30% were assigned “No Type.” Alpha diversity within the community types was calculated using Simpson reciprocal index (1/SI) and there were significant differences across these vaginal community types, Kruskal-Wallis P-value = 1.7e-04 and 6.9e-05 (Figure 2A and B for the non-pregnant and pregnant cohorts, respectively). Post-hoc tests revealed that Simpson reciprocal index for samples with vaginal community types characterized by Gardnerella vaginalis and Atopobium vaginae in the pregnant cohort was significantly higher than Lactobacillus crispatus dominated communities (P-values = 0.00092 and 0.008, respectively). Table 2 lists the prevalence of the major vaginal community types identified in this pregnant cohort in order of abundance, and compared with a non-pregnant cohort with similar demographic characteristics. Sixty-six percent (n = 43) of samples from pregnant women were dominated by Lactobacillus species while only 4% (n = 2) were designated “No type.” Within the non-pregnant cohort, 64% exhibited Lactobacillus-dominated vaginal community types while none were designated “No type.” Additionally, only 6% of pregnant women and 11% of healthy non-pregnant women had community profiles dominated by G. vaginalis.
Figure 1: Vaginal microbiome profiles of early pregnancy. Stacked barplot showing relative abundance of bacterial taxa representing ≥0.1% relative abundance in at least 5% of samples in A: VaHMP cohort (n = 81 participants) and B: MOMS-PI cohort (n = 66 participants) (x-axis). The profiles are grouped by the most abundant species and are ordered by decreasing proportion of the dominant bacterium. A color code showing the most abundant taxa is shown as an annotation bar on the top. VaHMP: vaginal human microbiome project.
Figure 2: Alpha diversity differs significantly across community types. Boxplot of Simpson reciprocal index (1/SI) across vaginal microbiome community types. A: VaHMP cohort. B: MOMS-PI cohort. Post-hoc pairwise comparisons are shown in Supplementary Table S2. MOMS-PI: multi-omic microbiome study-pregnancy initiative; VaHMP: vaginal human microbiome project.
Table 2: Abundance of vaginal microbiome profiles.
Abundance of specific taxa correlate with specific cytokines and lipids through pairwise analysis
To assess the interactions between microbiome, lipids, and cytokines, we used targeted lipidomics to quantitate 10 eicosanoids and 32 sphingolipid species and a multiplexed enzyme-linked immunosorbent assay to quantify 24 cytokines in the vaginal samples. These targeted lipidomic and immuno-proteomic datasets were evaluated in parallel with the 16S rRNA gene-based taxonomic surveys. Bacterial taxa present at relative abundance of at least 0.1% and present in at least 5% of samples were included in the analysis. Figure 3 shows the largest linear discriminant analysis effect sizes (LEfSe) for bacterial taxa associated with being above-the-mean or below-the-mean level for the cytokines tested. Only linear discriminant analysis (LDA) scores above 3.5 are reported. L. crispatus was strongly associated with lower levels of a number of cytokines, including ILs 12 and 17 (IL-12, IL-17), macrophage inflammatory protein (MIP1α), and vascular endothelial growth factor (VEGF). In contrast, BV-associated taxa, such as Sneathia, G. vaginalis, and BVAB1 were associated with very different cytokine profiles. For example, granulocyte colony stimulating factor (G-CSF) levels were lower when these taxa were abundant, but IL-9, IL-10, and TNF-α levels were higher. Figure 4 shows LEfSe results for bacterial taxa associated with being above-the-mean or below-the-mean lipid levels. Prevotella cluster 2 was associated with below-the-mean levels of a number of different lipids including the pro-inflammatory eicosanoid, PGE2, and sphingomyelins C24:1, C22; C18:1, and C14, and others. BV-associated species, G. vaginalis, A. vaginae, and Megasphaera were associated with above-the-mean levels of the eicosanoid precursor eicosapentaenoic acid (EPA).
Figure 3: Association between bacterial taxa and cytokine levels. Cytokine levels were dichotomized to above-the-mean (green, positive) and below-the-mean (red, negative). LEfSe analysis detected significant associations between the abundance of a number of bacterial taxa and cytokine levels. LEfSe: largest linear discriminant analysis effect sizes.
Figure 4: Association between bacterial taxa and lipid levels. Lipid levels were dichotomized to above-the-mean (green, positive) and below-the-mean (red, negative). LEfSe analysis detected significant differences in the abundance of a number of bacterial taxa and certain lipids. LEfSe: largest linear discriminant analysis effect sizes.
Bacterial taxa correlate with specific lipids and cytokines through multiple comparison analysis
Spearman correlation coefficients for lipid: cytokine pairs were calculated and are listed in Supplementary Table S1 (Supplemental Digital Content, https://links.lww.com/IMD/A1). To investigate the significant correlations between these lipids and cytokines in the context of community type, we assembled scatterplots in which scatter was based on lipid and cytokine levels and the line of best fit was based on community type. Results from this analysis are summarized in Table 3. Several lipid: cytokine pairs for which the community type predicted a correlation between the lipid and cytokine were detected. Significant correlations between C18-1 SM and G-CSF and MIP1α corresponded with predominance of A. vaginae. Similarly, significant correlation between EPA and MIP1α and C16 DH SM and VEGF corresponded with predominance of Lactobacillus gasseri and BVAB1, respectively. The data suggests there to be a deterministic association between community type and lipid:cytokine signatures.
Table 3: Community types in which lipid: cytokine pairs correlated significantly.
To further explore natural patterns of variation across participants’ samples and better understand the contribution of combining multiple omic platforms, we performed a principal component analysis (PCA)-based multiple co-inertia analysis (MCIA). MCIA identifies relationships between multiple high-dimensional datasets.36,37 The first and second axes of MCIA (Figure 5A) explained ∼49% of variation across multi-omic datasets of samples from all 58 participants. Most participant samples were grouped near the origin of the sample projection plot; although several individuals were markedly divergent. There was no apparent clustering of participant samples by vaginal microbiome community type (Figure 5A). However, variable projections on the first and second MCIA axes revealed interesting groupings. While variables were mostly clustered by omic technology, several groupings consisted of variables from two or more omics (Figure 5B). Pro-inflammatory cytokines VEGF and IL-12 appeared to be grouped near taxon Coriobacteriaceae OTU27 as well as the eicosanoids Prostaglandin J2 and PGE2. Likewise, pro-inflammatory cytokines IL-6 and MIP-1β formed a second group with taxa Clostridiales OTU22, L. gasseri, Lachnospiraceae OTU33, Peptoniphilus lacrimalis, Campylobacter ureolyticus, Actinomycetales OTU158, Finegoldia magna, and Dialister propionifaciens. These clustering patterns are suggestive of relationships between variables from different omic datasets.
Figure 5: Multiple co-inertia analysis identifies relationships among multi-omic variables. Sample (A) and variable (B) projection from multiple co-inertia analysis (MCIA) using PCA to integrate the omics datasets. A: Samples are colored by vaginal microbiome community type. B: Variables are colored by omic dataset. Clusters of variables from three omic datasets are circled and labeled. PCA: principal component analysis.
Discussion
Cytokine and lipid mediators have critical roles in the maintenance of pregnancy and the inflammatory processes characteristic of labor in normal pregnancy. As such, dysregulation of these key mediators could have negative impacts on the outcome of pregnancy and on the timing of labor onset. While the impact of the vaginal immune status on pregnancy is not yet entirely clear, dendritic cells from the vagina can migrate to the nodes draining the uterus and can thereby modulate the immune status of the upper reproductive tract.38 Evidence suggests that signaling cascades initiated during BV and/or intrauterine infections can lead to the increased levels of inflammatory cytokines such as IL-1β, IL-6, IL-8, and TNFα that have been shown to be risk factors for preterm labor and birth.39,40 However, the impact of specific vaginal bacterial taxa on the host is difficult to study. Tissue culture models are compromised by their over-simplicity in that they usually employ a single cell type and they fail to factor in all of the compounds, such as hormones and other signaling factors, to which host tissues are exposed. Animal models are suboptimal because animals, even primates, have a very different vaginal microbial profile than humans, and most of the bacterial taxa associated with the human vagina do not colonize efficiently. Additionally, human vaginal microbiomes are never composed of a single bacterial taxon and it is unethical to introduce a potential pathogen into a human. Thus, to garner useful information about the impact of bacteria on the human host, increased sample size is necessary to power discovery of biologically meaningful associations between bacterial taxa and host responses, and omics techniques are required to cover a larger swath of microbial and host signals. The multi-omic strategy of the Integrative Human Microbiome Project effort seeks to use such large-scale approaches to explore relationships between the microbiome and the host response. In this study, using data from the MOMS PI, we characterized the vaginal microbiomes of a cohort of 66 women and the microbiomes, lipidomes, and cytokine profiles of a subset of 58 of these women within the first trimester of pregnancy to obtain a snapshot of the microbial and host components at play in early pregnancy that may provide insight into the mechanisms and processes associated with adverse pregnancy outcomes. Moreover, insight into these mechanisms early in pregnancy may permit early intervention to prevent adverse outcomes.
We observed similar abundance of Lactobacillus-dominated community types among our pregnant and non-pregnant cohorts. This may be because the women were sampled early in pregnancy, because prior studies suggest that Lactobacillus-dominated community types are more common among pregnant women41–43 and evidence suggests that Lactobacillus dominance may promote a positive pregnancy outcome.44–47 The H2O2-producing species of lactobacilli are regarded as “healthy” vaginal bacteria due to their ability to maintain a low vaginal pH environment and inhibit growth of other bacterial taxa.48 Indeed, other studies have shown,49,50 and data from this study supports, that diversity levels vary significantly across vaginal community types, with those dominated by lactobacilli exhibiting low diversity and those dominated by G. vaginalis and A. vaginae exhibiting higher diversity.
In addition to regulating the vaginal microbiome, lactobacilli may also play a role in regulating the host inflammatory response. LEfSe analysis revealed the large effect size of the association between L. crispatus and below-the-mean levels of IL-12, IL-17, MIP1α, and VEGF. Studies have failed to show an association between IL-12 and PTB, but elevated tissue levels of IL-12 have been associated with recurrent early pregnancy loss.51,52 IL-17 is a potent pro-inflammatory cytokine that is associated with Th17 T cells and chronic inflammatory conditions.53 MIP1α is a chemokine produced by activated macrophage that induces neutrophil infiltration. VEGF plays important roles in both angiogenesis and in T cell trafficking and both of these roles impact implantation and pregnancy, although the role of VEGF in the vagina is not clear.54,55 Thus, healthy lactobacilli like L. crispatus may abrogate cellular immune infiltration within the vagina.
BV-associated bacteria appear to modulate the genital tract inflammatory response in a very different way than healthy lactobacilli.56–58 Vaginal microbial communities of high diversity are associated with pro-inflammatory immune responses including increased levels of pro-inflammatory cytokines.59,60 Analysis of our pregnant cohort revealed associations between specific BV-associated taxa and above-the-mean levels of pro-inflammatory cytokines including IL-9, MIP1α, and TNFα. They were also associated with decreased levels of other cytokines, including G-CSF and IL-1ß. IL-9 and MIP1α are primarily produced by Th2 cells including Natural Killer cells and activated macrophages, respectively, and increased levels of these cytokines may represent infiltration of these immune cells. However, the roles of IL-9 and MIP1α in BV and during pregnancy have not been well characterized although increased TNFα levels have been detected in association with spontaneous PTB.61
Though most bacteria, with the exception of species of Sphingobacterium, Sphingomonas, Bacteroides, and Bdellovibrio stolpii, do not contain sphingolipids and cannot synthesize eicosanoids due to lack of cyclooxygenase activity, many taxa can alter host lipid metabolism through expression of microbial determinants such as bacterial sphingomyelinase (SMase) and phospholipase A2.62–64 Furthermore, in cervico-vaginal fluid of women with BV, acid SMase activity is increased and directly correlated with levels of IL-1β.63 As such, we expected vaginal microbiome communities to be associated with distinct lipid profiles. Prevotella cluster 2 stood out among taxa as being associated with below-the-mean levels of a number of lipids. Most reads within this cluster display highest identity to Prevotella timonensis. P. timonensis has not been well characterized, although analysis of the shotgun sequencing project available through National Center for Biotechnology Information (NCBI), revealed a number of putative esterases. Thus, it is plausible that this species produces a SMase or other esterases capable of modifying host lipids. BV-associated taxa, G. vaginalis, A. vaginae, and Megasphaera were associated with above-the-mean levels of the EPA and 6-Keto-PGF1α. EPA is an eicosanoid precursor and is thus related to lipid-mediated inflammatory processes. The lipid 6-Keto-PGF1α is a metabolite of and a marker of prostacyclin, a PG eicosanoid that mediates vasodilation and inhibits platelet aggregation. Therefore, BV-associated bacteria may influence lipid-mediated inflammation. While the mechanism and biologic significance of these associations between taxa and lipids were not explored in this study, it is the first to compare these two omics in the vagina and demonstrates that specific associations exist.
Linear regression of scatter plots grouped by community type detected communities dominated by major taxa that correlated with levels of certain cytokines and lipids while MCIA was able to detect correlations between minor taxa, lipids, and cytokines. This demonstrates the utility of multiple pairwise analyses and supports the use of different analytic methods to examine multi-omic data sets.
The sub-cohort for this study was assembled early in the MOMS-PI project as a pilot study, and consequently, pregnancy outcomes were not known at the time of the analyses. The relatively small size of the cohort also precluded associations between pregnancy outcome and taxa, cytokines, and lipids. Another limitation of the study is the single time-point at which samples were collected from each subject. Development of a predictive tool for adverse outcomes early in pregnancy might be best informed by longitudinal sampling. In fact, our recent work suggests that microbial signatures may be more evident in early stages of pregnancy.32,65 The MOMS-PI pregnant cohort and VaHMP non-pregnant cohort were sequenced using different platforms66 and were subject to batch effects, thus we cannot rule out the possibility that technical sources of variation may have affected this comparative analysis. However, the DNA extraction method and the template-specific portion of the 16S rRNA gene-specific primers were conserved between studies and a variety of protocols have observed robust detection of Lactobacillus Community State Types across a wide range of vaginal microbiome surveys.67 Finally, a caveat of all such studies is that there may be biases in the technologies; for example, preferential lysis of bacteria, primer bias, and so on, that impact the results. However, we validated our protocols using mock samples containing panels of bacteria commonly found in the female reproductive tract and attempted to quantify inherent bias.68
The results of this study support the hypothesis that the vaginal microbiome is an important modulator of the host immunological response, and that its impact can be assessed early in pregnancy. Future studies integrating multi-omic datasets from larger cohorts with substantial numbers of negative pregnancy outcomes may reveal associations between microbial potential and host cytokine and lipidome profiles that comprise signatures associated with adverse outcomes, and ultimately promote women's health during pregnancy.
Materials and methods
Participants
The MOMS-PI, a part of the National Institutes of Health Integrative Human Microbiome Project Consortium (http://hmp2.org/), was designed to investigate the dynamic changes in vaginal and related microbiomes and associated host responses during normal pregnancy, complicated pregnancy, and early infancy using a multi-omic strategy.33 The MOMS-PI cohort includes women at least 15 years of age who were not incarcerated or surrogate mothers at their first presentation in the Virginia Commonwealth University (VCU) Health Center Women's Clinics. Participants willing to provide samples throughout their pregnancy and delivery were enrolled following written, informed consent, and asked to complete comprehensive questionnaires and health histories, including information about their race and ethnicity, level of education, employment, sexual history, past pregnancies, health, and dietary habits. A subset of the MOMS-PI cohort, composed of 66 women between the ages of 18 and 37 from whom vaginal samples had been collected before 14 weeks gestation, and for whom 16S rRNA gene taxanomic profiles were complete and yielded ≥1000 reads, were selected. Of this subset, lipidomics datasets and cytokine profiles were complete for 58 women.
The VaHMP included participants recruited from outpatient clinics at the VCU Medical Center and the Virginia Department of Health following written, informed consent from 2009 to 2013. Participants filled out a detailed questionnaire that included questions about ethnicity, education, employment, health habits, dietary habits, and sexual history. Clinicians also filled out a diagnosis form at the time of each visit that included information about the purpose of each visit, and any diagnoses. Inclusion criteria for VaHMP included women age 18–44 years who were able to provide informed consent, who were visiting the clinic for annual exam, and who were willing or already scheduled to undergo a vaginal examination using a speculum. A cohort of 81 non-pregnant women who were matched to the MOMS-PI subset in regards to age, race, and BMI, were selected to compare the prevalence of community state types among pregnant (MOMS-PI) versus non-pregnant (VaHMP) women.
The Institutional Review Board (IRB) for Human Subjects Research at VCU reviewed and approved these studies, IRB HM15527 MOMS-PI and IRB HM12169 The Vaginal Microbiome: Disease, Genetics, and the Environment. All methods were carried out in accordance with the guidelines and regulations outlined in the approved protocols and informed consent was obtained from all participants.
Sampling and sample processing
During a new obstetric visit, four consecutive sets of vaginal swab (double-tipped CultureSwab EZ [Becton Dickinson, Franklin Lakes, NJ]) samples were collected by a physician from the vagina without a speculum, or by self-sampling. For collection, physicians and participants were instructed to spread the skin around the vagina, insert the swab ∼2 inches, and rotate it against the vaginal wall for 5 seconds. Swabs from the first set were designated for DNA extraction and swabs from the third set were used for lipid and cytokine analysis. Remaining swabs were collected for other analyses.
Swabs for DNA extraction, cytokine, and lipid analysis were immediately suspended in 750 μL of Powersoil Bead Solution (MoBio Laboratories, Inc., Hilden, Germany), 500 μL 100 mM Tris-HCl pH 7.5, or 500 μL 0.01% butyl hydroxy toluene in phosphate buffered saline (PBS), respectively, and processed within 1 hour of collection. Samples designated for DNA extraction and cytokine analysis were stored at −80°C until further processing. Lipid samples were analyzed within 2 hours. DNA was extracted in a 96-well plate using the Powersoil® DNA Isolation Kit (MoBio) according to the manufacturer's instructions. Total protein concentration from cytokine swabs and post-amplification DNA concentration from DNA swabs were compared by collection type (self vs clinician) and are reported in S1 Figure.
16S rRNA gene-based taxonomic survey
The V1–V3 hypervariable regions of the bacterial 16S rRNA gene was polymerase chain reaction (PCR) amplified using barcoded primers.69–71 Briefly, for each reaction, 2 μL of extracted DNA was combined with 12.5 μL 2X Phusion Hot Start II High-Fidelity PCR Master Mix (Invitrogen/Thermo Fisher Scientific, Carlsbad, CA), 3% dimethylsulfoxide and 100 nM each of forward and reverse primers. The primers included an Illumina linker adaptor (used for binding and sequencing), a unique index sequence followed immediately by a variable sequence spacer (0–6 bases) and 16S rRNA gene primers (see S1 Appendix). The forward primer was a mixture (4:1) of primers Fwd-P1 and Fwd-P2. The PCR was carried out in a 25 μL reaction in a Thermal Cycler (Applied BioSystem GeneAmp PCR system 9700) with the following parameters: initial denaturation at 98°C for 30 seconds, followed by 30 cycles of 98°C for 15 seconds, 58°C for 15 seconds, and 72°C for 15 seconds with a final extension at 72°C for 1 minutes. The PCR was performed in 96-well format with two PCR controls, a negative water control and a positive MOMS-PI Mock Community control. All amplicons were quantified using Picogreen (Invitrogen/Thermo Scientific) and a spectrofluorimeter (Biotek; a part of Agilent, Winooski, VT). Equal amounts of amplicon were combined followed by removal of unincorporated primers, salts, and enzymes using Agencourt AMpure XP beads. The DNA concentration of this concentrated pool was verified by qPCR using the KAPA Library Quantification Kit for Illumina platform using Thermo Fisher Scientific ViiA™ 7 Real-Time PCR System. The library pool was diluted and denatured according to the Illumina MiSeq library preparation guide. The sequencing run was conducted on the Illumina MiSeq sequencer using the 2 × 300 PE reagent kit 3. Sequence reads were de-multiplexed using an in-house Python script. The raw overlapping paired-end sequence data were merged and quality-filtered using MeFiT72 with meep (maximum expected error rate) cut-off of 1.0. Non-overlapping high-quality paired end reads were retained for downstream analysis by linking them artificially with 15 N's. High-quality species-level taxonomic assignments were performed using STIRRUPS,73 employing the USEARCH algorithm74 combined with a curated vaginal 16S rRNA gene database.
Targeted lipid analysis
For eicosanoid and sphingolipid quantification, an equal volume of ethanol containing 10 ng of eicosanoid internal standards or methanol containing 50 pmol of sphingolipid internal standards were added to clarify vaginal swab contents dispersed in PBS. Eicosanoid internal standards consisted of 30 deuterated analytes including, (d4) 6-keto-PGF1α, (d4) PGF2α, (d4) PGE2, (d4) PGD2, (d4) LTB4, (d4) TXB2, (d4) LTC4, (d5) LTD4, (d5) LTE4, (d8) 5-HETE, (d8) 15-HETE, (d8) 14,15 epoxyeicosatrienoic acid, (d8) arachidonic acid, and (d5) EPA. Sphingolipid internal standards consisted of d17 sphinogsine, sphinganine, sphingosine-1-phosphate, sphinganine-1-phosphate and d18:1/12:0 ceramide-1-phosphate, sphingomyelin, ceramide, and monohexosylceramide (Avanti). After centrifugation at 12,000g for 20 minutes, the resultant mixture was subjected to UPLC ESI-MS/MS analysis using a hybrid triple quadrupole linear ion trap mass analyzer (ABSCIEX 6500 QTRAP®) via multiple-reaction monitoring. Detailed separation, elution and ionization conditions have been previously described75 and are summarized in S1 Appendix. Spectral data were analyzed using MultiQuant software (AB SCIEX) and quantitation was carried out via comparison against known quantities of internal standards.
Cytokine analysis
The BIO-RAD Bio-Plex Pro™ Human Cytokine 27-Plex Assay was employed with a Luminex® 100/200™ System to quantify cytokine levels in the vaginal samples. Frozen vaginal swab samples suspended in 500 μL 100 mM Tris-HCl pH 7.5 were thawed on ice and centrifuged at 10,000g for 10 minutes at 4°C. The Bio-plex assay was conducted on samples and serial dilutions of standards in duplicate according to the manufacturer's instructions. Values were analyzed using a five-parameter logistic non-linear regression curve model. Cytokine values deemed out of range were assigned the upper or lower limit of detection for the specific cytokine. IL-2, IL-5, and IL-15 levels were below the limit of detection in >50% vaginal samples and were not included in subsequent analyses. Though biological replicates could not be performed due to the effects of multiple freeze-thaw cycles on sample viability, intra-assay coefficient of variation was found to be acceptable at <7% for all cytokines. Total protein concentration for each sample was determined using a Bradford Assay.76 Cytokine concentration (pg/mL) was divided by total protein concentration (mg/mL) to yield normalized cytokine concentration (pg cytokine/mg protein) for each sample. Samples for which total protein could not be determined were not included in the analysis.
Statistical analyses
Total above threshold reads (reads with average ≥97% identity) were used to calculate relative abundances. Cytokine values relative to total protein (pg/mg total protein) and lipid values (ng/μL) were normalized to Z-scores. For unsupervised analyses, MCIA using PCA for table ordination was performed to integrate scaled and centered 16S rRNA-based taxonomic, eicosanoid, sphingolipid, and cytokine datasets across all subjects. MCIA projects several datasets into lower-dimensional spaces and transforms diverse sets of variables or features onto the same scale.36,37 MCIA uses the same optimization criterion and constraints as the closely related method consensus principal component analysis.77 In a two-step process, MCIA first transforms individual datasets into comparable lower-dimensional spaces using a table ordination method (such as PCA) followed by a generalized co-inertia analysis that projects multiple datasets into the same maximally covariant dimensional space.36,37 As such, samples and variables sharing similar trends are closely projected in the same hyperspace, which can highlight the lack or presence of co-structure between datasets and permit identification of biomarkers and clusters of samples.36,78 All statistical analyses were performed with R.79 R package ade480 was used for MCIA.
LEfSe applies a Kruskal-Wallis rank-sum test for each bacterium, then uses linear discriminant analysis to estimate effect size.81 The effect size is the contribution of a variable to the ability to distinguish two different groups. To discover relationships between bacterial taxa and cytokines or lipids, cytokine and lipid levels were dichotomized into above-the-mean and below-the-mean. Bacterial taxa associated with being above-the-mean were given a positive LDA score, while the LDA score was converted to a negative number for taxa associated with being below-the-mean.
Data availability
All sequence, lipid and cytokine data generated from the MOMS-PI project will be hosted and made available to the public through the Integrative Human Microbiome Project Data Coordination Center (http://ihmpdcc.org). MOMS-PI data will be accessible through NCBI BioProject PRJNA326441. Vaginal microbiome, cytokine and lipid profiles used in this study can also be found on the Vaginal Microbiome Consortium website (http://vmc.vcu.edu/resources/momspi).
Acknowledgments
In addition to named authors, we acknowledge additional members of the Vaginal microbiome Consortium (VMC) who contributed to this work including Abdallah A. Abdelmaksoud, Niels C. Asmussen, Joseph F. Borzelleca, Jamie L. Brooks, Samuel Boundy, Charles E. Chalfant, Molly R. Dickinson, Jennifer I. Drake, Robert A. Duckworth, Joe El Khoury, Abigail L. Glascock, Michael G. Gravett, Karen D. Hendricks-Muñoz, Nicole R. Jimenez, Ana M. Lara, Vladimir Lee, Andrey V. Mateyev, Luiz V Ozaki, Victor V. Pokhilko, Ronald M. Ramus, Sarah K. Rozycki, Craig E. Rubens, Amber S. Sexton, Stephany C Vivadelli, Niran R. Wijesooriya, Jie Xu. Sequence analysis was performed in the Genomics Core of the Nucleic Acids Research Facilities at VCU. Lipidomic analysis was provided by the VCU Lipidomics/Metabolomics Core Facility. The staff of the Bioinformatics Computational Core Laboratories at VCU provided Bioinformatics analysis.
References
1. Haahr T, Zacho J, Bräuner M, Shathmigha K, Skov Jensen J, Humaidan P. Reproductive outcome of patients undergoing in vitro fertilization treatment and diagnosed with bacterial vaginosis or abnormal vaginal microbiota: a systematic PRISMA review and meta-analysis.
BJOG Int J Obstet Gynaecol 2019; 126 (2):200–207.
2. Donders GG, Van Calsteren K, Bellen G, et al. Predictive value for preterm birth of abnormal vaginal flora, bacterial vaginosis and aerobic vaginitis during the first trimester of pregnancy.
BJOG Int J Obstet Gynaecol 2009; 116 (10):1315–1324.
3. Dingens AS, Fairfortune TS, Reed S, Mitchell C. Bacterial vaginosis and adverse outcomes among full-term infants: a cohort study.
BMC Pregnancy Childbirth 2016; 16 (1):278.
4. Nelson DB, Hanlon AL, Wu G, Liu C, Fredricks DN. First trimester levels of BV-associated bacteria and risk of miscarriage among women early in pregnancy.
Matern Child Health J 2015; 19 (12):2682–2687.
5. Ugwumadu AHN. Bacterial vaginosis in pregnancy.
Curr Opin Obstet Gynecol 2002; 14 (2):115–118.
6. Carey JC, Klebanoff MA, Hauth JC, et al. Metronidazole to prevent preterm delivery in pregnant women with asymptomatic bacterial vaginosis. National Institute of Child Health and Human Development Network of Maternal-Fetal Medicine Units.
N Engl J Med 2000; 342 (8):534–540.
7. Haahr T, Ersbøll AS, Karlsen MA, et al. Treatment of bacterial vaginosis in pregnancy in order to reduce the risk of spontaneous preterm delivery – a clinical recommendation.
Acta Obstet Gynecol Scand 2016; 95 (8):850–860.
8. Lamont RF, Keelan JA, Larsson PG, Jørgensen JS. The treatment of bacterial vaginosis in pregnancy with clindamycin to reduce the risk of infection-related preterm birth: a response to the Danish Society of Obstetrics and Gynecology guideline group's clinical recommendations.
Acta Obstet Gynecol Scand 2017; 96 (2):139–143.
9. Wilczyński JR. Th1/Th2 cytokines balance--yin and yang of reproductive immunology.
Eur J Obstet Gynecol Reprod Biol 2005; 122 (2):136–143.
10. Mor G, Cardenas I, Abrahams V, Guller S. Inflammation and pregnancy: the role of the immune system at the implantation site.
Ann N Y Acad Sci 2011; 1221 (1):80–87.
11. Keelan JA, Blumenstein M, Helliwell RJA, Sato TA, Marvin KW, Mitchell MD. Cytokines, prostaglandins and parturition--a review.
Placenta 2003; 24: (Suppl A): S33–S46.
12. Bowen JM, Chamley L, Mitchell MD, Keelan JA. Cytokines of the placenta and extra-placental membranes: biosynthesis, secretion and roles in establishment of pregnancy in women.
Placenta 2002; 23 (4):239–256.
13. Markovic D, Bari MF, Lu B, Vatish M, Grammatopoulos DK. Corticotropin-releasing hormone interacts with interleukin-1β to regulate prostaglandin H synthase-2 expression in human myometrium during pregnancy and labor.
J Clin Endocrinol Metab 2013; 98 (7):2864–2875.
14. Markovic D, Vatish M, Gu M, et al. The onset of labor alters corticotropin-releasing hormone type 1 receptor variant expression in human myometrium: putative role of interleukin-1beta.
Endocrinology 2007; 148 (7):3205–3213.
15. Grammatopoulos DK, Hillhouse EW. Basal and interleukin-1beta-stimulated prostaglandin production from cultured human myometrial cells: differential regulation by corticotropin-releasing hormone.
J Clin Endocrinol Metab 1999; 84 (6):2204–2211.
16. Roh CR, Oh WJ, Yoon BK, Lee JH. Up-regulation of matrix metalloproteinase-9 in human myometrium during labour: a cytokine-mediated process in uterine smooth muscle cells.
Mol Hum Reprod 2000; 6 (1):96–102.
17. Choi S-J, Jung K-L, Oh S, Kim J-H, Roh C-R. Cervicovaginal matrix metalloproteinase-9 and cervical ripening in human term parturition.
Eur J Obstet Gynecol Reprod Biol 2009; 142 (1):43–47.
18. Geng J, Huang C, Jiang S. Roles and regulation of the matrix metalloproteinase system in parturition.
Mol Reprod Dev 2016; 83 (4):276–286.
19. Sykes L, MacIntyre DA, Teoh TG, Bennett PR. Anti-inflammatory prostaglandins for the prevention of preterm labour.
Reprod Camb Engl 2014; 148 (2):R29–R40.
20. Denison FC, Kelly RW, Calder AA, Riley SC. Secretory leukocyte protease inhibitor concentration increases in amniotic fluid with the onset of labour in women: characterization of sites of release within the uterus.
J Endocrinol 1999; 161 (2):299–306.
21. Kaneko-Tarui T, Zhang L, Austin KJ, et al. Maternal and embryonic control of uterine sphingolipid-metabolizing enzymes during murine embryo implantation.
Biol Reprod 2007; 77 (4):658–665.
22. Ballou LR, Chao CP, Holness MA, Barker SC, Raghow R. Interleukin-1-mediated PGE2 production and sphingomyelin metabolism. Evidence for the regulation of cyclooxygenase gene expression by sphingosine and ceramide.
J Biol Chem 1992; 267 (28):20044–20050.
23. Kim JI, Jo EJ, Lee H-Y, et al. Sphingosine 1-phosphate in amniotic fluid modulates cyclooxygenase-2 expression in human amnion-derived WISH cells.
J Biol Chem 2003; 278 (34):31731–31736.
24. Damirin A, Tomura H, Komachi M, et al. Sphingosine 1-phosphate receptors mediate the lipid-induced cAMP accumulation through cyclooxygenase-2/prostaglandin I2 pathway in human coronary artery smooth muscle cells.
Mol Pharmacol 2005; 67 (4):1177–1185.
25. Brocklehurst P. Infection and preterm delivery.
BMJ 1999; 318 (7183):548–549.
26. Goldenberg RL, Hauth JC, Andrews WW. Intrauterine infection and preterm delivery.
N Engl J Med 2000; 342 (20):1500–1507.
27. Thurman AR, Kimble T, Herold B, et al. Bacterial vaginosis and subclinical markers of genital tract inflammation and mucosal immunity.
AIDS Res Hum Retroviruses 2015; 31 (11):1139–1152.
28. Giraldo PC, de Carvalho JBJ, do Amaral RLG, da Silveira Gonçalves AK, Eleutério J, Guimarães F. Identification of immune cells by flow cytometry in vaginal lavages from women with vulvovaginitis and normal microflora.
Am J Reprod Immunol N Y N 19892012; 67 (3):198–205.
29. Mattsby-Baltzer I, Platz-Christensen JJ, Hosseini N Rosén. IL-1beta, IL-6, TNF alpha, fetal fibronectin, and endotoxin in the lower genital tract of pregnant women with bacterial vaginosis.
Acta Obstet Gynecol Scand 1998; 77 (7):701–706.
30. Cauci S, Guaschino S, De Aloysio D, et al. Interrelationships of interleukin-8 with interleukin-1beta and neutrophils in vaginal fluid of healthy and bacterial vaginosis positive women.
Mol Hum Reprod 2003; 9 (1):53–58.
31. Srinivasan S, Morgan MT, Fiedler TL, et al. Metabolic signatures of bacterial vaginosis.
MBio 2015; 6 (2):e00204–e00215.
32. Fettweis JM, Serrano MG, Brooks JP, et al. The vaginal microbiome and preterm birth.
Nat Med 2019; 25 (6):1012–1021.
33. Integrative HMP (iHMP) Research Network Consortium The Integrative Human Microbiome Project: dynamic analysis of microbiome-host omics profiles during periods of human health and disease.
Cell Host Microbe 2014; 16 (3):276–289.
34. Fettweis JM, Alves JP, Borzelleca JF, et al. The vaginal microbiome: disease, genetics and the environment.
Nat Preced 2010; (713):
35. Fettweis JM, Brooks JP, Serrano MG, et al. Differences in vaginal microbiome in African American women versus women of European ancestry.
Microbiol Read Engl 2014; 160 (Pt 10):2272–2282.
36. Meng C, Kuster B, Culhane AC, Gholami AM. A multivariate approach to the integration of multi-omics datasets.
BMC Bioinformatics 2014; 15 (1):162.
37. Chessel D, Hanafi M. Analyses de la co-inertie de K nuages de points.
Rev Stat Appliquée 1996; 44 (2):35–60.
38. Saito S, Nakashima A, Shima T, Ito M. Th1/Th2/Th17 and regulatory T-cell paradigm in pregnancy.
Am J Reprod Immunol N Y N 2010; 63 (6):601–610.
39. Agrawal V, Hirsch E. Intrauterine infection and preterm labor.
Semin Fetal Neonatal Med 2012; 17 (1):12–19.
40. Wennerholm UB, Holm B, Mattsby-Baltzer I, et al. Interleukin-1alpha, interleukin-6 and interleukin-8 in cervico/vaginal secretion for screening of preterm birth in twin gestation.
Acta Obstet Gynecol Scand 1998; 77 (5):508–514.
41. MacIntyre DA, Chandiramani M, Lee YS, et al. The vaginal microbiome during pregnancy and the postpartum period in a European population.
Sci Rep 2015; 5:8988.
42. Walther-António MRS, Jeraldo P, Berg Miller ME, et al. Pregnancy's stronghold on the vaginal microbiome.
PloS One 2014; 9 (6):e98514.
43. Romero R, Hassan SS, Gajer P, et al. The composition and stability of the vaginal microbiota of normal pregnant women is different from that of non-pregnant women.
Microbiome 2014; 2 (1):4.
44. DiGiulio DB, Callahan BJ, McMurdie PJ, et al. Temporal and spatial variation of the human microbiota during pregnancy.
Proc Natl Acad Sci U S A 2015; 112 (35):11060–11065.
45. Huang Y-E, Wang Y, He Y, et al. Homogeneity of the vaginal microbiome at the cervix, posterior fornix, and vaginal canal in pregnant Chinese women.
Microb Ecol 2015; 69 (2):407–414.
46. Hyman RW, Fukushima M, Jiang H, et al. Diversity of the vaginal microbiome correlates with preterm birth.
Reprod Sci 2014; 21 (1):32–40.
47. Elovitz MA, Gajer P, Riis V, et al. Cervicovaginal microbiota and local immune response modulate the risk of spontaneous preterm delivery.
Nat Commun 2019; 10 (1):1305.
48. Petrova MI, Lievens E, Malik S, Imholz N, Lebeer S. Lactobacillus species as biomarkers and agents that can promote various aspects of vaginal health.
Front Physiol 2015; 6:81.
49. Ravel J, Gajer P, Abdo Z, et al. Vaginal microbiome of reproductive-age women.
Proc Natl Acad Sci U S A 2011; 108: (Suppl 1): 4680–4687.
50. Gajer P, Brotman RM, Bai G, et al. Temporal dynamics of the human vaginal microbiota.
Sci Transl Med 2012; 4 (132):132ra52.
51. Cift T, Uludag S, Aydin Y, Benian A. Effects of amniotic and maternal CD-146, TGF-β1, IL-12, IL-18 and IFN-γ, on adverse pregnancy outcome.
J Matern-Fetal Neonatal Med 2013; 26 (1):21–25.
52. Comba C, Bastu E, Dural O, et al. Role of inflammatory mediators in patients with recurrent pregnancy loss.
Fertil Steril 2015; 104 (6):1467–1474.e1.
53. Baeten DLP, Kuchroo VK. How cytokine networks fuel inflammation: Interleukin-17 and a tale of two autoimmune diseases.
Nat Med 2013; 19 (7):824–825.
54. Andraweera PH, Dekker GA, Roberts CT. The vascular endothelial growth factor family in adverse pregnancy outcomes.
Hum Reprod Update 2012; 18 (4):436–457.
55. Reinders MEJ, Sho M, Izawa A, et al. Proinflammatory functions of vascular endothelial growth factor in alloimmunity.
J Clin Invest 2003; 112 (11):1655–1665.
56. Mitchell C, Marrazzo J. Bacterial vaginosis and the cervicovaginal immune response.
Am J Reprod Immunol N Y N 2014; 71 (6):555–563.
57. Cherpes TL, Marrazzo JM, Cosentino LA, Meyn LA, Murray PJ, Hillier SL. Hormonal contraceptive use modulates the local inflammatory response to bacterial vaginosis.
Sex Transm Infect 2008; 84 (1):57–61.
58. Masson L, Mlisana K, Little F, et al. Defining genital tract cytokine signatures of sexually transmitted infections and bacterial vaginosis in women at high risk of HIV infection: a cross-sectional study.
Sex Transm Infect 2014; 90 (8):580–587.
59. Borgdorff H, Gautam R, Armstrong SD, et al. Cervicovaginal microbiome dysbiosis is associated with proteome changes related to alterations of the cervicovaginal mucosal barrier.
Mucosal Immunol 2016; 9 (3):621–633.
60. Anahtar MN, Byrne EH, Doherty KE, et al. Cervicovaginal bacteria are a major modulator of host inflammatory responses in the female genital tract.
Immunity 2015; 42 (5):965–976.
61. Immaculate Mbongo L, Yamunah Devi A, Zain S, Omar SZ, Mohamed Z. Protein profiling of women with spontaneous preterm birth.
Pharmacology 2015; 96 (1–2):44–48.
62. Flores-Díaz M, Monturiol-Gross L, Naylor C, Alape-Girón A, Flieger A. Bacterial sphingomyelinases and phospholipases as virulence factors.
Microbiol Mol Biol Rev 2016; 80 (3):597–628.
63. Hemalatha R, Ramalaxmi BA, KrishnaSwetha G, et al. Cervicovaginal inflammatory cytokines and sphingomyelinase in women with and without bacterial vaginosis.
Am J Med Sci 2012; 344 (1):35–39.
64. Heung LJ, Luberto C, Del Poeta M. Role of sphingolipids in microbial pathogenesis.
Infect Immun 2006; 74 (1):28–39.
65. Serrano MG, Parikh HI, Brooks JP, et al. Racioethnic diversity in the dynamics of the vaginal microbiome during pregnancy.
Nat Med 2019; 25 (6):1001–1011.
66. Schloss PD, Gevers D, Westcott SL. Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies.
PloS One 2011; 6 (12):e27310.
67. Brooks JP, Buck GA, Chen G, et al. Changes in vaginal community state types reflect major shifts in the microbiome.
Microb Ecol Health Dis 2017; 28 (1):1303265.
68. Brooks JP, Edwards DJ, Harwich MD, et al. The truth about metagenomics: quantifying and counteracting bias in 16S rRNA studies.
BMC Microbiol 2015; 15 (1):66.
69. Fadrosh DW, Ma B, Gajer P, et al. An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform.
Microbiome 2014; 2 (1):6.
70. Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform.
Appl Environ Microbiol 2013; 79 (17):5112–5120.
71. Bartram AK, Lynch MDJ, Stearns JC, Moreno-Hagelsieb G, Neufeld JD. Generation of multimillion-sequence 16S rRNA gene libraries from complex microbial communities by assembling paired-end illumina reads.
Appl Environ Microbiol 2011; 77 (11):3846–3852.
72. Parikh HI, Koparde VN, Bradley SP, Buck GA, Sheth NU. MeFiT: merging and filtering tool for illumina paired-end reads for 16S rRNA amplicon sequencing.
BMC Bioinformatics 2016; 17 (1):491.
73. Fettweis JM, Serrano MG, Sheth NU, et al. Species-level classification of the vaginal microbiome.
BMC Genomics 2012; 13: (Suppl 8): S17.
74. Edgar RC. Search and clustering orders of magnitude faster than BLAST.
Bioinforma Oxf Engl 2010; 26 (19):2460–2461.
75. Simanshu DK, Kamlekar RK, Wijesinghe DS, et al. Non-vesicular trafficking by a ceramide-1-phosphate transfer protein regulates eicosanoids.
Nature 2013; 500 (7463):463–467.
76. Bradford MM. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding.
Anal Biochem 1976; 72 (1-2):248–254.
77. Bersanelli M, Mosca E, Remondini D, et al. Methods for the integration of multi-omics data: mathematical aspects.
BMC Bioinformatics 2016; 17 (2):S15.
78. Meng C, Zeleznik OA, Thallinger GG, Kuster B, Gholami AM, Culhane AC. Dimension reduction techniques for the integrative analysis of multi-omics data.
Brief Bioinform 2016; 17 (4):628–641.
79. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2015.
http://www.R-project.org/.
80. Thioulouse J, Chessel D, Dole’dec S, Olivier J-M. ADE-4: a multivariate analysis and graphical display software.
Stat Comput 1997; 7 (1):75–83.
81. Segata N, Izard J, Waldron L, et al. Metagenomic biomarker discovery and explanation.
Genome Biol 2011; 12 (6):R60.
