Factors associated with the composition and diversity of the cervical microbiota of reproductive-age Black South African women: a retrospective cross-sectional study

Ethics statement

This study was approved by the Human Research Ethics Committee of the University of Cape Town, South Africa (references 258/2006 and 580/2014). All participants provided written informed consent to participate in the study and use of their stored samples for future studies.

Study population and study design

This was a retrospective cross-sectional study based on data and baseline cervical DNA samples from the HPV Couples Cohort Study (Mbulawa et al., 2009) that examined the transmission of genital HPV among Black heterosexual couples in Gugulethu, Cape Town, South Africa. Details of enrollment, recruitment and sample collection for the HPV transmission study have been described previously (Mbulawa et al., 2009).

In brief, speculum examination was performed and excess mucus around the cervical area was cleared using a filamented swab. This was followed by collection of two cytobrush samples from the cervix. The first sample was for Papanicolaou (Pap) smear or T cell assay. It was collected by inserting the cytobrush into the mouth of the cervix and rotating the cytobrush at 360° thrice. For the Pap smear, the sample was smeared immediately onto the frosted glass slide, quickly fixed using Cytofix spray and then stained with Pap stain. BV was identified on Pap smears by using the Bethesda criteria for reporting cervical/vaginal cytologic diagnoses (Kurman & Solomon, 1994). Smears showing clue cells with coccobacilli (mostly Gardnerella vaginalis) and/or any shifts in bacterial flora suggestive of BV (noticeable absence of lactobacilli) on wet microscopy were considered as having findings suggestive of BV. All smears read by cytotechnologists were reviewed. The second sample from the cervix was for HPV genotyping and herpes simplex virus (HSV) testing (and subsequent analyses such as characterization of microbiome). It was collected by inserting a second cytobrush (Digene cervical sampler) into the cervix and rotating it thrice (360°) inside the mouth of the cervix. This sample was then stored in Digene specimen transport medium (Digene Corporation, Gaithersburg, MD, USA) at −80 °C until nucleic acid extraction.

Nucleic acids were extracted from the cervical samples as previously described (2009) (Mbulawa et al., 2009) using the MagNA Pure Compact System and the MagNA Pure Compact Nucleic Acid Isolation kit (Roche Molecular Diagnostics, Mannheim, Germany). HPV typing was also performed as previously documented (Mbulawa et al., 2009) using the Roche Linear Array HPV genotyping test (Roche Molecular Diagnostics, Mannheim, Germany) that detects 37 HPV genotypes. These include 12 oncogenic high-risk, 8 probable oncogenic high-risk, and 17 non-oncogenic low-risk HPV types as listed elsewhere (Mbulawa et al., 2009). Only samples with positive human beta (β)-globin (a housekeeping gene) hybridization results (a measure of sample adequacy) were included in this study. Roche Linear Array HPV genotyping test measures sample adequacy by relying on two endogenous β-globin positive controls (high and low) run concurrently with samples. The primers targeting β-globin are different from those that target the HPV genome (polymorphic L1 region). A valid genotyping result (negative or positive for HPV genotype) is one where both the β-globin probe lines are positive. A result (negative or positive for at least one HPV genotype) is considered invalid if the sample is negative for one or both β-globin control(s). This is suggestive of inadequate cellular material, poor storage and processing (extraction), presence of PCR inhibitors and/or completion with a high titer HPV target.

To be eligible for the present study, only cervical DNA specimens from human immunodeficiency virus (HIV)-seronegative women aged 18-44 years were considered. These samples had to have information on the HPV status and sufficient volume for microbiota analysis (≥15.0 µl of the extracted DNA). Exclusion criteria included being a woman aged <18 or >44 years, self-reported menstruation or pregnancy at the time of sampling, and HIV-seropositivity. Participants’ metadata including demographics, sexual history, smoking, contraceptive use, and clinical characteristics were abstracted from the HPV Couples Cohort Study (Mbulawa et al., 2009).

Bacterial V4 hypervariable region (16S ribosomal rRNA) library preparation and sequencing

The hypervariable V4 region of the 16S ribosomal rRNA (rRNA) gene was amplified using the universal polymerase chain reaction (PCR) primers 515f (5′-GTGCCAGCMGCCGCG GTAA-3′) and 806r (5′-GGACTACHVGGGTWTCTAAT-3′) (Caporaso et al., 2011). Each PCR contained 1x Ex Taq buffer (Takara Bio Inc., Japan), 0.025 U Ex Taq polymerase, 0.8 mM deoxynucleotide triphosphate (dNTP) mixture, 0.56 mg/ml bovine serum albumin (BSA), 400 nM each primer and 100 ng template. Each sample (and no template PCR control, i.e., nuclease free water) was amplified in 3 replicate reactions. PCR conditions were 98 °C for 2 min, followed by 30 cycles of 98 °C for 20 s, 50 °C for 30 s and 72 °C for 45 s, and a final elongation step at 72 °C for 10 min. The triplicate samples were pooled and purified using the Agencourt AMPure XP system (Beckman Coulter, Germany) according to the manufacturer’s instructions. Amplicon sizes were confirmed by electrophoresis on 1.5% Tris Borate EDTA (TBE) agarose gels and imaging with the ultraviolet transilluminator (UVT) GelDoc-It™ system. The amplicons were quantified using the Quant-iT^® PicoGreen dsDNA assay (Thermo Fisher Scientific, USA) with FLUOstar OPTIMA Multi-Mode Micro Plate Reader (BMG Labtech, Germany). Sequencing libraries were prepared using the KAPA Library Preparation kit and barcoded Adaptor kits for Ion Torrent™ platforms (KAPA Biosystems, Wilmington, MA, USA). Purified barcoded amplicons were pooled (24 samples per pool) in equimolar amounts and the final library sizes and concentrations assessed on a Bioanalyzer High Sensitivity Chip (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed on the Ion Torrent Personal Genome Machine (PGM) (Life Technologies, Beverly, MA, USA) at the Central Analytical Facilities (CAF) at Stellenbosch University (Stellenbosch, South Africa).

Bacterial V4 hypervariable region (16S rRNA) data analysis using bioinformatics tools

The qualities of the raw sequenced reads were visually inspected using FastQC v0.11.2 (Andrews, 2010). Quantitative Insights Into Microbial Ecology (QIIME) v1.8.0 (Caporaso et al., 2010b) with imported UPARSE (usearch7.0.1090) (Edgar, 2013), was used to analyze and interpret the nucleotide sequence data from the cervical microbiota. In the initial sequence pass, reads were quality-filtered and demultiplexed in QIIME using the user-defined parameters in Table S1. Briefly, reads with lengths outside the 200-400-bp range, with a quality score of <25 (sliding window 50) and without barcodes or with any mismatches in the barcode sequences were discarded. A second quality filter was performed in UPARSE with user-defined parameters (Table S1). Sequences were dereplicated followed by abundance sorting and discarding singletons. Operational taxonomic unit (OTU) clustering was performed by UPARSE-OTU method that uses a greedy clustering algorithm, with binning of reads with 97% pairwise identity. This step was performed simultaneously with representative sequence picking and de novo chimera filtering. Representative sequences from each unique OTU cluster were picked using abundance algorithm. Additional chimeras were removed by UCHIME algorithm (Edgar et al., 2011). Taxonomy was assigned using the Ribosomal Database Project (RDP) Naïve Bayesian Classifier (Wang et al., 2007), with the Greengenes database (gg13_8 Release) (De Santis et al., 2006). Phylogeny was inferred by aligning representative sequences to Greengenes core set using Python Nearest Alignment Space Termination (PyNAST) (Caporaso et al., 2010a). A phylogenetic tree was built using FastTree. Other parameters used for our analyses are defined in Table S1.

Diversity, rarefaction, and sample ordinations were computed in QIIME. Multiple rarefactions at different sequencing depths were performed, and rarefactions (collector’s) curves plotted to evaluate the completeness of the sampling efforts. Alpha diversity was computed by chao1, observed_species, Shannon, Simpson, and PD_whole_tree metrics. Beta diversity was estimated using weighted and unweighted UniFrac distances, and Bray-Curtis dissimilarity metric. The strength and statistical significance of sample clustering (beta diversity) was computed using permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001), with 999 permutations. Other biodiversity metrics including Dominance and Shannon Equitability indices were calculated using an in-house script in RStudio v1.1.447 (RStudio Team, 2016). An all-by-all pairwise distance matrix of UniFrac distances were generated and used to hierarchically cluster and ordinate samples. The ordinations were performed using Principal Coordinate Analysis (PCoA).

Identification and comparison of community state types

Hierarchical clustering and heatmap generation were performed in RStudio v1.1.447 (RStudio Team, 2016). Hierarchical clustering on Bray-Curtis dissimilarity (Vegan package v2.5-5 (Philip, 2003)) was done using the average neighbour algorithm. The heatmap was generated using the heatmap3 v1.1.6 package (Zhao et al., 2014).

Correlational analyses of 60 cervical bacterial OTUs

To assess bacterial positive and negative relationships, Spearman’s correlation of log₂-transfomed counts of all OTU pairs were tested on metagenomeSeq v1.12.1 (Paulson et al., 2014; Paulson et al., 2016) in RStudio v1.1.447 (RStudio Team, 2016), and the topmost OTUs (n = 60, 28.0%) with the greatest variance displayed on a correlogram.

Statistical methods

Statistical analyses were carried out using GraphPad Prism v6.01 (San Diego, CA, USA). Mann–Whitney unpaired nonparametric and Chi-square/Fisher’s exact tests (with two-tailed p-value) were used to examine the association of continuous and categorical variables with CSTs, as appropriate. The alpha diversity metrics of the three CSTs were compared by Kruskal-Wallis test. Two-group comparison between the alpha diversity of CST, HPV, and BV groups was computed by Mann–Whitney unpaired nonparametric test.

Study cohort baseline characteristics

The demographic, sexual and smoking, behavioural and clinical information of the 62 heterosexual Black South African women included in this study are summarized in Table 1. All the women were sexually active, with the majority (72%) of hormonal contraceptives users being on Depo-Provera. HPV and high-risk HPV infections were detected in 37.1% and 29.0% of the women, correspondingly. Thirteen of the women (22.0%) had abnormal cervical cytology. A few women had experienced vaginal discharge (16.1%) and genital ulceration (3.2%) in the last six months. Twenty two women (35.5%) had findings suggestive of BV. A majority of the women (79.0%) had never smoked cigarettes.

Taxonomic composition of the cervical microbiota

A total of 1,392,562 high-quality non-spurious sequencing reads from 62 samples were included in the final analysis with a median of 16,453 reads per sample (range: 5,343–235,897 per sample). The Ion Torrent PGM raw sequence data and metadata have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA472137 (SRP148486) with accession numbers (SRX4103412 –SRX4103416).

A heatmap of the relative abundances of the predominant genera (>0.5% relative abundance) and their respective phyla identified in the cervical microbiota of the 62 women is shown in Fig. S1. Above the 0.5% threshold, only 9 phyla with 43 genera were observed. Lactobacillus (46.0% mean relative abundance and 100.0% prevalence), Gardnerella (19.5% and 93.5% (58/62)), Prevotella (9.7% and 95.2% (59/62)), and Sneathia (9.5% and 87.1% (54/62)) were the most predominant genera in the phyla Firmicutes, Actinobacteria, Bacteriodetes, and Fusobacteria, respectively. Other bacteria such as Dialister (1.2%), Atopobium (0.5%), and Clostridium (0.1%) commonly found in women with BV (Srinivasan et al., 2012) were not abundant.

A total of 221 unique OTUs (potential species) were identified in the 62 cervical microbiota and ranged from 16-104 OTUs per sample. The 30 most abundant OTUs represented 97% of all the reads. The most abundant OTU was classified as L. iners, representing 40.6% of all the reads. The prevalence and abundance of this OTU together with those of the other most abundant OTUs are shown in Table S2.

Prevotella was the most diverse genus with 26 different OTUs identified. Six and eight different OTUs belonging to the genera Gardnerella and Lactobacillus, respectively, were detected. Apart from L. iners, the other Lactobacillus OTUs included Lactobacillus coleohominis, Lactobacillus mucosae, Lactobacillus ruminis, and four unclassified Lactobacillus spp. Species names could not be assigned to the four Lactobacillus OTUs due to insufficient taxonomic discrimination by the V4 hypervariable region of the 16S rRNA gene. Forty seven (75.8%) of the cervical microbiota had at least two Lactobacillus spp., but at unequal abundances.

Characterization of cervical community state types

Hierarchical clustering of the cervical microbiota based on the type and relative abundances of the bacterial taxa identified three distinct community state types, CSTs I–III (Fig. 1). CST I was dominated by L. iners and found in 24 women (38.7%). CST II was dominated by an unclassified Lactobacillus (Lactobacillus.4) and present in only three women (4.8%). CST III was the most common CST occurring in 35 women (56.5%). This CST was characterized by a diverse and complex array of facultative and strictly anaerobic BV-associated bacteria (Gardnerella, Prevotella, Sneathia, Shuttleworthia, Clostridium, Atopobium, Dialister, and a consortium of low-abundant bacteria) and very low numbers of Lactobacillus, including L. iners and three unclassified Lactobacillus (Lactobacillus.1, Lactobacillus.3, and Lactobacillus.4). While there was a continuum of relative abundances of the bacterial taxa in CST III, four sub-clusters were evident (Fig. 1), three sub-clusters were dominated by Shuttleworthia (n = 4 women, 11.4%), Gardnerella (n = 8, 22.9%), and Sneathia ( n = 5, 14.3%), and one with mixed taxa (n = 18, 51.4%). Next, the metadata for the women in each CST were compared.

Comparison of the community state types by participants’ metadata

The demographic, sexual, smoking, and clinical characteristics of the women assigned to each of the three CSTs are shown in Table 2. The metadata for women in CST I and CST III were compared, while CST II was excluded from statistical comparisons due to the small sample size.

A significantly greater number of women with cervical microbiota from CST I reported hormonal contraceptive use compared to women with CST III (15/22 (68.2%) versus 9/31 (29%), p = 0.005). Findings suggestive of BV on smears were significantly more frequent (p = 0.001) in women with CST III than CST I. The other participants’ variables, including HPV status, were not significantly different between the women in CST I and III.

Comparison of alpha diversity across CSTs, BV, HPV, and hormonal contraceptive use

Alpha diversity in the cervical microbiota was estimated using a variety of indices, including Simpson, Dominance, Shannon Diversity and Shannon Equitability (Fig. 2). A higher Dominance, Shannon, and Shannon Equitability index value, and lower Simpson index value designates greater alpha diversity. Based on rarefaction curves of alpha diversity metrics (Shannon and Simpson indices, Fig. S2), we chose 5,000 reads per sample as a sufficient subsampling depth to accurately assess microbial diversity.

When grouped by CST (Fig. 2A), CST I (L. iners-dominated group) and CST III (diverse group) were significantly different for all the alpha diversity indices (p < 0.0001). Thus, bacterial diversity in CST III was significantly greater than CST I.

The diversity of the cervical microbiota of women with findings suggestive of BV was significantly greater than that of women without findings of BV (Shannon index: 2.1 (1.8–2.4) versus 0.5 (0.2–1.9), p < 0.0001), Fig. 2B. Alpha diversity was also shown to be significant by Simpson, Dominance, Shannon, and Shannon Equitability metrics.

When grouped by HPV status (Fig. 2C), no significant difference in alpha diversity was observed between the HPV-negative and HPV-positive groups, with Shannon index of 1.9 (0.3–2.2) and 1.0 (0.3–2.2), p = 0.483, respectively. Alpha diversity was also shown to be significant by Simpson, Dominance, Shannon, and Shannon Equitability metrics.

Lastly, hormonal contraceptive users (Shannon index: 0.9 (0.2–1.9)) had a lower microbial diversity than non-users (2.3 (0.6–2.3)), p = 0.025 (Fig. 2D). When stratified by the type of contraceptive used, the diversities based on Shannon index were significantly different (p = 0.028); microbial diversity of Depo-Provera (0.7 (0.2–1.7)) and norethisterone enanthate (0.5 (0.2–1.4)) was significantly lower than that of non-users (2.0 (0.6–2.3)) (Fig. S3).

Comparison of beta diversity across CSTs, BV, HPV, and hormonal contraceptive use

Beta diversity analysis (PCoA of weighted UniFrac distances) of the 62 samples showed that each of the established CST (I-III) represented a highly distinct bacterial community, p = 0.001 (Fig. 3A). This result was supported by the Jackknife replicates that were used to estimate the uncertainty in PCoA plots and hierarchical clustering of the cervical microbiota. The majority of the samples with Lactobacillus-dominated cervical microbiota (CST I and CST II, 25/27, 92.6%) exclusively clustered together in the upper right quadrant (Fig. 3A). Two samples from CST I did not cluster with this group likely due to presence of other bacterial taxa, e.g., Gardnerella, Prevotella, and Aerococcus, in these cervical microbiota. These samples were from women with BV.

Further, beta diversity analysis showed that the clustering of the samples was dependent on the findings suggestive of BV, p = 0.001 (Fig. 3B). The 25 samples that clustered together in the upper right quadrant consisted mostly of women without findings suggestive BV (22/25, 88.0%). Samples from women with findings suggestive of BV were spread over a greater area in the plot due to their high and varying bacterial diversity.

The weighted UniFrac distances of the cervical microbiota showed that there was no apparent influence of HPV infection on beta diversity, p = 0.373 (Fig. 3C). Of the samples that clustered together, 48.0% (12/25) were HPV-positive.

The majority of the samples from women on hormonal contraception (16/25, 64.0%), mostly Depo-Provera (12/18, 75.0%), clustered together in the upper right quadrant, p = 0.019 (Fig. 3D). Whereas 80.0% (4/5) of the samples from women on norethisterone enanthate clustered together with the 75.0% samples from women on Depo-Provera, all the samples from women on oral contraceptives (2/2, 100.0%) did not (p = 0.018).

Co-occurrence and co-exclusion patterns of cervical bacterial OTUs

Pairwise correlations were calculated for all pairs of OTUs identified in the cervical microbiota. A correlation matrix of the pairwise correlations between the 60 most abundant cervical bacterial OTUs is shown in Fig. 4.

From the dendrogram, two major bacterial correlation clusters, Cluster-A (mostly with OTUs classified as Lactobacillus and Streptococcus) and Cluster-B (mostly with BV-associated bacteria) were observed. Each of these clusters had two sub-clusters: Cluster-A1 and Cluster-A2 for Cluster-A, and Cluster-B1 and Cluster-B2 for Cluster-B. OTUs in Cluster-A had an inverse correlation with OTUs in Cluster-B. For example, Lactobacillus OTUs (in Cluster-A) had negative correlations with Gardnerella and Prevotella OTUs (in Cluster-B). OTUs in the same Cluster-A1 had stronger positive correlations with one another than with OTUs in another sub-cluster from the same cluster (e.g., Cluster-A2). There was some overlap in the interaction of bacteria in Cluster-A2 and Cluster-B1. Interactions between these sub-clusters appeared very low to moderate. We noted that strong positive correlations were very common between phylogenetically related bacterial OTUs, e.g., Lactobacillus spp., but the extent of these interactions varied.