Participant enrollment
We report the results from 31 participants enrolled at the Gynecologic Division, Mayo Clinic, Rochester, MN under an IRB approval protocol (12–004445). The inclusion criteria were the following: 18 years of age or older; women undergoing hysterectomy by any standard surgical approach; undergoing hysterectomy for benign disease, hyperplasia, or any stage of endometrial cancer. Patients with any of the following criteria were excluded from our study: women who were pregnant or nursing; had taken antibiotics within two weeks preceding surgery; surgeon using morcellation during the hysterectomy procedure, due to the size of the uterus or for any other reason. Upon enrollment the participants were requested to fill out an optional questionnaire about sexual and reproductive health and history. The metadata from the questionnaires was stored at REDCap [14]. Cancer participants were also requested to provide a stool sample for the search for putative endometrial cancer signatures.
Sample collection
Vaginal and cervical samples
All participants were requested not to douche with betadine on the day of surgery or the day immediately preceding it. All the vaginal and cervical swabs and scrapes were collected by the surgeon (with guidance on site by the research team) immediately after the administration of anesthesia and immediately preceding the standard pre-surgical betadine douche. Both the vaginal and cervical swabs were performed with three sterile Dacron swabs each and placed in a sterile tube with 1 mL of Tris-EDTA (TE) buffer kept on dry ice until storage at –80 °C. One of the vaginal swabs was used for immediate on-site vaginal pH measurement with a Hydrion measuring pH tape. The scrapes were performed using sterilized (autoclaved at 121 °C for 20 min) pap smear spatulas and placed in sterile tubes with TE buffer kept in dry ice until storage at –80 °C.
Uterine, Fallopian, and ovarian samples
Once removed, the uterus, Fallopian tubes, and ovaries were handed by the surgeon to the instrumentalist nurse who placed them inside a sterile transport bag and into a closed sterile container. The research team then transported the container to the pathology lab (within the same clean area) where the organs were handed to a pathologists’ assistant (PA) to be processed under sterile conditions. The grossing station where the specimen was processed was sterilized by the research team, including all the tools needed by the PA for handling. The PA used surgical gloves and mask when handling the specimen. The PA performed a bilateral cut of the uterus and splayed it. The research team advanced to the collection of the uterine swabs (Dacron) and scrapes (sterilized pap smear spatulas) and documentation (by placement of push pins in sampled locations and digital photograph). The PA then proceeded to the aseptic collection of samples needed for the diagnosis and, once complete, the research team collected the uterine, Fallopian, and ovarian biopsies (approximately 4 mm of tissue was collected per biopsy by the use of a pair of sterile tweezers, scalpel, and surgical ruler). Each collected sample was placed in a sterile tube with 1 mL of TE buffer and kept on dry ice until storage at –80 °C. A petri dish with Lysogeny broth (LB) was kept open on the grossing station during sample collection to detect any possible airborne contamination of the specimen. The LB was swabbed and the swab was stored in a tube with 1 mL of TE and kept on dry ice until storage along with all the other samples.
Sample processing
Once thawed, the swab and scrape samples were vortexed to bring the collected material into solution. The biopsy samples were macerated by the use of sterile pestles. The swab and scrape samples were centrifuged for 10 min at 10,000 g to collect the bacterial cells and the supernatant was discarded. All genomic DNA extractions were performed by using the MoBio PowerSoil Kit (MoBio Laboratories, Inc., Carlsbad, CA, USA) as described by the manufacturer; however, instead of vortexing, an MP FastPrep (MP Biomedicals, Solon, OH, USA) was used instead, for 60 s at 6.0 m/s, to obtain a more effective and rapid lysis of the cells. After extraction the DNA content was measured using High Sensitivity Qubit (Life Technologies Corporation, Carlsbad, CA, USA). The V3-V5 region of the 16S rDNA was then amplified through a polymerase chain reaction (PCR) as follows: 25 μL of Kapa HiFi (Kapa Biosystems, Woburn, MA, USA), 1.5 μL (10 uM) forward primer, 1.5 μL (10 uM) reverse primer, 50 ng of DNA with the remaining volume being added by molecular grade water (up to a final volume of 50 μL per reaction). The forward primer was the universal primer 357 F (5’GTCCTACGGGAGGCAGCAG3’) with the added construct on the 5’ end of the 5’ Illumina Adapter (5’AATGATACGGCGACCACCGAGATCTACAC3’) + Forward Primer Pad (5’TATGGTAATT3’) to a total sequence: 5’AATGATACGGCGACCACCGAGATCTACACTATGGTAATTGTCCTACGGGAGGCAGCAG3’ and the universal bacterial reverse primer was 926R (5’CCGTCAATTCMTTTRAGT3’) with an added construct on the 5’ end of the reverse complement of 3’ Illumina adapter (5’CAAGCAGAAGACGGCATACGAGATGCCGCATTCGAT3’) + Barcode (12 base pairs) to a total sequence: 5’CAAGCAGAAGACGGCATACGAGATGCCGCATTCGATXXXXXXXXXXXXCCGTCAATTCMTTTRAGT3’. The barcode introduced in the reverse primer construct was unique to each sample, functioning as a genetic ID for sequencing. The PCR cycle was the following: 95 °C for 3 min, 98 °C for 20 s, 70 °C for 15 s, 72 °C for 15 s, cycle repeated 34 times, and 72 °C for 5 min. The products of the amplification were verified by a TapeStation D1K Tape (2200 TapeStation Instrument, Agilent Technologies, Santa Clara, CA, USA) to be free of contamination and to contain the expected amplification size, approximately 700 base pairs. If the amplification was unsuccessful, the parameters of the reaction or cycle were adjusted in repeated attempts. In some cases (mostly biopsy samples) the amplification was not successful even after repeated attempts. The reduced number of microorganisms present in the upper reproductive tract is likely to justify this result and attests for the success of the sterile collection of the samples. In samples that failed 16S rDNA amplification, NEBNext Microbiome DNA Enrichment Kit (New England Biolabs Inc., Ipswitch, MA, USA) was used to separate the microbiome from the human DNA to increase the odds of a successful amplification from samples naturally enriched with human DNA (mostly tissue samples). Controls of both the DNA extraction and Microbiome Enrichment processes were performed and are shown in Supplement 5. Upon verification the PCR products were purified using Agencourt AMPure (Beckman Coulter, Brea, CA, USA). After purification the concentrations were measured using Qubit High Sensitivity. The 16S rDNA sequencing was performed by the MGF (Medical Genome facility at Mayo Clinic, Rochester) using a high-throughput next-generation Illumina MiSeq (San Diego, CA, USA) sequencing platform.
Sequence analysis
Sequence reads were aligned with our own custom multiple alignment tool known as the Illinois-Mayo Taxon Operations for RNA Dataset Organization (IM-TORNADO) that merges paired end reads into a single multiple alignment and obtains taxa calls [15]. IM-TORNADO then clusters sequences into operational taxonomic units (OTUs) using AbundantOTU+ [16].
Sequencing outcome
A total of 16,366,472 sequence reads (17,657–828,181 reads per sample) were obtained (mean of 199,591 ± 190,153 reads) after quality control. Further processing for visualization was performed using QIIME [17] and METAGENassist [18].
Data analysis
α-diversity and β-diversity analysis
To compare the microbiota composition between cohorts, we summarized the data using both α-diversity and β-diversity. α-diversity reflects species richness and evenness within bacterial populations. Two α-diversity metrics, the observed OTU number and the Shannon index, were investigated. Rarefaction curves were used to compare the α-diversity measures. The observed OTU number reflects species richness, whereas the Shannon index measures both species richness and evenness. β-diversity reflects the shared diversity between bacterial communities in terms of ecological distance between samples; different distance metrics provide distinctive views of community structure. Two β-diversity measures (unweighted and weighted UniFrac distances) were calculated using the OTU table and a phylogenetic tree (“GUniFrac” function in the R package GUniFrac) [19]. The unweighted UniFrac reflects differences in community membership (i.e. the presence or absence of an OTU), whereas the weighted UniFrac captures this information and also differences in abundance. Rarefaction was performed on the OTU table before calculating the distances.
To assess the association with α-diversity, we fitted a linear mixed effects model (LME) to the α-diversity metrics with a random intercept for each subject (“lme” function in R package “nlme”), adjusting for covariates if necessary. Wald test was used to assess the significance. To assess the association with β-diversity measures, we used a variant of PERMANOVA procedure (“adonis” function in the R “vegan” package), which is a multivariate analysis of variance based on distance matrices and permutation [20]. To retain the within-subject correlation, we used a block-permutation scheme, where samples from the same participant were assigned a different subject ID. Significance was assessed by 1000 permutations and the covariate was adjusted if necessary. Ordination plots were generated using non-metric multidimensional scaling (NMDS) as implemented in R (“metaMDS” function in the R “vegan” package).
To test for the correlation between organs, we used a permutation test based on Bray-Curtis distance with the test statistic calculated as the distance between the organs from different participants minus the distance between the organs from the same participant. We next permuted each participant for the same organ type using the same block-permutation scheme as above. The p value was calculated as the percentage of permutations that produce a test statistic more extreme than what is observed. To identify the taxa shared by both organs, we used a taxon-specific Euclidean distance, defined based on the presence and absence of a given taxon, and applied the same permutation test. To test whether the distance from cohort 1 to cohort 2 is greater than the distance from cohort 1 to cohort 3, we used a permutation test with the test statistic as the difference between these two distances and block-permutation was used for assessing the significance.
Differential abundance analysis
We conducted differential abundance analysis at phylum, family, and genus levels and filtered rare taxa with prevalence less than 20 % to reduce the number of the tests. We fit a generalized mixed-effects model to the taxa count data using the PQL method, assuming a random intercept for each participant to account for within-subject correlation (“glmmPQL” in R “MASS” package). We fitted an overdispersed Poisson to the counts if the zero proportion is less than 25 % and an overdispersed Binomial model (presence/absence) otherwise. For the overdispersed Poisson model, we included the log of library size as an offset to account for variable sequencing depth. In the overdispersed Binomial model, the log of library size was included as a covariate to account for potential dependence of occurrence probability with sequencing depth. We used the winsorized data (97 % upper quantile) to reduce the potential impact of outliers upon the parameter estimates. To improve power to detect differential taxa, which show consistent change in both the uterus and lower tract microbiome, we pooled the uterus and lower tract data and included the sampling site (uterus/lower tract) as a covariate in the model. The same analyses were also repeated for both datasets separately to confirm the source of the identified signals using pooled data. Statistical significance was assessed based on the Wald test. False discovery rate (FDR) control (B-H procedure, “p.adjust” in standard R packages) was used for correcting for multiple testing, and FDR-adjusted p values or q values will be reported. All statistical analyses were performed in R 3.0.2 (R Development Core Team, Vienna, Austria). The receiver operating characteristic (ROC) curve and area under the curve (AUC) were generated using the median of the replicates with the software generated by Johns Hopkins. (http://www.rad.jhmi.edu/jeng/javarad/roc/).