The cervical microbiome over 7 years and a comparison of methodologies for its characterization

Abstract

Background: The rapidly expanding field of microbiome studies offers investigators a large choice of methods for each step in the process of determining the microorganisms in a sample. The human cervicovaginal microbiome affects female reproductive health, susceptibility to and natural history of many sexually transmitted infections, including human papillomavirus (HPV). At present, long-term behavior of the cervical microbiome in early sexual life is poorly understood.

Methods: The V6 and V6-V9 regions of the 16S ribosomal RNA gene were amplified from DNA isolated from exfoliated cervical cells. Specimens from 10 women participating in the Natural History Study of HPV in Guanacaste, Costa Rica were sampled successively over a period of 5-7 years. We sequenced amplicons using 3 different platforms (Sanger, Roche 454, and Illumina HiSeq 2000) and analyzed sequences using pipelines based on 3 different classification algorithms (usearch, RDP Classifier, and pplacer).

Results: Usearch and pplacer provided consistent microbiome classifications for all sequencing methods, whereas RDP Classifier deviated significantly when characterizing Illumina reads. Comparing across sequencing platforms indicated 7%-41% of the reads were reclassified, while comparing across software pipelines reclassified up to 32% of the reads. Variability in classification was shown not to be due to a difference in read lengths. Six cervical microbiome community types were observed and are characterized by a predominance of either G. vaginalis or Lactobacillus spp. Over the 5-7 year period, subjects displayed fluctuation between community types. A PERMANOVA analysis on pairwise Kantorovich-Rubinstein distances between the microbiota of all samples yielded an F-test ratio of 2.86 (p<0.01), indicating a significant difference comparing within and between subjects' microbiota.

Conclusions: Amplification and sequencing methods affected the characterization of the microbiome more than classification algorithms. Pplacer and usearch performed consistently with all sequencing methods. The analyses identified 6 community types consistent with those previously reported. The long-term behavior of the cervical microbiome indicated that fluctuations were subject dependent.

Figure 1. Primer design.

Panel A shows the primer design for amplifying the V6–V9 16S rRNA gene region analyzed by cloning and Sanger sequencing and 454 pyrosequencing. Regular (non-encoded) 10 bp barcodes were added to the 5′end of the forward PCR primer. Panel B shows the primers used to amplify the 16S V6 region, analyzed by Illumina sequencing. Hamming barcodes (8 bp in length) and padding sequences were introduced to the 5′ ends of the forward and reverse PCR primers, different for each strand, so that reads from each strand could be distinguished. Note: the reverse primer sequences shown are the actual oligonucleotide sequences used in PCR amplification (i.e., the reverse complement of the 5′–3′ target DNA sequence).

Figure 2. Flowchart of sequencing technologies, methodological pipelines and associated software.

Sequencing files in FASTQ or FASTA, and QUAL formats underwent the following steps shown in the indicated panels: (panel 1) Quality filtering, where short and low quality reads were discarded and chimeric sequences were detected and removed; (panel 2) Read demultiplexing was performed where reads were assigned names according to the clinical sample from which they originated based on each unique barcode; (panel 3) Read identification was performed using (subpanel A) usearch, (subpanel B) RDP Classifier, and (subpanel C) pplacer. For usearch and pplacer, classification involved multiple processing steps and format modifications (panel 4) to allow for direct comparison between methodological configurations. The data standardizing scripts yielded tables containing the counts for each detected genus (rows) and clinical sample (columns). Some taxa appeared multiple times in the initial tables, therefore the counts for these taxa were pooled. Filtering was also applied to discard any counts that constituted <1% of the total sample composition. Taxa that were empty of counts across all samples after this low-pass filtering were discarded. Finally, to allow direct comparison, all nine classification-tables were formatted such that the numbers of rows and columns in each table were equal and contained a union of all taxa and samples.

Figure 3. Community compositions of cervical samples at the genus level as determined by 9 different methodological configurations.

Heat-maps show the log₁₀ (proportional abundance) of each bacterial genus detected in each clinical sample for each methodological combination. In the 3×3 grid of heat-maps, the sequencing method is indicated at the far left and the classification software is indicated at the top of the panels. The cladograms to the left of the genus names indicate the approximate evolutionary relationships between genera.

Figure 4. Microbiome reclassification by different methodological configurations.

Panel A shows a boxplot of the percentage of total reads reclassified as different genera for all samples, for each pairwise comparison between methodologies. Outliers beyond the interquartile range are shown as points. Panel B shows the percentage of total classified reads assigned to each genus. The distributions reflected by the boxplots indicate the variability of the classification percentages between methodological configurations.

Figure 5. A comparison of the Shannon diversity indices for each methodology.

The Shannon diversity index (H’) was calculated based on the genus-level classification tables produced by each combination of sequencing method and software pipeline. The boxplots show the distribution of H’ values across all samples. For a given sequencing method, the Shannon diversity index appears consistent across classification software, except for the Illumina and RDP Classifier combination, where a large increase in apparent diversity occurs.

Figure 6. Categorical community types by squash clustering and prevalence of species.

Reads from the 454 (panel A) and Illuimna (panel B) platforms were classified at the species level by pplacer and guppy, and clustered using squash clustering . The figure shows the distributions of reads between species for each clinical sample as heat-maps, on a logarithmic scale, arranged according to the squash clustering. The tree produced by the clustering algorithm is shown at the top of the heat-map, with community type designations appearing below; the type names are in accord with those proposed by Ravel et al. .

Figure 7. Time courses and distribution of microbiome community types.

Panels A and B show the proportions of samples assigned to each community type using 454 and Illumina, respectively, for the whole study population (10 women) across all time points (5±1). Experimental replicates are excluded. Community types III and IV constitute over half of the cervical microbiome from these women. Panels C and D show the microbiome community types over time, as characterized by 454 and Illumina, respectively, when using pplacer and guppy to classify and cluster the reads.