Replication and refinement of a vaginal microbial signature of preterm birth in two racially distinct cohorts of US women

Study Population and Sampling Procedure.

Two cohorts of pregnant women were enrolled using the same study design, procedures, and methods at separate sites within the United States between 2013 and 2015, and analyzed under a case-control design. One cohort, the “Stanford” cohort, was enrolled from women presenting to the obstetrical clinics of the Lucille Packard Children’s Hospital at Stanford University for prenatal care. The underlying population is predominantly Caucasian and Asian and has a low risk for PTB (<10%); the Stanford cohort analyzed here was selected to enrich for PTB outcomes. Another cohort, the “UAB” cohort, was enrolled from women referred to the UAB for intramuscular 17α-hydroxyprogesterone caproate therapy due to a history of prior preterm delivery. For both cohorts, participants were included in the analysis if they met the following inclusion criteria: 18–45 y of age, singleton gestation, not immunosuppressed, ability to perform the study procedures, and ability to provide informed consent. Both cohorts were enrolled using the same study protocol, but each cohort was analyzed in a separate case-control analysis (using identical bioinformatics approaches). Cases were defined as women who delivered preterm (<37 wk of gestation); controls were women who delivered at term. The analyzed Stanford cohort consisted of 39 women, and the UAB cohort comprised 96 women; demographic and clinical characteristics of each cohort, stratified by enrollment site and by delivery outcome, are presented in Table S1. Features of the 50 women who delivered preterm are presented in Table S2.

Study participants self-collected vaginal swabs weekly by using sterile Catch-All Sample Collection Swabs (Epicentre Biotechnologies) to obtain material from the midvaginal wall. All clinical specimens were placed immediately after collection at −20 °C until transport to the laboratory for storage at −80 °C until processing. Information related to socioeconomic status, diet, medical conditions and symptoms, medication and dietary supplement use, and stress were captured in a detailed questionnaire completed by each study participant upon enrollment; any significant changes in these parameters were documented in a follow-up questionnaire completed at periodic clinical visits.

DNA Extraction, 16S rRNA Gene Amplification, and Amplicon Sequencing.

Whole genomic DNA was extracted from each clinical specimen by means of the PowerSoil DNA isolation kit (MO BIO Laboratories) according to the manufacturer’s protocol except for the inclusion of a 10-min incubation at 65 °C immediately after the addition of solution C1. The V4 hypervariable region of the 16S rRNA gene was amplified by PCR. The forward PCR primer (5′ AAT GAT ACG GCG ACC ACC GAG ATC TAC ACG CTN NNN NNN NNN NNT ATG GTA ATT GTG TGY CAG CMG CCG CGG TAA 3′) was a 74-nucleotide (nt) fusion primer consisting of the 32-nt Illumina adapter (designated by bold), a unique 12-nt barcode to label each amplicon (designated by the N’s), a 9-nt forward primer pad, a 2-nt linker (GT), and the 19-nt broad-range bacterial primer 515F (designated by underlining). The 56-nt reverse primer (5′ CAA GCA GAA GAC GGC ATA CGA GAT AGT CAG CCA GCC GGA CTA CNV GGG TWT CTA AT 3′) consisted of the 24-nt Illumina adapter (designated by bold), a 10-nt reverse primer pad, a 2-nt reverse primer linker (CC), and the 20-nt broad-range bacterial primer 806R (designated by underlining).

Triplicate 25-µL PCRs were carried out by using 1× HotMasterMix (5 PRIME), 0.4 µM concentrations of each commercially synthesized primer, and 3 µL of prepared DNA template. Thermal cycling conditions consisted of an initial denaturing step of 94 °C for 3 min, followed by 30 cycles of 94 °C for 45 s, 50 °C for 60 s, and 72 °C for 90 s, with a final extension step of 72 °C for 10 min. Upon completion of the PCRs, the corresponding triplicate reaction mixtures were pooled and purified by using the Ultra-clean-htp 96-well PCR clean-up kit (Mo Bio Laboratories) according to the manufacturer’s protocol. DNA concentrations from each triplicate pool were quantified using the Quant-iT High-Sensitivity dsDNA Assay Kit (Invitrogen) and combined in equimolar ratios into a single tube. The resulting amplicon mixture was concentrated by ethanol precipitation and resuspended in 100 µL of molecular biology-grade water (Life Technologies). The resuspended amplicon mixture was gel purified and recovered using a QIAquick gel extraction kit (Qiagen). Recovered PCR products were sequenced on an Illumina HiSeq 2500 instrument (Illumina) at the W. M. Keck Center for Comparative Functional Genomics at the University of Illinois, Urbana–Champaign, IL.

Bioinformatics.

Bioinformatics processing largely followed the DADA2 Workflow for Big Data (benjjneb.github.io/dada2/bigdata_paired.html). Forward/reverse read pairs were trimmed and filtered, with forward reads truncated at 245 nt and reverse reads at 235 nt, no ambiguous bases allowed, and each read required to have less than two expected errors based on their quality scores. The relationship between quality scores and error rates was estimated for each sequencing run to reduce batch effects arising from run-to-run variability. ASVs were independently inferred from the forward and reverse of each sample using the run-specific error rates, and then read pairs were merged. Chimeras were identified in each sample, and ASVs were removed if identified as chimeric in a sufficient fraction of the samples in which they were present.

Taxonomic assignment was performed against the Silva v123 database using the implementation of the RDP naive Bayesian classifier available in the dada2 R package (45, 46). Lactobacillus species were assigned by hand via BLAST against sequences from cultured Lactobacillus strains.

Gardnerella Genome Analysis.

The multiple genome alignment for the phylogenetic tree was built with progressiveMauve (47). Briefly, progressiveMauve finds initial local alignments used as seeds, which are progressively extended to global alignments (the pangenome) using a guide tree. The output of progressiveMauve, an extended multifasta file (.xmfa format), was converted to fasta format using a publicly available script (https://github.com/kjolley/seq_scripts/blob/master/xmfa2fasta.pl). The multiple genome alignment was used to construct the phylogenetic tree using FastTree Version 2.1.8 with parameters –nt –gtr –gamma. Orthologous genes were obtained from the Mauve alignment based on at least 70% alignment coverage.

G. vaginalis genomes that contained a full-length 16S rRNA gene were classified as sequence variants G1 (seven genomes), G2 (nine genomes), and G3 (one genome) by exact matching, and were used for gene orthology analysis (Figs. S1 and S2). Presence or absence of genes unique to each strain were not evaluated because only 3 of the 34 Gardnerella genomes are closed. From a total of 3,750 genes predicted in all 17 strains, 1,448 were orthologous genes present in at least half of the genomes in each variant group. Hierarchical clustering on the Euclidean distances of the presence/absence of the 1,448 orthologous genes was done with Cluster 3.0 and Average Linkage as clustering method (48), and visualized with Java Treeview (49). The relative contribution of genes in functional categories was estimated from the number of genes in each functional category divided by the total number of orthologous genes in each cluster.

Accession nos. for Gardnerella vaginalis genomes used in this study are as follows: SAMN02436832, SAMN02472074, SAMN02472073, SAMN02393773, SAMN02393775, SAMN00138210, SAMN02393781, SAMN02393782, SAMN02393783, SAMN02393779, SAMN02393784, SAMN02393777, SAMN02393778, SAMN02393774, SAMN02472072, SAMN02471014, SAMN02436904, SAMN02436712, SAMN02436831, SAMN02436711, SAMN02436912, SAMN02436773, SAMN02436710, SAMN02436911, SAMN02436830, SAMN02436772, SAMN02436771, SAMN02436910, SAMN02436829, SAMN02436909, CP001849, CP002104, CP002725, SAMN02472071.

Low Frequency Variants and Contamination.

Our focus here was on testing previously hypothesized PTB–microbiota associations, and then refining those associations with higher-resolution methods. However, we additionally performed an exploratory analysis of all possible PTB-genus associations to generate new hypotheses while controlling the false discovery rate to be less than a common standard for reporting new associations (FDR < 0.1).

One of the most notable features of this analysis, and one not usually highlighted, is that we believe many of these newly detected associations result from contamination rather than a legitimate biological signal. Of particular note are the Yersinia and Tumebacillus genera, which had prevalence in vaginal and control samples that was strongly correlated across runs (Fig. 4B), suggesting these taxa are contaminants in the vaginal samples. We did not detect such strong run-to-run correlations for Abiotrophia (present in 11 PTB, one term) or Stomatobaculum (16 PTB, four terms), nor when considering other known potential batch effects such as technician or cohort. The correlation coefficients for all of the most significant associations are reported in Table S2.

We are presenting these results without attempting to remove likely contaminants, in part to demonstrate the critical impediment that contaminants present to fulfilling the potential for metabarcoding methods to probe low frequencies. In principle, metabarcoding methods can detect variation at frequencies as low as one over the library size, which here is about a frequency of 10⁻⁵, and when averaging over a gestation can be as low as 10⁻⁶. However, the prevalence of contaminants at study-wide frequencies of up to 10⁻⁴ made identifying real biological patterns at those low frequencies very difficult. This is more problematic in lower biomass samples, such as vaginal swabs, and represents a serious barrier in identifying rare but highly impactful bacterial taxa.

High correlations between prevalence in negative controls and prevalence in vaginal samples of many of the PTB-associated taxa reported here strongly suggests that contaminants are driving many of the potential associations. The problematic taxa are all relatively rare (<5 × 10⁻⁴ frequency), and we do not expect contamination to materially affect inference among more abundant taxa. Unfortunately, identifying contaminating sequences is not trivial; 94/100 of the most abundant genera in our study appear in negative controls because of cross-contamination, eliminating the possibility of simply removing anything that appears in a negative control. The size of our study allows us to detect some contaminants by looking at prevalence correlations between the natural batches of sequencing runs, but even so there are some taxa for which it is difficult to make a determination.

Run-specific contaminants should minimally interfere with exploratory testing (beyond decreasing power by increasing the number of hypotheses tested) if samples are randomly distributed among sequencing runs. Unfortunately, the samples in this study were not completely randomized among the sequencing runs, in part because of the need to give sequencing precedence to the samples from a nonrandom subset of the subjects that were also included in a second collaborative study. This likely contributed to the number of spurious associations we found, and our results highlight the importance of robust randomization procedures for exploratory studies, especially those which span multiple sequencing runs and that intend to test low-abundance taxa.

Comparison with Results in Kindinger et al.

A recent study by Kindinger et al. (37) also examined the relationship between the vaginal microbiota and PTB in a similarly sized cohort (n = 161) of predominantly Caucasian and Asian women with a prior history of preterm birth. Using a CST-based analysis, they corroborated an association between a reduction in PTB and an increased relative abundance of L. crispatus during pregnancy. They also found a weaker or absent association between the vaginal microbiota and PTB in women of black race. However, their study differed from ours in finding a significant positive association between L. iners-dominated communities (CST3) and PTB, and no significant association between diverse BV-like communities (CST4) and PTB. As a result, the authors argued that L. iners might promote PTB, rather than members of the diverse BV-associated community as has been previously suggested and as was supported by our study. We believe the differences in our findings may reflect methodological differences in our 16S rRNA gene sequencing protocols and analysis methods, rather than biological differences.

Kindinger et al. used a V1-V3 primer set that is sensitive to Lactobacillus sequences, but may be insensitive to some other important taxa. In particular, Gardnerella is essentially absent (<0.01%) in their amplicon sequencing data, while it accounts for almost 5% of our study-wide sequencing reads. We do not believe this reflects structural differences in the vaginal communities between our cohorts given the widespread observation of significant Gardnerella fractions in previous population studies of the pregnant vaginal microbiota (11, 12). Instead, we suspect that the three mismatches in their forward primer prevented amplification of the Gardnerella 16S rRNA gene and, thus, Gardnerella remained essentially absent.

If it is the case that the primer set used by Kindinger et al. is systematically less sensitive to non-Lactobacillus species than ours, the divergences between the studies may yet reflect a common biological reality. Our study showed clearly that L. iners and BV-associated taxa like Gardnerella often exist at comparable frequencies. An amplicon sequencing protocol that poorly amplified BV-associated taxa would then cause amplicon sequencing data from BV-like communities to contain a preponderance of reads from the well-amplified L. iners populations with which they coexist. This effect would be compounded by assigning communities to CSTs, which largely tracks the identity of the most frequent taxa in amplicon reads. Thus, even if BV-like communities were promoting PTB, as we believe our findings suggest, the classification of BV-like communities (CST4) as L. iners-dominated communities (CST3) based on amplicon data depleted of Gardnerella sequences could create the appearance of an association between communities classified as CST3 and PTB that would otherwise be absent.

The differences between the specifics of the positive associations detected in these two studies should not diminish the areas of agreement: Both studies found a significant protective role for L. crispatus, a significant positive association between states of the vaginal microbiota not dominated by L. crispatus and PTB in predominantly Caucasian populations, and weaker or absent associations in African-American women. Together, we believe these studies have strengthened the evidence for a role of the vaginal microbiota in the pathophysiology of PTB, and point toward L. iners and G. vaginalis as taxa of particular importance for further study. Other taxa that often cooccur with G. vaginalis, such as the genera Prevotella, Dialister, Mobiluncus, and Atopobium, among others, also merit additional investigation.