Characterization of the Vaginal Microbiota among Sexual Risk Behavior Groups of Women with Bacterial Vaginosis

Background The pathogenesis of bacterial vaginosis (BV) remains elusive. BV may be more common among women who have sex with women (WSW). The objective of this study was to use 454 pyrosequencing to investigate the vaginal microbiome of WSW, women who have sex with women and men (WSWM), and women who have sex with men (WSM) with BV to determine if there are differences in organism composition between groups that may inform new hypotheses regarding the pathogenesis of BV. Methods Vaginal swab specimens from eligible women with BV at the Mississippi State Department of Health STD Clinic were used. After DNA extraction, 454 pyrosequencing of PCR-amplified 16S rRNA gene sequences was performed. Sequence data was classified using the Ribosomal Database Program classifer. Complete linkage clustering analysis was performed to compare bacterial community composition among samples. Differences in operational taxonomic units with an abundance of ≥2% between risk behavior groups were determined. Alpha and beta diversity were measured using Shannon’s Index implemented in QIIME and Unifrac analysis, respectively. Results 33 WSW, 35 WSWM, and 44 WSM were included. The vaginal bacterial communities of all women clustered into four taxonomic groups with the dominant taxonomic group in each being Lactobacillus, Lachnospiraceae, Prevotella, and Sneathia. Regarding differences in organism composition between risk behavior groups, the abundance of Atopobium (relative ratio (RR)=0.24; 95%CI 0.11-0.54) and Parvimonas (RR=0.33; 95%CI 0.11-0.93) were significantly lower in WSW than WSM, the abundance of Prevotella was significantly higher in WSW than WSWM (RR=1.77; 95%CI 1.10-2.86), and the abundance of Atopobium (RR=0.41; 95%CI 0.18-0.88) was significantly lower in WSWM than WSM. Overall, WSM had the highest diversity of bacterial taxa. Conclusion The microbiology of BV among women in different risk behavior groups is heterogeneous. WSM in this study had the highest diversity of bacterial taxa. Additional studies are needed to better understand these differences.


Mycoplasma
hominis, Bacteroides spp., Peptostreptococcus spp., Fusobacterium spp., Mobiluncus spp., Atopobium vaginae, etc.) bacteria [11]. It remains controversial as to whether BV results from acquisition of one organism (i.e. Gardnerella vaginalis) as the founder organism leading to the complex microbial community that characterizes BV or whether BV is transmitted by a polymicrobial consortium of micro-organisms [12][13][14]. BV is most commonly diagnosed using a combination of physical exam and laboratory findings (Amsel criteria) [2] although the diagnosis is more rigorously defined by Gram staining of vaginal fluid to determine the Nugent score [15].
Identification of specific bacterial species present in vaginal fluid can be done using either cultivation-dependent or cultivation-independent molecular methods [12,[36][37][38]. Cultivation-independent molecular methods are increasingly being used because of their ability to detect fastidious microorganisms that cannot be cultivated by conventional methods. A landmark study by Fredricks et al used a combination of several broad-range PCR techniques to identify bacterial species present in samples of vaginal fluid from women with BV and without BV [38]. In this study, women with BV had a mean of 12.6 bacterial species per sample compared to women without BV who had a mean of 3.3 bacterial species per sample. Nevertheless, despite the use of both cultivationdependent and cultivation-independent molecular methods, the cause of the complex vaginal flora consistent with BV has remained elusive as Koch's postulates for establishing disease causation have not been fulfilled for any bacterium or group of bacteria associated with BV [39].
Since cultivation-independent molecular methods have become increasingly available, researchers have continued to investigate the vaginal flora of women with and without BV [37,[40][41][42][43]. However, there is little cultivation-independent data investigating the complex vaginal flora among different sexual risk behavior groups of women with BV. Fethers and colleagues have recently shown with species-specific primers and qPCR that Megasphaera type I was 91.3% (95% CI 70.4-98.4) sensitive for BV in WSW but only 72.2% (95% CI 61.2-81.2) sensitive for BV in non-WSW [44]. The authors concluded that BV in WSW has a distinct microbiota compared with other women and acknowledged that additional studies are needed to further explore any potential differences that may be present among these groups of women. The objective of this study was to use 454 pyrosequencing to investigate the vaginal microbiome of WSW, women who have sex with women and men (WSWM), and women who have sex with men (WSM) with BV to determine similarities and differences in organism composition between sexual risk behavior groups that may inform new hypotheses regarding the pathogenesis of this common vaginal infection.  [45] or a Mycoplasma genitalium cohort study between October 2006 and January 2009 (WSM) (data unpublished) were used in this study. All women enrolled in each of these studies were of reproductive age, not pregnant, and reported a history of sexual activity (oral, vaginal, and/or anal) with women only (WSW), with both women and men (WSWM), or with men only (WSM) during the 12 months preceding enrollment.

Study Population and Sample Collection
All women had undergone a pelvic examination when enrolled in their respective study. Two vaginal swabs were collected by a study clinician for saline microscopic examination and Gram stain for Nugent score determination [15]. Three out of four Amsel criteria were necessary for the clinical diagnosis of BV and included (1) homogeneous, thin, grayish-white vaginal discharge (2), vaginal pH >4.5 measured by the clinician using pH test strips (EM Science, Gibbstown, NJ) (3), positive whiff-amine test, and (4) clue cells present on a wet mount of vaginal fluid [2]. A Nugent score of 0-3 on vaginal Gram stain was considered indicative of Lactobacilluspredominant normal vaginal flora, a score of 4-6 of abnormal, intermediate, vaginal flora, and a score of 7-10 of complex vaginal flora consistent with BV [15]. A third vaginal swab was collected and stored, with participant consent, at -80°C in individual tubes containing 3 mL of phosphate-buffered saline (PBS, pH 7.4) for future research related to the vaginal microbiome. Testing for sexually transmitted infections (STIs) and other vaginal infections was performed for women enrolled in each these studies, as previously described [45] (data unpublished for the M. genitalium cohort study).
Symptomatic women with BV were treated with metronidazole 500 mg orally twice daily for seven days, as recommended in the 2006 Centers for Disease Control STD Treatment Guidelines [46]. For the purposes of this study, stored vaginal swab specimens collected from women with three or four Amsel criteria and a Nugent score of 7-10 were used for DNA extraction and pyrosequencing.

DNA Extraction and Metagenomic Sequencing of Vaginal Microbiota
DNA extraction was performed on eligible vaginal swab specimens using the QIAamp ® DNA mini kit and following the manufacturer's instructions (QIAGEN, Valencia, CA). Bacterial tag-encoded FLX amplicon pyrosequencing (454 pyrosequencing) was performed by the Research and Testing Laboratories (Lubbock, TX). Briefly, 16S universal Eubacterial primers 28F: 5' -GAGTTTGATCNTGGCTCAG -3' and 519R: 5' -GTNTTACNGCGGCKGCTG -3' were used to amplify an approximately 500 base pair genomic locus which included the V1-V3 region of the 16S rRNA gene. Secondary PCR using tagged primers was done before massively parallel pyrosequencing on a Genome Sequencer FLX (Roche, Nutley, NJ).

Sequence Data Analysis and Composition
The QIIME data analysis package was used for 16S rRNA data analysis [47]. Quality filtering on raw sequences (994,670 reads) was performed and sequences which did not fulfill the following criteria were discarded: read length between 200-1000 bp, no ambiguous bases, mean quality score ≥25, and sequences matching the barcode. After quality filtering, a total of 956,793 sequences were used in the final analysis. The number of reads per sample ranged from 2,499 to 25,274 with an average of 8,467 reads per sample. The average read length was 492 bases across all samples; the read length mode was 509 bases. To remove sequence bias specific to pyrosequencing sequences, sequences were "denoised" using the flowgram clustering algorithm present in the QIIME package.
Sequences were grouped into operational taxonomic units (OTUs) using the clustering program UCLUST [48] at a similarity threshold of 0.97. In the dataset, 1,480 OTUs were identified in a total of 112 samples. The Ribosomal Database Program (RDP) classifier [49] was used to assign taxonomic category to all OTUs at confidence threshold of 80% (0.8). The RDP classifier uses the 16S rRNA RDP database which has taxonomic categories predicted to the genus level [50]. An OTU table (Table S1) was then generated combining all 1,480 OTUs, their taxonomic identification, and their abundance information. Since the RDP classifier provides taxonomic identification of OTUs only to the genus level, an attempt was made to predict species information for several highly abundant OTUs using BLAST (Basic Local Alignment Search Tool) searches against a non-redundant database at the National Center for Biotechnology Information website (available at http://blast.ncbi.nlm.nih.gov). If the lowest common taxonomy for a top blast hit was at the species level, then the species name is mentioned throughout the manuscript.
Two tables (Tables S2 and S3) were generated from the OTU table (Table S1) for further analysis. For the first table (Table S2), the abundance of OTUs present across all samples was normalized to account for different read depth across samples. OTUs whose abundance was not ≥0.1% in at least one of the 112 samples were subsequently filtered. This allowed the removal of "rare taxa" (in the context of sampling depth) [42] leading to a total of 342 most abundant OTUs.  Table S3 using the "heatmap.2" R package, version 2.10.1 (available at http://CRAN.R-project.org/ package=gplots). Sexual risk behavior groups, pH (categorized at an interval of 0.5), and Nugent scores were included in the heatmap. Groups of vaginal bacterial communities were identified based on clustering patterns and abundance of taxa. Alpha diversity within the samples was measured using Shannon's Index implemented in QIIME and beta diversity among the samples was measured using Unifrac analysis [51].

Statistical Analysis
Statistical analysis was performed using SAS software v9.2 (SAS Institute, Cary, NC). Table S2 was used to summate similar taxonomic groups. Summary statistics (sum, percentages, mean, standard deviation (SD), minimum, and maximum) were calculated for the overall cohort and for each sexual risk behavior group. Taxonomic groups in Table S2 with an abundance of ≥2% across the entire sequence library were considered for further analysis. Generalized linear models with negative binomial distributions and log link functions were used to estimate means and relative ratios (RR) of bacterial concentrations for comparing sequence reads between the sexual risk behavior groups. Huber-White robust standard errors and 95% confidence intervals (CIs) were calculated. Alpha was adjusted for multiple-comparison contrasts between sexual risk behavior groups using Bonferroni corrections. A Forest plot [52] was constructed for visual comparison of the relative ratios of bacterial concentrations between sexual risk behavior groups for taxonomic groups with an abundance of ≥ 2% across the entire sequence library.

Description of the Sample
A total of 112 women with BV (based on Amsel criteria and Nugent score) had stored vaginal specimens that were included in this study: 33 WSW, 35 WSWM, and 44 WSM ( Table 1). Mean age ± standard deviation (SD) of women in each behavioral group was: 25.1 ± 5.0 (WSW), 23.4 ± 4.4 (WSWM), and 23.3 ± 3.8 (WSM). All women in the WSW and WSWM groups were African American (AA); 93.2% of women in the WSM group were AA and 6.8% were Caucasian. The mean vaginal pH ± SD from each risk behavior group was 5.8 ± 0.4 (WSW), 5.6 ± 0.5 (WSWM), and 5.6 ± 0.4 (WSM) and was not statistically significant among sexual risk behavior groups. In addition, there were no significant differences among the Amsel or Nugent scores between groups.
Cluster analysis revealed four major groups of bacterial communities (I -IV) (Figure 2) with the dominant taxonomic group in each being Lactobacillus (I), Lachnospiraceae (II), Prevotella (III), and Sneathia (IV). Three samples were considered as outliers with each having a major proportion of unidentified bacteria, Enterobacteriaceae, or Serratia. The mean Nugent score ± SD of each cluster was 8.4 ± 1.1 (I), 9.6 ± 0.8 (II), 9.1 ± 1.1 (III), and 8.7 ± 1.1 (IV) (data not shown).  Figures 2 and 3 depict the cluster abundance among the sexual risk behavior groups. Fisher's exact test was used to determine if there were any significant differences in the overall cluster distribution between the sexual risk behavior groups however no statistically significant differences existed (p=0.09). Statistical comparisons between cluster distribution and individual sexual risk behavior groups (i.e. WSW and WSM) were not possible due to the small sample size in each group. Figure 1 depicts the most abundant taxonomic groups in the bacterial communities, stratified by sexual risk behavior group. Overall, Lachnospiraceae, which primarily includes BVAB1, was the most abundant taxonomic group detected in all sexual risk behavior groups.

Microbiota Distribution in the Sexual Risk Behavior Groups
Differences in these taxonomic groups between sexual risk behavior groups are represented on a Forest plot (Figure 4). With regards to comparisons between WSW/WSM groups, the abundance of Atopobium (RR=0.24; 95% CI 0.11-0.54) and Parvimonas (RR=0.33; 95% CI 0.11-0.93) were significantly lower in the WSW group compared to the WSM group. When comparing WSW/WSWM groups, the abundance of Prevotella (RR=1.77; 95% CI 1.10-2.86) was significantly higher in the WSW group than the WSWM group. Finally, when comparing WSWM/WSM groups, the abundance of Atopobium (RR=0.41; 95% CI 0.18-0.88) was significantly lower in the WSWM group compared to the WSM group. There were no significant differences with regards to the abundance of Lachnospiraceae, Sneathia, Lactobacillus, Megasphaera, or Dialister between sexual risk behavior groups.
The alpha diversity measurement, depicting diversity within the samples, demonstrated that the WSM group had the highest diversity of bacterial taxa compared to the WSW (p=0.04) or WSWM (p=0.01) groups ( Figure 5). In addition, the  unweighted beta diversity Unifrac analysis (Figure 6), depicting diversity between the samples (principle components PC1, PC2, and PC3, accounting for the highest percentage of diversity seen, captured 27% variation), showed a distinct cluster within the WSM group. A gtest (which determines whether OTU presence or absence is associated with a category) was performed between the WSM group and the other groups (WSW, WSWM) which revealed that several OTUs were present significantly more often in the WSM group (Bonferonni corrected, p<0.05) (data not shown). The WSM group was more likely to have the taxonomic groups Atopobium (OTU 519, OTU 313, OTU 19, OTU 492) and Peptoniphilus (OTU 1120) and less likely to have the taxonomic group Plantibacter (OTU 319) than the WSW and WSWM groups.

Discussion
This study investigated the vaginal microbiome of sexual risk behavior groups of women with BV. Lachnospiraceae (primarily BVAB1) was the most abundant taxa across all sexual risk behavior groups (Figure 3). BVAB1, a member of the Clostridiales order of anaerobic Gram-positive rods, is a fastidious bacterial taxon which has not yet been cultivated in vitro [38]. Its presence has been found to be highly specific for BV [38] although not all studies have found such a high abundance of BVAB1 among women with BV as our study [43]. It is possible that the abundance of BVAB1 varies based on the characteristics of the population of women with BV studied (i.e. race, ethnicity, sexual practices, etc.); further research is necessary to better understand these differences. In addition, it was interesting to note that, although BV is characterized by a lack of H 2 O 2 -producing lactobacilli, Lactobacillus (predominately L. iners) was still detectable among women in each of the sexual risk behavior groups in this study. L iners, found to be present in both normal and BV microflora [42,53,54] as well as the dominant vaginal bacterium after BV treatment [12,55,56], is hypothesized to facilitate transitions between BV and non-BV states [37,57].
This study found several significant differences in the vaginal flora of predominately African American WSW, WSWM, and WSM with BV. Atopobium (primarily A. vaginae) was significantly more abundant among WSM than WSWM or WSW (Figures 3, 4). This finding was corroborated by the principal coordinate analysis (Figure 6) which showed a distinct cluster within the WSM group with Atopobium being one of the major OTUs in this cluster. It could be hypothesized that A. vaginae is not critical to the BV syndrome as WSW (who frequently have a higher prevalence of BV than WSM) have a much lower abundance of A. vaginae than WSM, as shown in this study. It may be that A. vaginae's primary niche is in males, that it is sexually transmitted, and that it survives best in the presence of other BV microbiota. In addition, Parvimonas was also significantly more abundant among WSM than WSW. Parviomonas is a well-known oral pathogen known to be associated with odontogenic infections [58]. Finally, Prevotella spp. (primarily P. bivia str, P. amnii, P. timonensis, and P. buccalis), were significantly more common among WSW than WSWM. Prevotella spp. have been isolated from the oral cavity, the upper respiratory tract, and the urogenital tract in humans [59,60]. P. bivia has been noted in women with BV, endometritis, pelvic inflammatory disease, and peri-rectal abscess [61,62] and has been found to produce high lipopolysaccharide concentrations in vaginal fluid which may be The number of women from each sexual risk behavior group is indicated in parentheses. Note that 3 samples from the WSWM group contained taxonomic groups that were considered as outliers (i.e., had a major proportion of unidentified bacteria, Enterobacteriaceae, or Serratia) and were not included in this analysis.
an important factor in the pathogenesis of BV and subsequent complications of infection [62].
Taken together, these data seem to suggest that certain Gram-positive anaerobes such as Atopobium and Parvimonas were significantly more common among WSM with BV than among women reporting sex with women. Additionally, WSM had the highest diversity of bacterial taxa among the sexual risk behavior groups and formed a distinct principal coordinate analysis cluster significantly distant from the other two groups. Potential explanations for these results are not immediately obvious although it could be hypothesized that differential participation in certain types of at-risk sexual practices could be contributing to differences noted in organism composition and BV prevalence among sexual risk behavior groups. Assuming that BV is sexually transmitted, it is currently not well known whether one type of sexual practice (i.e. receptive oralvulvovaginal sex, receptive oral-anal sex, digital-vaginal sex, digital-anal sex, use of sex toys, etc.) is more important in the pathogenesis of BV than another (i.e. penile-vaginal sex, penile-anal sex, etc.). It is interesting to note that, among WSW participating in a prospective cohort study investigating the risks for BV acquisition, a direct dose-response relationship between increasing numbers of episodes of receptive oralvulvovaginal sex and BV has been noted [63]. In another cohort of WSW, vaginal insertive use and sharing of sex toys were associated with decreased quantities of vaginal lactobacilli and higher risk of colonization with G. vaginalis [64]. Alternatively, studies of primarily heterosexual women have shown an increasing risk of BV in women who have more frequent penile-vaginal sex [19,65]. We did not have detailed   sexual practices data (oral, vaginal, and/or anal) with female and/or male sexual partners available for all women in this study and thus cannot make any correlations between BVassociated bacteria and sexual practices at this time. Future studies of the vaginal microbiome among sexual risk behavior groups of women with BV should include this pertinent information. Additionally, it is not well known whether certain BV-associated bacteria prefer the vaginal fluid environment to the male genitourinary environment or whether transmission of infected vaginal fluid is a more efficient mechanism for the acquisition of BV among WSW than transmission from the male genitourinary tract among WSM. Each of these factors could potentially influence differences in BV prevalence noted among sexual risk behavior groups of women.
This study has several limitations. First, due to feasibility issues and lack of follow-up in the parent studies from which the WSW, WSWM, and WSM specimens were collected, we did not have paired samples to assess the microbial diversity in the same women pre-and post-resolution of BV. Second, because this was a pilot study with budgetary constraints, we were not able to include an age-matched control group of WSW, WSWM, and WSM with normal vaginal flora to reveal the baseline microbiome in each of these groups of women. Third, the influence of the menstrual cycle and hormoneinduced changes in the microbial flora were not able to be controlled for in this study as specimens were obtained from women participating in two different parent studies; future studies in these groups of women should control for phase of the menstrual cycle. Fourth, DNA sequences for amplification of the V1-V3 region of the 16S rRNA gene were generated from a single set of primers which are considered to be conserved among most Eubacteria, however, some taxa such as Gardnerella may have divergent sequences which did not efficiently amplify with the primers used here. Previous studies using primers for conserved sequences in 16S rRNA genes have also underrepresented or not consistently amplified 16S rRNA sequences from Gardnerella [54,66]. This is an important limitation as G. vaginalis has been the bacterium most extensively studied in relationship to BV pathogenesis [14] and is present in 95-100% of clinically diagnosed BV cases using cultivation-dependent methods [34,67,68]. It has recently been detected in the oral and anal cavities of WSW who subsequently developed incident BV, leading to the hypothesis that it may be acquired vaginally from preexisting reservoirs at extravaginal sites [69]. Future studies among sexual risk behavior groups of women with BV should include primers that detect G. vaginalis in order to further explore its association with specific sexual practices. Other as yet unknown bacterial taxa associated with BV may have also been missed in this study by using restricted primer sets for DNA amplification. Of note, a recent study by Srinivasan et al which targeted the V3/V4 region of the 16S rRNA gene for broad-range PCR and sequenced the V4 region allowed for high resolution phylogenetic analyses of bacterial communities in women with BV [37]. In addition to not having data on specific sexual practices, we did not have data on recent antibiotic use, hormone use, or vaginal product use which could also have impacted the vaginal microbiomes (with regards to inflammation, pH, and microbial community content) of sexual risk behavior groups of women with BV in this study. Finally, the majority of women in our study were of African American race. Differences in the composition of vaginal microbial communities of women of different racial and ethnic groups with and without BV have been noted [37,42,54] and our results may not be generalizable to other racial and ethnic sexual risk behavior groups of women with BV. Overall, further studies with well-controlled prospective samples, age-matched non-BV controls, and robust questionnaire data are needed to understand whether the differences in vaginal flora identified in this study are valid.
Despite these limitations, the results of this study demonstrate that the microbiology of BV among women in different sexual risk behavior groups is heterogeneous. How this variation arises is unclear, but it is possible that vaginal, oral, and male urethral secretions are niches for unique communities of micro-organisms which may be transmitted sexually and play important roles in the pathogenesis of BV. Teasing out the specific contributions of sexual practices and sexual partners to the altered vaginal flora of BV should be a priority for future research. Additionally, because treatment of BV is directed at anaerobic organisms, it is critical to define the major anaerobic components of the vaginal flora in BV, as some BV-associated bacteria such as Atopobium spp. have been shown to be resistant to commonly used agents such as metronidazole [70]. It is possible that different sexual behaviors will define groups of women with BV that respond differently to antibiotic therapy. Table S1.