tified eight nucleotide positions that showed high variation in the V4–V6 region of G. vaginalis 16S rRNA gene (Figure 1). The variable locations that emerged from this analysis coincided with 511st, 612th, 661st, 835th, 988th, 989th, 990th and 991st nucleotide positions of the 16S rRNA gene from the genome sequence of G. vaginalis strain 409-05 (GenBank accession number: CP001849). None of these positions were associated with homopolymer regions, and nucleotide variations at these locations were also observed in some of the full-length G. vaginalis 16S rRNA gene sequences found in the RDP database. For each sequence in the tag library, we retained nucleotides only from those highly variable nucleotide positions and merged them into eight nucleotide oligomers, and used these oligomers to label individual G. vaginalis ‘oligotypes’. To reduce the probability of including an oligotype containing a nucleotide that may have been randomly generated by a sequencing error, we used only those oligotypes that were present in at least two samples. The resulting 46 oligotypes were used to generate G. vaginalis oligotype profiles for individual samples.
Peaks in entropy indicate nucleotide variation at given locations. While the X-axis of the figure indicates the location of the given column in the full length G. vaginalis 16S rRNA gene sequence, bar at the top superimposes the approximate locations of hyper-variable regions V4 (557–662), V5 (800–861) and V6 (981–1027) on G. vaginalis sequences that were used in this study.
Analyzing correlations among oligotype profiles
We used SciPy, an open-source scientific computation library (available from http://scipy.org/) for Python programming language, to compute Pearson correlation coefficients and p-values in order to explore linear relationships between sexual partners based on their G. vaginalis oligotype profiles. Pearson correlations were computed over the feature vectors that were constructed based on the percent abundance of oligotypes for every sample. The number of reads representing each oligotype was tallied for each sample to generate a 46-dimensional feature vector where each component of the vector reflected the percent abundance of the corresponding oligotype within the given sample. Pearson correlation analysis results are listed in Table 2.
Phylogenetic analysis of oligotypes and UniFrac clustering
Phylogenetic relationships among the oligotypes were assessed with Bayesian inference using MrBayes (version 3.1.2, http://mrbayes.sourceforge.net/) , . Analysis was initiated with random starting trees with representative sequences for each oligotype, and posterior probabilities were determined from two independent runs of one million generations of Markov chain Monte Carlo simulations, from which tree topologies were sampled every 100 generations. After discarding the first 25% of resulting trees, a consensus phylogenetic tree of oligotypes was estimated from remaining generations (Figure S1). The resulting tree was used as a common phylogeny to perform UniFrac analysis . Hierarchical clustering of oligotypes in vaginal (Figure 2), and penile skin and urethra samples (Figure S2) was performed based on distance matrices generated by the unweighted UniFrac analysis. Tree topology of the phylogenetic analysis and clustering results were visualized using the Interactive Tree of Life .
Samples were clustered (clustering significance: p<0.001, UniFrac significance: p = 0.016) based on the UniFrac distance metric. Red, orange and green colors indicate samples from BV, intermediate and normal Nugent categories, respectively.
The generation of oligotype profiles by merging nucleotides from variable locations of G. vaginalis tag sequences that are revealed by Shannon entropy analysis (Figure 1) made it possible to compare samples to each other based on their G. vaginalis oligotype compositions. This analysis showed extensive diversity within G. vaginalis sequences from different samples, as well as significant correlations between the oligotype profiles of many couples. The compositio