Cervicovaginal microbiome and natural history of HPV in a longitudinal study

1Persistent cervical infections by high-risk (HR) human papillomavirus (HPV) cause virtually all cervical cancers and their immediate precursor lesions. Most sexually active women have been infected with HPV at some point in their lives and in the vast majority the infection is cleared within a few months. However, a subset of women develop a persistent HPV infection that places them at high risk for cervical precancer and cancer. Fig 1 illustrates this paradigm canonically referred to as HPV natural history.

2Non-viral factors (HPV co-factors) associated with the outcomes of HR-HPV infections have not been fully elucidated. While smoking, hormonal contraceptive use, and parity are associated with developing precancer and cancer; systemic and local immune responses are thought to be important for clearance and control of infection (persistence vs. clearance). In addition, specific host immune regulatory alleles (e.g., human leukocyte antigen) are associated with risk of cervical cancer.

3The local, cervical microenvironment, including the microbiome, may also influence the natural history of HPV infection. Other studies have recently implicated the microbiome’s role in the natural history of other viral infections and a variety of cancers. The cervicovaginal microbiome (CVM) is of particular interested because it has been well characterized and specific features have been associated with gynecologic disease and reproductive health. The CVM has been categorized into community state types (CSTs) generally defined by a dominance of a specific Lactobacillus species (Lactobacillus crispatus, Lactobacillus iners, Lactobacillus gasseri or Lactobacillus jensenii), or a state of polymicrobialism. Transitions from Lactobacillus dominated CSTs have been linked to detrimental health outcomes including elevated risks for sexually transmitted infections, as well as higher incidences of preterm births.

4An association between increased CVM diversity and prevalence of HR-HPV infection and/or cervical abnormalities (vs. HPV negative) has been reported in several studies. Higher abundance of L. crispatus has been shown to be associated with lower HPV prevalence and increased detection of normal cytology. Long-term use of vaginal probiotics with Lactobacillus spp. has been associated with increase clearance of HPV compared to short-term use. However, evidence is conflicting on the association of CVM diversity and the severity of cervical neoplasia. Additionally, most studies looking at the natural history of HPV and the microbiome are cross-sectional and therefore lack the ability to draw potential causal links.

5For the current study, we leveraged longitudinal data and specimens from the placebo arm of a large randomized HPV vaccine trial to examine the impact of the CVM on the natural history of incident HR-HPV infections to study: 1) progression to cervical precancer, 2) viral persistence, and 3) viral clearance.

6A total of 273 women with an incident HR-HPV infection were included in the analyses, of whom 266 had a second sample with an average sampling interval of 1.5±0.9 years. Table 1 presents the sample demographic information and summarizes the bacterial and fungal sequencing results after taxonomic assignment for each infection outcome at baseline. There were no significant differences between groups in age (p = 0.13), 16S rRNA gene OTU clustered read counts (p = 0.33), or ITS1 OTU clustered read counts (p = 0.53).

7Fig 2 summarizes the bacterial Shannon diversity measures of the three ordered categorical HR-HPV outcomes (see Fig 1) within the CVT cohort at V1 and V2. Alpha diversity analysis did not reveal any significant differences at V1 in terms of bacterial Shannon alpha diversity (trend p = 0.52). At V2 there was a significant trend of rising diversity based on the Shannon diversity index (trend p = 0.024).

8To evaluate the overall structure of the cervicovaginal microbiome, we performed hierarchical clustering on all available samples (n = 539) (Fig 3A). This analysis revealed four distinct bacterial community state types (CSTs). Two CSTs were dominated by species of the genus Lactobacillus (Lactobacillus iners, 143/539, 26.5% and Lactobacillus crispatus, 83/539, 15.4%), one CST by Gardnerella vaginalis (94/539, 17.4%), and the other CST did not contain a major group but had a highly diverse microbiome (219/539, 40.6%).

9Fig 3B shows results of hierarchical clustering based on the detected fungal species. Candida albicans was the dominant fungal taxa. In terms of fungal clusters, there appears to be a single clade dominated by C. albicans (43/208, 20.7%), one dominated by an unidentified fungal species (12/208, 5.8%) and one that is composed of a diverse fungal community (153/208, 73.6%).

10Analysis of bacterial taxonomic categories of the microbiome associated with persistence/progression vs. clearance revealed a total of 24 taxa that were significant (LDA>2.0) at V1. G. vaginalis was the bacterial species with the highest positive correlation to progression (Fig 4A), while L. iners was the most positively associated taxon with clearance based on relative abundance. Amongst V2 samples there were a total of 13 significant taxa identified (three taxa associated with clearance and 10 with progression) (Fig 4B). At V2, bacteria that are commonly associated with bacterial vaginosis, such as Prevotella amnii and Anaerococcus prevotii, were significantly correlated with progression. We used a generalized linear model (GLM) to validate LEfSe biomarkers with adjustments for key covariates (age, smoking status, HPV16 and visit CST) (S1 Table).

11Fungal biomarker discovery revealed five fungal taxa associated with HPV progression (S1A Fig). The average of the five combined fungal taxa were detected at rates of 0.59%, 3.55% and 3.76% for the clearance, persistence and progression outcomes, respectively (p = 0.0080, S1B Fig). However in the adjusted GLM analysis, none of the fungal biomarkers were determined to be significant at either visit (S2 Table).

12To evaluate whether some common function of bacteria might be associated with HR-HPV outcomes, we used a functional analysis of gene groups imputed from the microbiome data as described in the methods. This analysis revealed 8 functional pathways at KEGG Level 2 significantly associated with progression to CIN2+ (S3 Table). Of the 8 identified pathways, only 2 had a mean read coverage of >1% and were considered for further analysis. The two identified pathways were “Xenobiotics Biodegradation and Metabolism” pathway, which was positively associated with progression (p = 0.0020) and the Cell Motility pathway, which was negatively associated with progression (p = 0.019). These two pathways were significantly correlated (Pearson correlation = -0.80, S2 Fig), and we chose to use Cell Motility in multivariate modeling since it produced a more stable GLM estimate (S3 Table).

13To investigate the contribution of components of the microbiome over time, we used a GLM in order to adjust for known covariates of CIN2+ progression that may influence the relationship of the CVM and progression to CIN2+ (e.g. age, smoking, HPV16 and CST). We utilized a GLM with a binary outcome (clearance/progression) and the significant microbial features identified in preceding sections as predictors. Specifically, we used the abundance of Gardnerella at V1, the abundance of Lactobacillus at V1, the Observed fungal species diversity at V1, the imputed Cell Motility pathway at V1 and the microbial diversity at V2. The model was adjusted for age, CST, smoking and HPV16 infection status. Fig 5 shows the model estimates of the resulting GLM analysis. The multivariate analysis revealed a significant protective effect of V1 Lactobacillus (genus) abundance, OR = 0.41 (0.22–0.79), V1 fungal species diversity, OR = 0.90 (0.82–1.00) and imputed V1 Cell Motility pathway OR = 0.75 (0.62–0.92). In addition, the model showed that the V2 microbial diversity was a significant risk factor for CIN2+ progression, OR = 1.17 (1.02–1.29).

14Following the multivariate analysis, we wanted to explore the reason for V1 Gardnerella being insignificant despite being the top microbial risk factor in differential abundance in all previous analyses (Fig 4, S3 Fig and S1 Table). Thus, we performed a mediation analysis to determine if V1 Gardnerella was acting by inducing the elevated diversity at V2 (Fig 6A and 6B). This analysis showed that after adjustment for V1 Gardnerella there was a significant association of V2 Shannon diversity with progression to CIN2+, p = 0.04. The Average Direct Effect (ADE) showed that V1 Gardnerella wasn’t significant after adjustment for the V2 Shannon diversity, p = 0.23 supporting our mediation hypothesis (Fig 6A).

15Power for detection of the effects of each microbiome component was performed using the lmSupport package. Specifically each of the microbiome components (i.e., V1 Gardnerella, V1 Lactobacillus, V1 Fungal Observed OTUs, V1 Cell Motility, and V2 Shannon diversity) was tested separately to assess power with adjustment for age, CST, HPV16 status and smoking. S4 Table shows the results of the power calculation. Given large effect sizes for V1 Lactobacillus, V1 Fungal Observed OTUs, V1 Cell Motility, and V2 Shannon diversity we calculated that these could be detected with a power >98%, while the V1 Gardnerella had a power of 82% at the 0.05 alpha level.

16Previous cross-sectional studies analyzing the association between the cervicovaginal microbiome and HPV infection outcomes have consistently identified Gardnerella as a key biomarker associated with CIN2+. This has been reported in studies that utilized both next-generation sequencing and other methods of microbiome analyses. We present data that Gardnerella is in fact associated with CIN2+ lesions, but rather than directly causing the CIN2+ lesion, appears to induce a higher diversity CVM over time as measured at V2, which in turn mediates the observed effect of Gardnerella in HR-HPV disease progression. Although it is not clear how a state of polymicrobialism in the presence of a persistent HR-HPV infection leads to the development of epithelial dysplasia, recent studies on the microbiome’s role in other cancers suggests that the answer lies in the establishment of a microbial microenvironment, perhaps a biofilm. For instance, it has been shown that certain cancers (e.g., colorectal cancer and prostate cancer) have distinct microbial communities at the tumor site that are associated with tumor development. In fact, other data from our lab using cervical biopsy tissue samples indicates that there are distinct microbial differences as cervical cancer progresses to advanced FIGO stages (manuscript in preparation). This idea is further supported by data indicating that cervical precancerous lesions that regress, compared to those that progress to cancer, harbor a different immune microenvironment. The local interplay between the microbiome and the local host immune system may be important to understanding the progression of HR-HPV infection to cervical cancer.

17The protective microbial biomarkers identified at V1 also suggest an association of the microbiome and host innate and acquired immunity in progression to CIN2+. Specifically, the protective effect of bacterial Cell Motility may be due to the known phenomena of bacterial flagella activating host immunity. This local activation may facilitate the innate immune system’s ability to clear an active HPV infection. Such stimulation may be critical in HPV control since cervical lesions have been shown to be associated with local immunosuppression through the reduction of factors such as IL-17. Despite the presence of studies to support this conjecture it should be noted that this type of immune activation needs to be confirmed with rigorous experimental precision.

18Gardnerella, as discussed above, continuously emerges as a risk factor for CIN2+ development and progression. Based on our findings and published data, the association may be tied to the ability of Gardnerella to be immunosuppressive in the cervicovaginal region. Whereas, it seems that the presence of commensal bacteria (e.g. Lactobacillus) with the ability to stimulate a local immune response may be contributing factors to the clearance of incident HR-HPV infections. Moreover, the presence of bacteria with immunosuppressive attributes, such as Gardnerella, may promote viral persistence and progression.

19Alternatively, there may be other explanations for the observed associations between the cervicovaginal microbiome and HPV’s natural history. One possible explanation is a host developed or inherited immune deficiency that is a common cofactor for both cervical cancer progression and microbial diversity. For example, elevation of a particular inflammatory cytokine may be both necessary for successful tumor growth and be a causal factor in increasing vaginal microbial diversity. Such a factor may also explain the consistent identification of Gardnerella, which is commonly identified as a biomarker for increased diversity in the CVM and a risk factor for CIN2+.

20We have identified distinct microbial biomarkers that either protect, or promote the progression of a HR-HPV infection to CIN2+ lesions. In the context of what is known about the cervicovaginal microbiome, it may be that these factors act to suppress (in the case of progression) or activate (in the case of clearance) a localized immune response, which in turn influences the natural history of a HR-HPV infection. However, additional prospective studies are needed to establish a causal link between the cervicovaginal microbiome, the immune system and the natural history of HPV. Nevertheless, our results suggest a marker for identifying women with persistent HR-HPV infection at risk for progression by monitoring the presence of Gardnerella and subsequent elevation in microbial diversity. If future studies support a causal role of the cervicovaginal microbiome and disease progression, then it might be possible to manipulate the CVM in a manner to activate a local immune response. It is possible that HPV vaccination might influence the CVM and future research will be needed to evaluate such potential changes.

21The strength of this study includes the prospective design and availability of a longitudinal cohort. In addition we used advanced epidemiological methods in a novel way to investigate potential causative factors in cervical intraepithelial neoplasia. Potential weaknesses in this study include the relatively small sample size, homogeneity of the population and the use of only two time points.

22In summary, through the use of longitudinal samples from the CVT cohort we investigated and identified key features of the cervicovaginal microbiome potentially associated with progression of HR-HPV infection (e.g., Gardnerella and subsequent increase vaginal microbial diversity). Additional studies are required to validate the model proposed in this report.

23The study of cervicovaginal microbiome and HR-HPV natural history was a nested analysis within the previously reported CVT (clinical trials registration NCT00128661). Written informed consent was obtained from all participants in CVT. The trial protocol can be obtained from the original trial publication.

24All CVT participants were adult women between the ages of 18–25 years. All participants were shown a video describing the study design and were then required to provide written consent to continue participating in the trial. Institutional review board approval was obtained for the informed consent forms at both the NCI and in Costa Rica. Registered with Clinicaltrials.gov NCT00128661

25Subjects for this nested study were selected from the placebo arm of a community-based clinical trial of the HPV 16/18 vaccine in Costa Rica that had enrolled women 18 to 25 years of age in 2004–2005. Women with an incident HR-HPV infection (HPV16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, or 59) were selected for analysis. Incident infections were classified based on outcomes from the well-established model of HR-HPV natural history including outcomes of clearance, persistence and progression. Specifically, outcomes related to the incident HR-HPV infection included women who developed a CIN2 or CIN3 (CIN2+) lesion (progression), women with an infection for 2 or more years with the same type in the absence of a CIN2+ diagnosis (persistent), or women who cleared their incident HR-HPV infections within 1 year (clearance). This analysis included 273 women of whom all had available samples at V1 (first visit positive for the studied HPV type) and 266 who had a second sample at V2 (for persistent, visit that was positive for the same type and at least 305 days after V1; progression, closest visit before diagnosis of CIN2+; clearance, following visit that was negative for that type); all had clinical follow-up data. Seven women, one with clearance, and six with persistence either did not have an available V2 sample or the sample failed in lab testing.

26DNA samples were shipped to the Burk Lab on dry ice where the microbiome analysis was performed. DNA had been extracted from cervical brush samples by DDL Diagnostic Laboratory (Voorberg, The Netherlands) where they had been tested for HPV as previously described. Laboratory procedures for the microbiome analyses were performed within a hood (AirClean Systems, Creedmoor, NC) in an isolated room to minimize environmental contamination and water-blank negative controls were used throughout the testing.

27Bacterial DNA was amplified using barcoded-primers 16SV4_515F (GTGYCAGCMGCCGCGGTA) and 16SV4_806R (GGACTACHVGGGTWTCTAAT) that amplify the V4 variable region of the 16S rRNA gene. This region has been demonstrated to accurately amplify and resolve vaginal bacteria. PCR reactions were performed with 17.75 μl of nuclease-free PCR-grade water (Lonza, Rockland, ME), 2.5 μl of 10X Buffer w/ MgCl2 (Affymetrix, Santa Clara, CA), 1 μl of MgCl2 (25 mM, Affymetrix, Santa Clara, California, USA), 0.5 μl of dNTPs (10 mM, Roche, Basel, Switzerland), 0.25 μl of AmpliTaq Gold DNA Polymerase (5 U/μl, Applied Biosystems, Foster City, California), 0.5 μl of HotStart-IT FideliTaq (2.5 U/μl, Affymetrix, Santa Clara, CA), 1 μl of each primer (5 μM), and 0.5 μl of sample DNA. Thermal cycling conditions consisted of initial denaturation at 95°C for 5 min, followed by 15 cycles at 95°C for 1 min, 55°C for 1 min, and 68°C for 1 min, followed by 15 cycles at 95°C for 1 min, 60°C for 1 min, and 68°C for 1 min, and a final extension for 10 min at 68°C on a GeneAmp PCR System 9700 (Applied Biosystems, Foster City, CA).

28The fungal DNA ITS1 region was amplified using barcoded-primers ITS1_48F (ACACACCGCCCGTCGCTACT) and ITS1_217R (TTTCGCTGCGTTCTTCATCG) as previously described. PCR reactions were performed with 8.25 μl of nuclease-free PCR-grade water (Lonza), 2.5 μl of 10X Buffer w/ MgCl2 (Affymetrix), 1 μl of MgCl2 (25 mM, Affymetrix), 0.5 μl of dNTPs (10 mM, Roche), 0.25 μl of AmpliTaq Gold DNA Polymerase (5 U/μl, Applied Biosystems), 0.5 μl of HotStart-IT FideliTaq (2.5 U/μl, Affymetrix), 1μl of each primer (5 μM), and 10 μl of sample DNA. Thermal cycling conditions consisted of initial denaturation of 95°C for 3 min, followed by 35 cycles at 95°C for 30 sec, 55°C for 30 sec, and 68°C for 2 min, followed by a final extension for 10 min at 68°C on a GeneAmp PCR System 9700 (Applied Biosystems).

29For both amplicon experiments, 20 negative controls were randomly mixed amongst samples. Negative controls were created using nuclease-free PCR-grade water (Lonza) as described above instead of extracted DNA.

30Barcoded-PCR products were combined for each amplicon type and the DNA fragments (~356 bp for 16S V4 and 400–600 for ITS1) were isolated by gel purification using a QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). Purified PCR products were quantified using a Qubit 2.0 Fluorometic High Sensitivity dsDNA Assay (Life Technologies, Carlsbad, CA) prior to library construction using a KAPA LTP Library Preparation Kit (Kapa Biosystems, Wilmington, MA). Size integrity of the amplicons was validated with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). High-throughput amplicon sequencing of 2x300 paired-end reads was conducted on an Illumina MiSeq (Illumina, San Diego, CA).

31Illumina reads were trimmed to remove bases that had a PHRED score of <25 using prinseq-lite V0.0.4. Quality trimmed reads were then demultiplexed using Novobarcode. Paired-end reads were joined using PANDAseq with default settings. The merged reads were processed through the VSEARCH quality-filtering pipeline to dereplicate the sequences, reduce noise and remove chimeric reads.

32For bacterial 16S V4 rRNA gene reads, closed-reference OTU picking was performed using VSEARCH with a custom database that contained sequences from the GreenGenes database the Human Oral Microbiome Database (HOMD) and cervicovaginal microbiome 16S reference sequences retrieved from NCBI. Representative sequences were aligned using PyNAST and taxonomy was assigned using VSEARCH.

33PICRUSt was used to impute microbial functional gene content and to collapse identified genes into functional pathways. Pathways that were associated with HR-HPV clearance were identified through the use of a generalized linear model (GLM) based on statistical significance (p<0.05) and relative abundance (1% or higher across all reads).

34For fungal ITS1 reads, open-reference OTU picking was performed using VSEARCH and the UNITE database. Taxonomy of representative fungal sequences was assigned using BLAST.

35Phyloseq was used to import BIOM data for 16S and ITS assays into R separately, followed by the determination of Shannon and Chao1 alpha diversity. For all analyses, bacterial data was subsampled for 2,500 reads. For fungal analyses subsampling was performed at 500 reads. Biomarker discovery analysis was performed using the LEfSe tool. Linear discriminant analysis (LDA) scores greater than 2.0 are considered to be significant. Microbial community state types (CSTs) were assigned on the basis of hierarchical clustering of the 20 most abundant OTUs. Prior to clustering, OTUs were agglomerated at the species level or the lowest identified taxonomic level. Clustering was performed using the wardD2 algorithm using Euclidian distances.

36R v3.4.2 was used for statistical analyses. The Kruskal-Wallis test was used to assess significance of continuous data. Linear regression was used to assess the significance of variables associated with the ordered outcome states of a HR-HPV infection (1). Logistic regression was performed using the GLM function and a binomial family generalized linear model in R. For categorical data, dummy variables were created and each individual factor level was tested in a univariate GLM analysis. Models were adjusted for age, smoking, HPV16 and CSTs. Smoking status was determined through a questionnaire and incorporated into the model as ordered categories: never smoked, former smoker and current smoker. Power of GLM results was computed using the lmSupport package.

37We performed a statistical mediation analysis to test whether V1 Gardnerella (an independent variable) could be acting by inducing a subsequent elevated microbiome diversity at V2 (mediator variable) that influences the outcome of HR-HPV progression using the package mediation. The outcome model we used was binary clearance/progression. Models were adjusted for age, CST, smoking status and HPV16 infection status. In the results we present the mediation effect (average causal mediation effects (ACME)), which is the total effect of V2 Shannon diversity and V1 Gardnerella minus the direct effect of V1 Gardnerella. Additionally, we estimate the direct of effect (presented using the average direct effect (ADE)) of V1 Gardnerella on the binary outcome clearance/progression, minus the effect of the V2 Shannon diversity mediator; the total effect, which is the sum between the indirect effect of the V2 Shannon diversity and the direct effect of the V1 Gardnerella; and the proportion mediated which is the ratio of the ACME and total effect estimates.

38Where authors are identified as personnel of the International Agency for Research on Cancer/ World Health Organization, the authors alone are responsible for the views expressed in this article and they do not necessarily represent the decisions, policy or views of the International Agency for Research on Cancer / World Health Organization.

39Proyecto Epidemiológico Guanacaste, Fundación INCIENSA, San José, Costa Rica—Bernal Cortés (specimen and repository manager), Paula González (LTFU: co-principal investigator), Rolando Herrero (CVT: co-principal investigator), Silvia E. Jiménez (trial coordinator), Carolina Porras (co-investigator), Ana Cecilia Rodríguez (co-investigator).

40United States National Cancer Institute, Bethesda, MD, USA—Allan Hildesheim (co-principal investigator & NCI co-project officer), Aimée R. Kreimer (LTFU: co-principal investigator & NCI co-project officer), Douglas R. Lowy (HPV virologist), Mark Schiffman (CVT: medical monitor & NCI co-project officer), John T. Schiller (HPV virologist), Mark Sherman (CVT: QC pathologist), Sholom Wacholder (statistician).

41Leidos Biomedical Research, Inc., Frederick National Laboratory for Cancer Research, Frederick, MD, USA (HPV Immunology Laboratory)—Ligia A. Pinto, Troy J. Kemp

42Georgetown University, Washington, DC, USA—Mary K. Sidawy (CVT: histopathologist)

43DDL Diagnostic Laboratory, Netherlands (HPV DNA Testing)—Wim Quint, Leen-Jan van Doorn, Linda Struijk.

44University of California, San Francisco, CA, USA—Joel M. Palefsky (expert on anal HPV infection and disease diagnosis and management), Teresa M. Darragh (pathologist and clinical management)

45University of Virginia, Charlottesville, VA, USA—Mark H. Stoler (QC pathologist)