Characterization of vaginal microbiota in Thai women

Subject selection and sample collection

Twenty-five women attending the gynecologic out-patient clinics for cervical cancer screening, Srinagarind Hospital, Faculty of Medicine, Khon Kaen University, were enrolled in this study. The inclusion criteria used for subject selection specified non-pregnant women aged 20–45 years with regular menstruation, with normal vaginal examination and without any clinical symptoms on examination by one gynecologist. The exclusion criteria were women with the following conditions: having menstruation or sexual intercourse within the previous 24 h, use of antibiotics or vaginal antimicrobials within the last month, presence of any intravaginal product, use of a vaginal douche in the past week, and active infection or diagnosis of any defined sexually transmitted disease. The subjects selected by the above criteria were defined as normal women. All participants provided written informed consent and filled out a questionnaire detailing their level of education, employment, sociodemographic characteristics, dietary habits, reproductive health habits, and sexual behaviors (including number of children), income per month, alcohol consumption per week, exercise per week, sexual activity per week, number of sexual partners in the past, regularity of their menstrual cycles in the last six months, current hormone contraceptive use, douching practice in the past month, abnormal vaginal discharge in the past month and female sterilization. This study was approved by the Ethics Committee of Khon Kaen University (HE601017).

A total of four points at mid right and left, anterior and posterior vaginal wall were swabbed with sterile cotton-tipped applicators by one gynecologist. The swabs were placed in sterile containers with phosphate-buffer saline (PBS) and stored at −80 °C until analysis.

Genomic DNA extraction from vaginal swabs

Genomic DNA was extracted using the method described by Ravel et al. (2011) with slight modification. In briefly, the vaginal swabs were thawed on ice and adhered material resuspended by vigorous vortexing for 5 min. Bacterial cells were lysed by incubation for 1 h at 37 °C in TE50 buffer (10 mM Tris HCl, 50 mM EDTA, pH 8.0) containing 10 mg/mL of lysozyme (Merck, Darmstadt, Germany), 25,000 U/mL of mutanolysin (Sigma Aldrich, St. Louis, MO, USA) and 4,000 U/mL of lysostaphin (Sigma-Aldrich), followed by beating with 40–400 µm glass beads (NucleoSpin^® Bead Tubes Type B; Macherey-Nagel, Düren, Germany) at speed 10 for 2 min with air cooling using the Bullet Blender^® Blue Homogenizer (BBX 24B; Next Advance, Inc., Troy, NY, USA). After bead-beating, the crude lysates were processed for purification of genomic DNA using QIAamp^® DNA Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s recommendations. Integrity, size and concentration of the purified genomic DNA was determined using Fragment Analyzer™ (Advanced Analytical Technologies, Inc., Ankeny, IA, USA).

16S rRNA gene sequencing by NGS

16S rRNA gene sequencing was done as previously described (Barb et al., 2016) with minor modifications. The seven 16S hypervariable regions were amplified with two primer sets, one set for amplifying the V2, V4, V8 hypervariable regions and the other set for amplifying the V3, V6-7, V9 hypervariable regions. Primers were supplied with the Ion 16S™ Metagenomics Kit (Life Technologies, Grand Island, NY, USA) and used according to the manufacturer’s recommendations. The combined PCR products were processed to make the DNA library by using the Ion Plus™ Fragment Library Kit (Life Technologies, Carlsbad, CA, USA) with sample indexing using the Ion Xpress™ Barcodes Adapters 1-32 Kit (Life Technologies). The adapter-ligated and nick-repaired DNA was amplified with the following steps: 1 hold at 95 ° C for 5 min; 5–7 cycles of 95 °C for 15 s, 58 °C for 15 s, 70 °C for 1 min; hold at 4 °C. The PCR products were then purified using Agencourt^® AMPure^® XP Reagent (Beckman Coulter, Inc, Atlanta, GA, USA). The processed libraries were quantified using a Bioanalyzer^® instrument (Agilent Technologies, Santa Clara, CA, USA) with the use of the Agilent^® High-Sensitivity DNA Kit. Equal volumes of all 25 samples were pooled and processed with the Ion PGM™ Hi-Q™ OT2 Kit in the Ion OneTouch™ 2 System (Life Technologies) according to the manufacturer’ instructions. Sequencing was performed in an Ion 318™ Chip (Life Technologies) using the Ion PGM™ Sequencing 400 Kit on the Ion PGM™ System (Life Technologies).

Sequence data analysis

Base calling and run demultiplexing were performed using Torrent Suite™ Software version 5.0.5 (Life Technologies) with default parameters. All sequences were trimmed according to the quality. To pass, a sequence read had a following criteria: (1) a perfect match to a barcode sequence and the primer; (2) was at least 200 bp in length; (3) had no more than two undetermined bases; and (4) had an average quality score greater than Q20. Sequences were analyzed using Ion Reporter™ Software version 5.6 with the 16S Metagenomics workflow module (Life Technologies). Clustering of operational taxonomic units (OTUs) and taxonomic classification were performed based on two-step Basic Local Alignment Search Tool (BLAST) that maps to two separate reference libraries. In the first step, reads were aligned to the MicroSEQ^® 16S Reference Library v2013.1 database (Woo et al., 2003). Subsequently, any unaligned reads subject to second alignment to the Greengenes v13.5 database (DeSantis et al., 2006). At least 10 of unique reads were valid and ≥90% for the alignment coverage needed between hit and query. Genus- and species-level identifications were accepted at ≥97% and ≥99% sequence identity, respectively (Drancourt et al., 2000). In case of any unusual microorganism, the sequence reads were reanalyzed by BLAST with the National Center for Biotechnology Information (NCBI) taxonomy database. Any OTUs represented by only one sequence (singletons) and those with fewer than 10 reads in the sample were excluded from further analysis. Alpha- (within sample) and beta- (among samples) diversities were calculated using Quantitative Insights Into Microbial Ecology (QIIME) software (Caporaso et al., 2010). Alpha diversity metrics, consisting of the Shannon diversity index, observed species and Chao1 index, were graphed using GraphPad Prism version 5.01 (GraphPad Software Inc., San Diego, CA, USA). The beta diversity metric was calculated according to Bray-Curtis dissimilarity index (Bray & Curtis, 1957) and displayed through a principal coordinates analysis (PCoA) plot. Heat map and hierarchical clustering tree were generated using R Studio (RStudio Team, 2018) employing R version 3.4.1 (R Core Team, 2017). Raw sequences were deposited in the NCBI Sequence Read Archive (SRA) accession number SRP158176.

Statistical analysis

Statistical analysis was performed using IBM SPSS Statistics for Windows, version 19.0 software (SPSS Inc., Chicago, IL, USA). Demographic characteristics between two groups were compared using the chi-squared or Fisher’s exact test. Mann–Whitney U test for non-parametric data was used for bacterial abundance and diversity comparisons between two groups. P-values less than 0.05 were considered as statistically significant.

Participant characteristics

The median age and body mass index of the enrolled subjects were 39 years (ranging from 31 to 45) and 21.2 kg/m² (ranging from 16.0 to 35.8), respectively. Subjects fell into two groups: lactobacilli-dominated (LD) and non-lactobacilli dominated (NLD) groups (see below). No significant association between sociodemographic characteristics and VMB group (LD or NLD) was found, as shown in Table 1.

Analysis of VMB in the studied subjects

High-throughput sequencing data of all vaginal samples derived from Ion Torrent PGM yielded a total of 4,775,332 raw sequences with an average of 191,013 reads per sample (ranging from 33,969 to 382,574 reads). We discarded any reads shorter than 150 bp and those OTUs represented fewer than 10 times in the sample. This left a total of 4,189,265 valid reads with an average of 167,570 reads per sample (ranging from 28,825 to 333,480 reads) for analysis. The average length of included reads was 243 bp (ranging from 237 to 253 bp). A total of 3,553,612 mapped reads were clustered into OTUs and given taxonomic identification. Forty-one genera and 72 species were identified (Tables S1 and S2). VMB were clustered into two distinct groups: LD (n = 14) and NLD (n = 11). In the LD group, the average of relative abundance of Lactobacillus species is 90% and the relative abundance of this species among the LD individuals is ≥60% and could be further divided into two subgroups: one dominated by L. iners (n = 10) and the second subgroup (n = 4) was dominated by L. crispatus, L. gasseri, L. jensenii and L. johnsonii. The NLD group could similarly be divided into two subgroups: one dominated by Gardnerella vaginalis (n = 3) and the other (n = 8) dominated by anaerobic bacteria of a mixture of several species including G . vaginalis, Atopobium vaginae and Pseudomonas stutzeri . Figure 1 represents the hierarchical clustering tree (Fig. 1A) and the heat map (Fig. 1B) of bacterial taxa identified from 25 vaginal samples. Table S3 indicates the mean percentage representation of each bacterial genus present in these two groups.

In this study, alpha-diversity for determination of bacterial diversity within individual women used three metrics: Shannon diversity index (estimated evenness and richness), observed species (observed bacterial richness) and Chao1 (estimated bacterial richness). The Shannon diversity index differed significantly between the two groups (P < 0.001) (Fig. 2). The observed and estimated numbers of OTUs in individuals in the NLD group were significantly higher than those in the LD group (both P < 0.001) (Fig. 3). All these results indicated that the bacterial communities in NLD individuals have more species diversity than those in LD individuals. In beta-diversity analysis, the PCoA plot based on taxonomic profiles showed that some distinct clustering appeared to overlap between two groups (Fig. 4). This result suggested that the VMB of both LD and NLD individuals did not tightly cluster by group. In addition, the relative distances based on the Bray-Curtis dissimilarity index showed significantly greater distances among LD individual samples than among NLD individual samples (0.82 ± 0.10 versus 0.72 ± 0.09, P = 0.021). Since the distance measurement indicates the degree of similarity among samples (Van de Wijgert & Jespers, 2016), this result indicates that VMB in LD individuals resemble one another more closely than do those in NLD individuals.

The number of species unique to either LD or NLD, or common to both groups, were 14 (19.44%), 41 (56.95%) and 17 (23.61%), respectively (Fig. 5). The 17 common bacterial species found in both groups were A . vaginae, Finegoldia magna, G . vaginalis, L. crispatus, L . gasseri, L . iners, L . johnsonii, L . kefiranofaciens, L . vaginalis, Megasphaera spp., Peptoniphilus harei, Peptostreptococcus anaerobius, Prevotella bivia, P . timonensis, Pseudomonas stutzeri, Streptococcus intermedius and Ureaplasma parvum. Fig. 6 represents the most abundant bacterial species detected in the vaginal samples of LD and NLD groups. G.vaginalis, A . vaginae and P. bivia were significantly associated with the NLD group (P < 0.05) while L. iners was significantly associated with the LD group (P < 0.05). The bacterial species with abundance less than 1% are F . magna, L. kefiranofaciens, L. vaginalis, Peptoniphilus harei, Peptostreptococcus anaerobius, S . intermedius and U. parvum.