Five adult female non-pregnant donors were enrolled in this study, after they signed an informed consent, approved by Sidra Institutional Review Board committee (IRB protocol #1709014294). The donors were self-reported as healthy. All methods were performed in accordance with the relevant guidelines and regulations. Each donor was asked to collect 4 vaginal samples using Copan ESwab (Copan Diagnostic Inc., Murrieta, CA, USA). The Copan ESwab is a collection and transportation system which incorporates a modified Liquid Amies transport medium and a flocked swab. The samples were self-collected by the volunteers as follows: using the non-dominant hand to open the labia, the swab is inserted into the vagina, then twisted several times while inside of the vagina. Once the swab was collected, it was placed immediately into the ESwab transport tube containing transport medium. Swabs were then transported directly to the laboratory (within 2 h) and stored at −80C before processing.

To extract microbial DNA from all vaginal samples collected, four different DNA extraction protocols were tested and evaluated:

Qiagen DNeasy Blood and Tissue (Qiagen, Venlo Netherlands) referred to as method #1 hereafter.

The first three protocols (methods #1–3) use a pre-centrifugation step (10 min at 7,500 rpm) in order to collect the pellet prior to cell lysis. All protocols are described in Figure 1.

The assessment of DNA quantity and quality was carried out using Nanodrop. The absorbance at 260 and 280 nm wavelengths was measured. Most of the isolated DNA samples had an OD260/OD280 ratio between 1.72 and 2.35. We then checked the DNA integrity using LabChip GX (PerkinElmer, Waltham, Massachusetts, United States). The Gel-Dye solution, DNA samples and DNA ladder were prepared according to the manufacturer's instructions; the run data was compared to the electropherogram of a typical high molecular weight ladder and assessed for quality. A genomic DNA (gDNA) quality score (GQS) was calculated for each sample. The GQS is derived from the size distribution of the gDNA and it represents the degree of degradation of a given sample, with a score of 5 corresponding to intact gDNA and a score of 0 corresponding to a highly degraded gDNA.

To evaluate the relative abundance of bacteria as compared to human DNA in samples tested, a quantitative real time PCR (qPCR) was performed to assess the presence and the relative quantity of microbial DNA. Universal bacterial primers amplifying the V3 region of 16S rDNA gene (341 F and 534 R) were used to amplify approximately 194 bp as described previously (Whiteley et al., 2012). Primers for the Firmicutes family were used as representative for the vaginal microbial DNA (Yang et al., 2015). STAT2 (RefSeq: NC_000012.12) primers were used as representative for the human DNA. The primers used are detailed in Table 1.

Each PCR reaction was performed in a 20 ul final volume containing 200 nM of forward and reverse primers, 10 μL of 2X KAPA Master Mix (Sigma-Aldrich, St. Louis, WI, USA). Input DNA was 50–100 ng/reaction. Cycling conditions were as following: initial denaturation at 95°C for 10 min; 30 cycles at 95°C for 30 s, 63°C for 30 s, and 72°C for 30 s; final step at 72°C for 10 min and hold at 4°C. Reactions were run on the LightCycler 480 System (Roche, Basel, Switzerland). To verify primers specificities, melting curves were generated at the end of each PCR reaction. Fluorescent data was acquired during the extension phase. After 30 cycles, a melting curve for each gene was generated by increasing the temperature from 60 to 95°C (1°C for each step), while the fluorescence was measured. Samples were run in triplicates. For each experiment a no-template reaction was included as a negative control. Two human DNA samples were used as controls. The qPCR data was analyzed using the 2^∧-(delta Ct) method (comparing Ct of genes of interest against STAT2).

The V1-V3 regions of the 16S rDNA were amplified using various forward primers: 27F with 12 bp golay barcodes containing a specific Illumina 5′ adapter for each sample and a common reverse primer 534 R (Zheng et al., 2015). In brief, PCR was performed in triplicate in a 50 μL reaction mixture containing 10 ng of template DNA and 2x KAPA HiFi HotStart Ready Mix. The following thermal cycling conditions will be used: 5 min of initial denaturation at 94°C; 30 cycles of denaturation at 94°C for 30 s, annealing at 62°C for 30 s, and elongation at 72°C for 30 s; and a last step at 72°C for 10 min. The amplified PCR products of ~650 bp in size from each sample were pooled in equimolar concentrations. This pooled PCR product was purified using the FlashGel™ DNA System (LONZA). High throughput sequencing was performed on an Illumina MiSeq 2 × 300 platform (Illumina, Inc. San Diego) in accordance with manufacturer's instructions. Image analysis and base calling were carried out directly on the MiSeq.

Demultiplexed sequencing data were analyzed using QIIME software v1.9.0 pipeline (Caporaso et al., 2010). FASTQ files were converted into FASTA files, and all demultiplexed files were concatenated into a single file. Further analysis was performed as previously reported (Murugesan et al., 2015; Garcia-Mena et al., 2016). Sequence alignments were done against the Greengenes core set (Desantis et al., 2006). Alpha diversity was measured by R software, using the phyloseq package (Mcmurdie and Holmes, 2013). Beta diversity was represented using weight UniFrac distance measure (Lozupone et al., 2007) and contributions to the differences in the beta diversity were presented as principle coordinate analysis as prposed in QIIME.

Data were analyzed using Medcalc 11 (Medcalc Software, Stata Software). Paired Student's t-tests were used to compare means. Alpha diversity measures such as Observed, Chao1, Shannon and Simpson indices were calculated using minitab17 (Minitab statistical software). Kruskal wallis tests were used to compare the statistical difference among diversities between the described DNA extraction protocols. P < 0.05 were considered statistically significant. We applied the Linear Discriminative Analysis (LDA) Effect Size (LEfSe) tool which uses Kruskal-Wallis and estimates the effect size of the comparisons to identify the significant shift in profile to distinguish the difference in extraction methods (Segata et al., 2011). Adonis and Distance based redundancy (db-RDA) method used to calculate distance matrix difference between the extraction methods included in this study using beta diversity parameters (Sobel, 1997).

After isolation of the DNA from the vaginal swabs, the quality and quantity of the DNA extracted were evaluated using Nanodrop. Most of the isolated DNA samples had a similar OD260/OD280 absorbance ratio of approximately 2 (Supplementary Figure 1). We further confirmed the quality of the DNA isolated by running all the samples on the LabChip and by assessing their GQS (a representative image of GQS scoring is reported in Supplementary Figure 2). As described in Figure 2A, significant differences in the DNA yield were observed when we compared the four extraction protocols; method #1 resulted in the highest DNA yield with an average of 5.96 ug, followed by methods # 4, 2, 3, respectively.

The size of the DNA extracted was ~40 kb. The GQS was calculated for each sample and plotted in Figure 2B. Method #1 resulted in the highest GQS, with an average of 4.24, followed by methods #2, 3, 4, respectively (Figure 2B).

The abundance of bacterial DNA was assessed by qPCR using universal primers for the bacterial 16S rDNA gene (341 F and 534 R) and primers for the Firmicutes family which represents one of the most predominant bacteria identified in the vaginal environment in healthy women of reproductive age (Ling et al., 2010). Data was calculated with the 2^∧-delta Ct method and plotted in a logarithmic scale (Figure 3).

Substantial variability in microbial DNA abundance was found across all the protocols tested. Overall, bacterial DNA was detected in all the extraction protocols tested, with some variabilities (Figure 3A). Method #2 gave the highest yield of bacterial abundance followed by methods #1, 3 and 4, respectively (Figure 3A). A positive correlation of Firmicutes and 16S rDNA gene abundances across all the 4 protocols was observed (Figures 3B–E), suggesting that indeed bacteria of Firmicutes family represent a good portion of the vagina flora.

All the DNA samples isolated from the four DNA extraction methods were processed for sequencing using Illumina Miseq platform after amplifying the V1-V3 regions of the 16S rDNA gene. Clustering and annotation of the operational taxonomic units (OTUs) was performed using the same pipeline for all the samples. Average number of sequences count per sample was also calculated (Supplementary Figure 3) and a significant difference was observed when comparing extraction method #1 and #2. No statistical significance was observed when we compared other extraction methods used. A comparison of OTUs from phylum to genus levels was performed. No statistical significance was observed in the various extraction methods used when we compared the total OTUs count (Supplementary Figure 3). However, our analysis revealed that regardless of the extraction protocol used, Firmicutes remains as expected the most abundant phylum with an average of 98% of all OTUs (Figure 4C). At the genus level, Lactobacillus species was the most abundant genera (Figure 4D). Prevotella spp., Streptococcus spp., Ruminococcaceae and Coriobacteraceae were the least abundant genera across the four DNA extraction methods (Figures 4A–D). Relative abundance of the vaginal microbiota profiles of each donor using different extraction methods showed that there is a shift observed in method #2, method #3, and method #4 as compared to method #1 (Figures 4A–D). Overall, Powersoil DNA extraction methods exposed more genera especially with method #3 and method #4 as compared to method #1 (Figure 4B).

LEfSe analysis revealed the significant differences in the abundances of Coriobacteriaceae, Atopobium spp., and Dialister spp., in method #2 and #4 when compared with method #1. Likewise, Method #3 showed significant difference in the abundance of Coriobacteriaceae and Atopobium spp., in comparison with method #1 (Figure 4E).

Diversity indices were calculated to describe the complexity of samples (alpha) and to differentiate between the samples assessed (beta). Alpha diversity was measured using four indices (Figure 5). Alpha diversity measures (Figure 5) indicated that significant increase was observed in both chao1 (P = 0.021, P < 0.05) and observed species richness (P = 0.021, P < 0.05) in the vaginal microbial community structure when method #4 was used, and significantly higher chao1 (P = 0.034, P < 0.05) and observed species richness (P = 0.034, P < 0.05) in method #3 as compared to method #1. This result is also concordant with the Shannon diversity index in method #3 (P = 0.021, P < 0.05) and in method #4 (P = 0.034, P < 0.05) (Figure 5).

Principle coordinate analysis (PCoA) of beta diversity was performed to examine whether each method clustered in a distinct pattern using both the taxonomic Bray-Curtis dissimilarities (Figure 6A) and phylogenetic weighted UniFrac distances (Figure 6B). The Adonis method was applied to calculate the variation of distances among different extraction methods. The Bray-Curtis method showed that method #1 was significantly different from other methods (Figure 6A). No significant clustering was observed when we compared the four DNA extraction methods using weighted UniFrac distances (Figure 6B). The PCoA of Bray-Curtis dissimilarities (Figure 6C) and phylogenetic weighted UniFrac distances (Figure 6D) of each donor also showed that there was no donor specific clustering irrespective of methods of extraction. A distance based redundancy analytic method was applied to calculate the variation of distances among different extraction methods. The Bray-Curtis method showed that method #1 was significantly different from other methods (Figure 7A). No significant clustering was observed when we compared the four DNA extraction methods using weighted UniFrac distances (Figure 7B).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.