Cervical Microbiota Associated with Higher Grade Cervical Intraepithelial Neoplasia in Women Infected with High-Risk Human Papillomaviruses

The human microbiome is gaining attention as a factor in health and disease, including cancer. The gut microbiome is the best characterized and implicated in carcinogenesis in the gut (1) and liver by translocation of bacterial metabolites or systemic inflammation (2). Microbial communities in other anatomical sites may alter cancer risk (3, 4). Studies have demonstrated that microbiota influence cancer susceptibility and progression by modulating inflammation, inducing oxidative stress, and promoting genomic instability of the host cells (5, 6).As the bacterial density is high in mucous membranes, microbiota may also play a role in the development of precancerous or cancerous lesions of the uterine cervix. High-risk (HR) human papillomaviruses (HPV) cause high-grade cervical intraepithelial neoplasia (CIN 2+), the precursor of invasive cervical cancer (ICC), but only a small proportion of women infected with HR-HPVs develop CIN 2+. It is important to identify the factors that may alter the natural history of HR-HPV infection.

As the microbiome influences acquisition and pathogenicity of viral infections in the gut (7, 8), the same may occur in the genital tract. Bacterial vaginosis (BV), a condition characterized by a decrease in Lactobacillus with a concomitant increase in anaerobic bacteria (e.g., Gardnerella, Prevotella, and Clostridiales; ref. 9), is associated with increased risk of acquisition, reactivation, or delayed clearance of HPV infection (9–11). Previous studies, however, have relied only on Gram stain to assess BV. Vaginal microbiome sequencing by quantitative PCR showed greater microbiome diversity and a higher prevalence of L. gasseri and Gardnerella vaginalis among HPV-positive women (12). Pyrosequencing of 16S rRNA genes revealed a lower proportion of Lactobacillus spp. and identified Fusobacteria, including Sneathia spp., as key bacteria among HPV-positive women (13). A recent report also using 16S rRNA sequencing found that a vaginal microbiome dominated by L. gasseri was associated with increased clearance of HPV (14).

The purpose of the current study was to examine the association between the cervicovaginal microbiome and CIN 2+ among women well characterized for HPV infection and histologically confirmed diagnoses of CIN lesions and other risk factors. We also assessed whether the cervical microbiome was associated with the presence of cervical cells positive for 8-hydroxy-2′-deoxyguanosine [8-OHdG], a marker of oxidative DNA damage.

We evaluated data on 430 women enrolled in two completed studies (R01 CA105448 and R01 CA102489). All patients were diagnosed with abnormal cervical cells by cytology at the Health Departments in Birmingham, AL, and were referred to a colposcopy clinic at the University of Alabama at Birmingham for further examination. Women were 19 to 50 years of age; had no previous history of cervical or other lower genital cancer; no previous hysterectomy, or destructive therapy of the cervix; and were not pregnant. A total of 340 were diagnosed with CIN 2+ [cases, including CIN 2 (n = 208), CIN 3 (n = 132)] and 90 with CIN 1 (non-cases). Both cases and non-cases were positive for HR-HPV (any one of 13 types of HR-HPV, HPV 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, and 68) based on the Roche Diagnostics Linear Array. Pelvic examinations and biopsies were carried out following established protocols. All research procedures were approved by the Institutional Review Board at the University of Alabama at Birmingham.

Age at entry, race, level of education, cigarette smoking, parity, and use of hormonal/oral contraceptives were ascertained via interview with study personnel. Height and weight measurements were obtained using standard protocols and used to calculate the body mass index (BMI). Cervical mucus samples were collected using Merocel ophthalmic sponges (Medtronic Xomed, Inc.) following speculum insertion prior to collecting exfoliated cervical cells. The sponge was placed in the cervical os, with care not to touch the vaginal wall and vulva, allowed to absorb secretions for 60 seconds and placed into a sterile cryovial containing 5 mL of phosphate-buffered saline (PBS) and then placed immediately on ice following collection. Exfoliated cervical cells collected from the transformation zone with a cervical brush were immediately rinsed in 10 mL of PBS and were kept cold. Both sample types were transported to the laboratory on ice within 2 hours of collection. In the laboratory, the sponges suspended in 5 mL of PBS were stored at 80°C. Cervical cell suspensions were separated into two aliquots and centrifuged, and supernatant was removed to obtain cell pellets. One cervical cell pellet used for HPV genotyping was stored in PreservCyt Solution at −20°C. The other cervical cell pellet was fixed by resuspending in 500 μL of 3% formalin and allowed to stand at room temperature for 15 minutes.

The formalin-fixed cervical cell suspension was centrifuged to get a pellet, and supernatant was removed. The cervical cell pellet was again resuspended in PBS and smeared on microscopic slides and allowed to dry out at 37°C for 20 minutes. After drying, slides were placed in 0.01 mol/L citric acid, pH 6.0, and boiled in a pressure cooker set at full power for 10 minutes. After the slides were cooled, they were rinsed in dH₂O and placed in 3.5 N HCl for 15 minutes at room temperature to open the DNA. The slides were rinsed with dH₂O and treated with 3.0% H₂O₂ for 5 minutes to quench endogenous peroxidase activity. Slides were incubated with preimmune goat serum (1%) for 1 hour at room temperature to suppress nonspecific staining and then subsequently incubated with 8-OHdG antibody (catalog NWA-8OHdG, 1:500 concentration) for 1 hour at room temperature. All cervical cells on each slide were counted to derive the percentages of cells positive for 8-OHdG.

Microbial DNA was isolated from cervical mucus samples using the Fecal DNA isolation kit from Zymo Research following the manufacturer's instructions. Once the sample DNA was prepared, PCR was used with unique bar-coded primers to amplify the variable region 4 (V4) region of the 16S rRNA gene to create an amplicon library from individual samples (15, 16).

The PCR product was ∼255 bases from the V4 segment of the 16S rDNA gene. Paired-end sequencing (2×250) was performed on the Illumina MiSeq (15, 16). FASTQ conversion of the raw data files was performed following de-multiplexing using MiSeq reporter. Quality assessment of the FASTQ files was performed using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and then quality filtering was done using the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). High-quality reads (reads where 80% bases have Q score > 20) were selected for analysis, and reads with unknown bases (“N”) were discarded. The remainder of the steps were performed with the Quantitative Insights into Microbial Ecology (QIIME) suite, version 1.8 as described below (15, 17, 18). Chimeric sequences were filtered with UCHIME using the Ribosomal Database Program (RDP) database (19, 20). Sequences were grouped into operational taxonomic units (OTU) using UCLUST at a similarity threshold of 97% (19). The RDP classifier trained using the Greengenes (v13.8) 16S rRNA database (21) was used to make taxonomic assignments for all OTUs at confidence threshold of 80% (22). OTUs whose average abundance was less than 0.005% were removed (23). Multiple sequence alignment was performed with PyNAST v1.2 (24) and a phylogenetic tree was constructed using FastTree v2.1 (25, 26). A total of 8,600,000 sequences clustered to 652 OTUs were obtained after quality filtering and subsampling to 20k reads. Alpha and beta diversity metrics were calculated as implemented in QIIME (19, 27, 28). Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt; ref. 29) was used to predict Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (i.e., functional gene counts) for each sample and to assign inferred gene counts to KEGG pathways.

Differences in demographic and clinical characteristics of participants according to CIN status were assessed by t tests or Wilcoxon rank-sum tests for continuous variables and χ² tests for categorical variables. Differences in the Chao1 index, Shannon index, and Faith's phylogenetic diversity (α-diversity metrics) according to CIN status were tested by the Wilcoxon rank-sum test. Principal coordinates analysis (PCoA) of the weighted and unweighted Unifrac metrics, as well as the Bray-Curtis, Sorensen, Jaccard, and Morisita-Horn indices (β-diversity), were used to visualize differences in microbial community structure according to CIN and 8-OHdG status. Permutational ANOVA of the distance matrices as implemented by the ADONIS function in the R package vegan (30) was used to test for differences in diversity measures. Non-metric dimensional scaling (NMDS) of the Bray-Curtis and Jaccard distance matrices followed by permutational ANOVA was used to assess differences in the metagenome according to CIN and 8-OHdG status. Although differences in measures of diversity provide evidence that microbial communities are different, it is possible for measures of diversity to be similar while the composition of microbial communities is different. We used a Dirichlet multinomial mixture (DMM) model to partition samples into metacommunities (cervical mucus community types, CT; ref. 31). The method provides a model-based approach to metagenomic clustering that recognizes the discrete nature of the data and uses an evidence-based clustering approach to estimate the number and composition of metacommunities. First, each specimen is considered the result of a sampling process whereby each taxa in the underlying microbial community has a specific probability of being detected, such that the likelihood of the observed sample is derived from a multinomial distribution with a parameter vector representing the probability that each sequence read is from a given taxa. Next, the model uses a Bayesian framework and the Dirichlet distribution as the natural conjugate prior to estimate the parameters of the multinomial sampling distribution of a common, underlying “metacommunity.” The novelty of the model lies in extending the Dirichlet prior to a mixture of Dirichlets and provides the means to view the samples as having been generated by a mixture of metacommunities, facilitating the partitioning of samples into clusters. For this analysis, the DMM model as implemented in Mothur v1.35 (32) was fit to an OTU table aggregated to the genus level to improve clustering. The optimal number of partitions is selected by minimizing the Laplace approximation to the negative log model evidence. The model assigns each sample to one of the CTs (based on its OTU composition) and estimates the weight (i.e., relative frequency, π) and the precision (i.e., a measure inversely related to the degree of dispersion with respect to OTU composition, θ) of each CT. We reviewed the OTU composition of the CTs identified by the DMM model to assess their biological significance, and selected one CT as the referent group, comparing other partitions to the referent. Odds ratios (OR) and 95% confidence intervals (CI) for CIN2+ and high 8-OHdG status according to cervical mucus CT were calculated using logistic regression. Multivariable models included adjustment for age, BMI, race (non-Hispanic black/non-Hispanic white), education (high school graduate yes/no), parity, hormonal contraceptive use (yes/no), smoking status (current smoker yes/no), and CIN status where appropriate. The indirect effect of CT on CIN 2+ mediated by 8-OHdG was assessed by decomposing the total effect into direct and indirect components and estimating the indirect effect and 95% CI using the PROCESS macro for SAS (33). Linear Discriminant Analysis Effect Size (LEfSe; ref. 34) and DESeq2 (35) were used to detect differences in OTUs according to CIN and 8-OHdG status. LEfSe was performed on the subsampled OTU table using the default values of α = 0.05 and LDA = 2.0. The negative binomial regression model implemented in DESEq2 was fit to the raw counts adjusting for factors contributing to CIN risk (covariates described above). Log₂fold changes for OTUs with false discovery rate (FDR) corrected P < 0.05 are reported. LEfSe was also used to assess enrichment in imputed genes and gene pathways according to CIN and 8-OHdG status.

Demographic and clinical characteristics of participants according to CIN status are provided in Table 1. The mean (±SD) age of participants was 26.1 (±4.7) years. The majority were non-Hispanic black (53%), had more than a high school education (81%), were not current smokers (53%), and were current users of hormonal contraceptives (84%). Mean BMI and parity were 28.0 (±8.0) and 1.3 (±1.0), respectively. CIN 2+ cases were younger, had a lower BMI, and higher rates of HPV 16 infection than CIN 1, non-cases (P < 0.05). The proportion of women classified to cervical mucosa CT, and with HPV 18 infection, did not differ. The proportion of 8-OHdG-positive cells was 43% in CIN 2+ cases and 35% in non-cases (P = 0.062).

The DMM model identified four CTs with weights π = (0.15, 0.41, 0.16, 0.28) and precision θ = (7.7, 89.8, 283.0, 55.6). Partition 1 had the lowest precision and was composed of diverse taxa, including unclassified Lactobacillus, L. iners, Bifidobacteriaceae, Clostridiales, and Allobaculum (Fig. 1). Partition 2 was the most prevalent, while Lactobacillus dominated, uniquely characterized by the high diversity and proportion of rare taxa. Partition 3 was the least varied and dominated by unclassified Lactobacillus (43%) and L. iners (40%) organisms. The percentages of sequence reads mapping to L. iners in partitions 1, 2, and 4 were 15%, 8%, and 6%, respectively. Partition 4 was the only CT lacking Lactobacillus dominance and primarily composed of Bifidobacteriaceae, Prevotella, Sneathia, and Megasphaera. Based on a priori considerations about the biological potential of these CTs, we considered partition 1 to be “low risk” and selected it as the referent group, and compared the other three partitions to the referent. Overall, Lactobacillus genera were dominant (represented >5% of all sequence reads) in 55% of samples. The proportion of African American women (67%) was higher in partition 4 compared with Caucasian American women (33%). Age, educational status, BMI, smoking status, parity, use of hormonal contraceptives, and the percentage of women positive for HPV 16 or 18 did not differ by partition (data not shown).

Alpha-diversity metrics did not differ in cases and non-cases (P > 0.39). Ordinations based on beta-diversity distance matrices and permutational ANOVA (e.g., ADONIS) did not support differences in community structure (Fig. 2). Similar findings were obtained when comparing CIN 3+ with CIN 1 (data not shown). Whereas overall CT partitioning did not significantly differ between CIN 2+ and CIN 1 (Table 1, P = 0.208), pairwise comparisons between “high-risk” CTs and the referent were significant.

In unadjusted models, samples clustered to partition 3 had a 2.5-fold higher odds [OR = 2.48 (1.01; 6.07)] of being diagnosed with CIN 2+ when compared with samples clustered to partition 1 (Table 2). The OR increased to 3.5 (1.27; 9.55) in models adjusted for factors contributing to CIN risk. Compared with partition 1, ORs for samples clustered to partition 2 were 1.63 (0.83; 3.21) and 1.84 (0.86; 3.93) for unadjusted and multivariable models, respectively. ORs for CIN 3+ were similar to those obtained for CIN 2+, but failed to achieve statistical significance. The odds of CIN 2+ and CIN 3+ for partition 4 were similar to those for partition 1. In the subset of samples (n = 156) for which 8-OHdG was assayed the odds ratio for the indirect effect of partition 3 operating through 8-OHdG on CIN 2+ was 1.12 (0.88; 2.44; ∼12.5% of total effect mediated; data not shown). Thus, a relatively small proportion of the association between CIN 2+ diagnosis and partition 3 appears to be mediated by oxidative DNA damage. The indirect effect of 8-OHdG was not testable for the other CT partitions.

LEfSe analyses assessing taxa discriminating on CIN status are shown in Fig. 2 (top 25) and Supplementary Table S1. Samples from CIN 2+ were enriched in Lactobacillaceae, Lactobacillus, L. reuteri and several sub–genus level Lactobacillus OTUs (effect size >2.0; P < 0.05). Bacterial families, including Bacteroidaceae, Porphyromonadacae, and Coxiellaceae; genera, including Bacteroides, Parabacteroides, Rickettsiella, and RFN20; and numerous OTUs (n = 20) were also enriched in CIN 2+ samples. Two OTUs classified to Ruminococcaceae and Clostridiales were enriched in samples obtained from CIN1 non-cases. Similar findings for OTUs were obtained from negative binomial regression after accounting for differences in sequencing depth and adjustment for factors contributing to CIN risk (Supplementary Table S2).

PICRUSt was used to infer the functional metagenomic content of each sample and to classify the inferred genes to KEGG pathways. NMDS of the Bray-Curtis and Jaccard distance matrices and permutational ANOVA did not support differences in the inferred metagenome by CIN status (P = 0.20; Supplementary Fig. S1). LEfSe analysis of the inferred metagenome identified five KO identifiers enriched in women diagnosed with CIN 1, including K02003 (ATP-binding protein), K09687 (antibiotic transport system ATP-binding protein), K02004 (ABC transport system permease protein), K06147 (ATP-binding cassette), and K07133 (unclassified; data not shown). Analysis of inferred functional pathways identified glycosyltransferase enzymes and pyruvate metabolism enriched in CIN 2+ samples and enrichment in alanine, aspartate, and glutamate metabolism; cysteine and methionine metabolism; and phenylalanine, tyrosine, and tryptophan biosynthesis in CIN1 (effect size >2.0; P < 0.05; Fig. 2).

The percentage of 8-OHdG-positive cervical cells was not associated with microbial diversity when examined as a continuous variable or dichotomized at the median (median = 38.8%; P < 0.23 all models; data not shown). Compared with partition 1, odds ratio for samples clustered to partition 2 and 3 were 2.23 (0.82; 6.70) and 2.61 (0.85; 8.03) in multivariable models, respectively (Table 2). Enrichment in Betaproteobacteria, Burkholderiales, and Rhizobiales; several bacterial families including Psuedomonadaceae, Oxalobacteraceae, Rhizobiacaeae; and several genera including Pseudomonas, Argobacterium, limnohabitans, and Phenylobacterium was found in samples from women with above median values for 8-OHdG (effect size >2.0; P < 0.05; Supplementary Fig. S2). Samples from women with low 8-OHdG levels were enriched in Clostridiaceae. Analysis of imputed metagenomes identified two genes enriched in non-cases including K08884 (serine/threonine protein kinase) and K01897 (long-chain acyl-CoA synthetase). Analysis of the imputed KEGG pathways identified enrichment in two component system, signal transduction, lipid metabolism, transcription factors, secretion system, metabolism of other amino acids, and aminobenzoate degradation in samples from women with a greater percentage of 8-OHdG–positive cells (Supplementary Fig. S2). Enrichment in alanine, aspartate, and glutamate metabolism; terpenoid backbone synthesis; thiamin metabolism; one carbon pool by folate; cysteine and methionine metabolism; translation proteins; amino acid–related enzymes; genetic information processing; and energy metabolism was found for samples from women below the median percentage cells positive for 8-OHdG.

We studied a well-characterized group of women to evaluate the association between cervical microbiota and CIN 2+ and to assess whether cervical microbiota are associated with oxidative DNA damage. Measures of alpha- or beta-diversity were not associated with either CIN severity or oxidative DNA damage. A cervical mucosal community type dominated by L. iners and unclassified Lactobacillus spp. was associated with CIN 2+. Sequence reads mapping to Lactobacillaceae, Lactobacillus, L. reuteri, and several sub–genus level Lactobacillus OTUs were also associated with CIN 2+. We found little evidence that DNA oxidative damage mediates the effect of the microbiome on the natural history of HPV infection and CIN severity.

Healthy vaginal flora is described as Lactobacillus dominated with low bacterial diversity and accompanying low pH, while “unhealthy” flora is often characterized by high diversity bacterial populations with increased anaerobic bacteria and lactobacilli present at lower levels (36). High diversity communities are associated with pro-inflammatory cytokine concentrations (37). Recent data from diverse human populations have suggested that this characterization is overly simplistic. Anahtar and colleagues recently documented that in a cohort of young, healthy South African women, only 37% had a lactobacillus-dominated vaginal microbiota (37). Although recent research has clarified host-commensal dynamics in the gastrointestinal tract and to a certain extent in the vaginal microbiome in relation to pregnancy complications and pelvic inflammatory diseases (38–40), much less is understood about the cervical microbiome in relation to cervical cancer risk. This gap must be filled if we want to explore the potential for microbiota-based interventions to prevent cervical cancer.

Microbial communities in the cervical mucosa have not been well studied, and comparing our findings with previous research is problematic. Because of the close proximity of the organs, however, given that the cervicovaginal microbiome plays a role in the acquisition and natural history of HPV infections, it is plausible that the microbial communities of the cervical mucosa modify the risk of developing HR-HPV–associated CIN 2+. No comprehensive studies have evaluated this hypothesis. A study compared Korean women diagnosed with CIN 1 (n = 55) and CIN 2+ (n = 15) with those with normal cytology or ASCUS (n = 50) and reported that a cervical cell microbiome with a predominance of A. vaginae, G. vaginalis, and L. iners and a concomitant paucity of L. crispatus, in combination with oncogenic HPVs, may be a risk factor for cervical neoplasia (41). The study is of limited value because CIN 1 is not a “true” pre-neoplastic lesion, and some of the women diagnosed with ASCUS are likely to have underlying CIN 2+ lesions. Further, the results reported were not adjusted by other risk factors.

We observed that a cervical mucosa CT dominated by L. iners and unclassified Lactobacillus was associated with CIN 2+. Smith and colleagues also observed that cervical microbiota clustered in six distinct CTs, one of which was characterized by a predominance of L. iners (42). L. iners, discovered by Falsen and colleagues (43), has escaped attention because it grows only on blood agar. L. iners–dominated samples have lower concentrations of D-lactic acid, which is protective against infection by pathogens (44). L. iners has also been reported to induce IL8 secretion in vivo and thereby potentially moderate pro-inflammatory activity in the cervix, which may influence CIN progression (36). We did not observe an association between L. iners and CIN 2+ when examining single OTUs, but did so when aggregating taxa to Lactobacillaceae and Lactobacillus, as well as for L. reuteri and several other sub–genus level Lactobacillus OTUs. Thus, if the association between L. iners and CIN 2+ is real, our findings suggest that the broader community structure may influence the effect of L. iners on CIN 2+ risk. Studies incorporating whole-genome shotgun sequencing, finer taxonomic resolution, and metatranscriptomics may provide greater insight into the nature of this association.

The limited precision of the ORs estimated for the cervical mucosa CTs may reflect the inherent challenge of attempting to reduce a high-dimensional community structure to a few clusters. Figure 1 shows clear clustering of samples into CTs, but the large within-CT spread suggests residual variation in community structure. Potential contamination of specimens during sampling, processing, or sequencing may have also reduced our ability to detect associations with CIN 2+. Primer selection and library preparation choices also represent a potential source of bias (45).

On the other hand, clustering may capture the interactions of specific microbes operating in complex community structures that may not be well represented using general diversity measures based on underlying distance matrices, and strengthens the power of our analyses. There is a growing body of literature in the field of microbial ecology pertaining to the use of Dirichlet multinomial models (DMM) to identify naturally occurring clusters of samples based on microbial composition. These methods have been used in numerous studies to cluster samples into "community types" or "enterotypes" (46–49) and readily implemented in commonly used software programs, including Mothur and R (16, 50). DMM clustering provides information that is distinctly different from ordinations performed on dissimilarity matrices, as it is performed directly on read counts (i.e., do not rely on a continuous approximation to the discrete nature of the data), takes into account sequence variation, does not rely on an underlying distance matrix, and is not limited to defining a fixed number of clusters. The DMM described by Holmes and colleagues (31) uses a Bayesian framework to assign samples to the community type for which it has the highest posterior probability of membership with the optimal number of clusters based on minimizing model error. The relative abundance of taxa in each cluster can then be used to facilitate interpretation of major organisms comprising the group. The DMM approach used here has been shown to reveal important CTs when applied to data from the Human Microbiome Project and superior performance to alternative clustering methods (51).

Whereas previous studies have suggested that BV may alter the natural history of HPV infection, in this study, a “BV-like” CT (partition 4) was not associated with CIN 2+. As we did not use conventional BV diagnostic tests, it is possible that women in partition 4 did not have clinical BV. The cervical microbiome of CIN 2+ also tended to be enriched in potentially beneficial microbes such as L. reuteri (52, 53), but also microbes found in the intestinal tract such as Bacteroides, which have been associated with colorectal cancer (54) and cause significant pathology when they escape the intestinal tract (55). In contrast, samples of non-cases were enriched with Ruminococcaceae, one of the most abundant gut Firmicute families associated with good gut health (56, 57) and Clostridia, common gut microbiota involved in the maintenance of gut function (58). In healthy women, microbiota in the gut and cervicovaginal area are confined to the mucosa by a multilevel barrier system that includes a mucus layer, secretion of soluble immune mediators, and intact epithelium with tight junctions. Failure of this barrier may lead to translocation of microbiota (59), possibly inducing inflammation, oxidative DNA damage, and ultimately the development of CIN 2+. Mucin-degrading enzymes secreted by these microbes may destroy the mucosal barrier and increase susceptibility to HPV (60). Many bacterial species colonize both the gastrointestinal and the reproductive tract, and the rectum serves as a reservoir of organisms that colonize the cervicovaginal area (61–63).

Whereas we did not find a statistically significant association between the percentage of cervical cells positive for 8-OHdG and microbial diversity or community structure, OR estimates were >2. We note that the highest percentage of women with above median 8-OHdG was found in partition 3, which was associated with CIN 2+. Thus, a biologically plausible association may exist despite our inconclusive findings. Oxidative DNA damage is an inevitable consequence of cellular metabolism and is increased following toxic insult. 8-OHdG is an established cancer biomarker of oxidative stress (64): elevated levels of 8-OHdG were reported in invasive ductal carcinoma (65), primary breast cancer (66), and colorectal cancer (67) compared with normal cells of the corresponding organ, and in precancerous lesions of the cervix compared with normal cervical epithelial cells (68, 69). Assessing the impact of the cervical microbiota on 8-OHdG levels is important to understand the role of microbial community composition and oxidative stress in the progression of HPV-related cervical lesions. The microbiome influences reactive oxygen species (ROS) production in the gut epithelium (70) and chronically high ROS levels can outpace DNA repair mechanisms, leading to DNA damage (71). In vitro studies have shown that ROS production is associated with 8-OHdG formation and increased DNA damage (72). ROS induce double-strand breaks in the host DNA and in the HPV episome, facilitating integration of the virus and carcinogenesis (63). There is accumulating evidence that microbial metabolites may induce cancer via DNA damage (73).

Glycosyltransferase enzymes and pyruvate metabolism were enriched in CIN 2+ samples. Although aberrant pyruvate metabolism (74) and glycosyltransferases (75) play a role in carcinogenesis, the impact of bacterial pyruvate metabolism on CIN is unclear. Some of the genes identified in non-cases and in women with lower 8-OHdG, such as methionine metabolism, play an important role in the pathway linking folate/vitamin B12 metabolism with DNA methylation and DNA repair: we have previously reported that a higher degree of DNA methylation is inversely associated with CIN 2+ (76). Because the gut microbiome supplies its host with up to a third of the recommended intake of B vitamins (77, 78), and because the rectum is a reservoir for organisms that colonize the vagina, our previous observations could be explained by differences in gut microbial profiles by CIN status. The folate biosynthesis pathway is present in nearly all Bacteroidetes genomes, as well as in most Fusobacteria and Proteobacteria while vitamin B12 synthesis of cobalamin is present in all Fusobacteria, but is rare in Actinobacteria and Proteobacteria (79). Like colonocytes, cervical epithelial cells contain transporters for these micronutrients (80, 81). Thus, cervical microbiota may influence the natural history of HPV infection and CIN by altering the status of these micronutrients in the cervical environment.

Some of the genes identified in women with lower 8-OHdG such as one-carbon pool by folate and methionine metabolism play an important role in the folate metabolic pathway, which contributes to DNA synthesis and repair. It is biologically plausible that pathways such as signal transduction are enriched in women with a greater percentage cells positive for 8-OHdG because DNA damage triggers a signal transduction pathway that appears to be a kinase cascade (82).

In summary, our study did not show an association between microbiome diversity and CIN severity or oxidative DNA damage, but provided suggestive evidence that a cervical microbiome characterized by a predominance of Lactobacillus and L. iners is associated with CIN 2+ in women infected with HR-HPVs. Future studies are needed to confirm our 16S rRNA sequencing results and to resolve associations for specific organisms and at finer taxonomic levels using whole-genome shotgun sequencing together with measures of cervical microbiome-induced molecular pathways that may alter the risk of CIN 2+. Our data suggest that only a small proportion of the effect of the microbiome on CIN may be mediated through oxidative DNA damage.

No potential conflicts of interest were disclosed.

Conception and design: C.J. Piyathilake, M. Macaluso

Development of methodology: C.J. Piyathilake, C.D. Morrow

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C.J. Piyathilake

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C.J. Piyathilake, N.J. Ollberding, R. Kumar, M. Macaluso, C.D. Morrow

Writing, review, and/or revision of the manuscript: C.J. Piyathilake, N.J. Ollberding, R. Kumar, M. Macaluso, R.D. Alvarez, C.D. Morrow

Study supervision: C.J. Piyathilake

The authors acknowledge the excellent technical assistance provided by Peter Eipers, Michelle Moses Chambers, and Suguna Badiga.

This work was supported by grants R21 CA188066 (C. Piyathilake), R01 CA102489 (C. Piyathilake), and R01 CA105448 (C. Piyathilake) of the National Cancer Institute. The following are acknowledged for their support of the Microbiome Resource at the University of Alabama at Birmingham: School of Medicine, Comprehensive Cancer Center (P30AR050948), Center for AIDS Research (5P30AI027767) and Center for Clinical Translational Science (UL1TR000165).The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.