Association of cervical microbial community with persistence, clearance and negativity of Human Papillomavirus in Korean women: a longitudinal study

The present study aimed to identify the cervical microbes that are associated with HPV negativity, HPV clearance and HPV persistence and to assess the microbes’ longitudinal associations as related to HPV infection dynamics among Korean women. We enrolled 41 women with 107 samples, and classified them according to the HPV infection dynamics: HPV negativity (21 samples, 10 subjects), HPV clearance (42 samples, 15 subjects), and HPV persistence (44 samples, 16 subjects). Cervical swabs were collected at the baseline and six-month-interval follow-up visits. HPV positivity was determined by HPV DNA HC2 assay, and the microbiome was analyzed using 16SrRNA pyrosequencing, linear discriminant analysis effect size and multivariate logistic analysis. In the multivariate logistic analysis results, Lactobacillus crispatus (multivariate OR (mOR) = 8.25, 95% CI 2.13~32.0) was predominant in the HPV-negative group. We observed that Eubacterium eligens (mOR = 11.5, 95% CI 1.31~101.4), Gardnerella vaginalis (mOR = 17.0, 95% CI 2.18–131.8), and Ureaplasma urealyticum (mOR = 7.42, 95% CI 1.3–42.46) had the strongest associations with HPV clearance, and Lactobacillus johnsonii (mOR = 16.4, 95% CI 1.77–152.2) with HPV persistence. Overall, greater diversity was observed in HPV-persistence than in HPV-negative women. Our findings suggest that the presence and prevalence of a specific cervical microbiome are factors involved in HPV dynamics.

Identification of cervical microbiome markers of HPV dynamics. The taxonomic groups that were relatively abundant in the HPV persistence, negativity and clearance groups were identified by linear discriminant analysis effect size (LEfSe with α = 0.05), an LDA score of at least 2, and a relative abundance greater than 0.1, both at the baseline and for total visits. The changes at the baseline are depicted in a cladogram (Fig. 1A). In the HPV-persistence group, significant over-representation of Haemophilus (Order: Pasteurellales; Family: Pasteurellaceae) (P = 0.0345) (Fig. 1A,B,E) was observed. By contrast, the HPV-clearance group was found to have significantly higher levels of Gardnerella vaginalis (Order: Bifidobacterialies; Family: Bifidobacteriaceae) (P = 0.0403) (Fig. 1A,B,D), and the HPV-negative group was determined to have significant overrepresentation of Lactobacillus crispatus (P = 0.0067) (Fig. 1A-C).
In the HPV-clearance group, Eubacterium eligens, Ureaplasma urealyticum, and Gardnerella vaginalis were over-represented ( Fig. 2A,B). The HPV-clearance group also was found to have significantly higher levels of Gardnerella vaginalis (P = 0.0028), Eubacterium eligens (P = 0.0068) and Ureaplasma urealyticum (P = 0.0112) ( Fig. 2D-F). We evaluated species richness (Chao 1 index), alpha-diversity (Shannon index) and beta-diversity of the cervical microbiome in baseline ( Fig. 3A-C) and total visit samples ( Fig. 3D-F). The results showed species richness (Chao1) was greater in the HPV-clearance and persistence subjects than in the HPV-negative subjects (p < 0.01) in total visits. Further, beta diversity differed between HPV negative, clearance and persistence group in total visit samples (bray p < 0.005). The baseline visit not showed any significance.

Discussion
In the present study, we assessed the longitudinal associations of cervical microbes in HPV-negativity, clearance and persistence women. The main findings showed that the cervical microbiome differs significantly by HPV group. We found that women with high proportions of Lactobacillus johnsonii, Haemophilus (genus) and Mycoplasmataceae (family) had the strongest association with HPV persistence, and that women with a cervical microbiome dominated by Eubacterium eligens, Gardnerella vaginalis and Ureaplasma urealyticum had the strongest associations with HPV clearance. Women with a high proportion of Lactobacillus crispatus were most   Based on these results of our two-year longitudinal study, we can suggest that bacterial dysbiosis is a factor associated with HPV dynamics (progression and regression). This study is, to the best of our knowledge, the first to examine the relationship between Haemophilus and Lactobacillus johnsonii in HPV-persistent women.
Our multivariate logistic analysis confirmed that Lactobacillus johnsonii was significantly associated (mORs = 16.4, 95% CI = 1.77-152.2) with HPV persistence among our Korean subjects. There is no previous study on Lactobacillus johnsonii in relation to cervical microbiomes, though it has been reported in saliva samples of HPV-positive and HPV-negative oropharyngeal cancer patients. A recent study by Guerrero-Preston et al. reported bacterial species Lactobacillus gasseri/johnsonii and Lactobacillus vaginalis in saliva of HPV-positive and HPV-negative oropharyngeal cancer patients 15 . However, the interaction between these microbes in saliva with HPV persistence and cancer risk was not discussed. The LEfSe analysis in the present study also confirmed enrichment of the genus Haemophilus and family Mycoplasmataceae in HPV-persistence women. The results for Mycoplasmataceae are in agreement with the finding of Abedamowo et al., that Mycoplasma hominis in cervical microbiota is significantly associated with hr-HPV infection 16 . It has also been reported that a few Mycoplasma species are efficient methylators and can promote cervical carcinogenesis through methylation of hr-HPV and cervical somatic cells 16 . One of the previous cervical microbiota studies (16-weeks samples) conducted by Brotman et al. 12 demonstrated that high proportions of Atopobium spp. (belonging to CST IV-B) were associated with the slowest HPV clearance rate and that the persistence effect may be due to the high proportion of Atophobium vaginae that can play a role in the disruption of the epithelial barrier. The plausible molecular mechanisms of HPV persistence infection have been proposed: the viral life cycle and immune system evasion (avoidance of triggering of an immune response by the host) 17 , the tethering mechanism on host mitotic chromosomes (to ensure that the viral genome is not lost during cell division) 17 , and E2 interaction with host receptors 17 . Then, 16 S rRNA gene sequencing revealed an increased proportion of L. johnsonii in HPV persistence and depletion of Lactobacillus crispatus relative to HPV-negative women; additionally, an alpha-diversity analysis provided further confirmation that HPV persistence in the total-visit samples showed increased species richness. Therefore, it can be suggested that changes in cervical microbiota with depletion in Lactobacillus crispatus and increased microbial diversity promote HPV infection and might be involved in HPV persistence 18 . The alpha-diversity results (Chao 1 index) showed increased diversity in HPV-persistence women relative to HPV-negative women in the total-visit samples. Therefore, high bacterial diversity in HPV persistence can represent an environment that increases the risk of progressing HPV persistence.
In the HPV-clearance women, multivariate logistic analysis and LEfSe analysis confirmed the significant enrichment of species G. vaginalis, E. eligens and Ureaplasma urealyticum. The clearance effect might be due to the innate immune response (Toll-like receoptors-TLR) that is involved in a cascade of events promoting HPV clearance. TLR 3 and TLR9 play a role in regulating the pro-inflammatory cytokine and anti-viral environments of the lower female genital tract during viral and bacterial infection 19 ; intense inflammatory response induced by G. vaginalis in viral clearance 20 , as well as viral proteins 21 and other intrinsic host factors (such as sexual behaviors) 21,22 responsible for HPV clearance. Another plausible explanation is that this clearance effect is due to an intense inflammatory response or that higher concentrations of inflammation markers (IP-10 and MIG chemokines) induced by bacteria in the cervical region assist in viral clearances 23 . One study conducted by Moscicki et al. 20 reported that the HPV viral clearance effect might be due to the serendipity effect of Neisseria gonorrhe. One previous study 24 reported the lowest clearance rate in women whose cervical microbiota were dominated by Atopobium and Gardnerella; therefore, further studies are needed to identify the mechanism or mechanisms whereby the composition of the cervical microbiota influence clearance. Among the HPV-negative women, the Lactobacillus crispatus population remained stable at the baseline and follow-up visits, and Lactobacillus crispatus was significantly over-represented. Interestingly, in the LEfSe analysis results, we observed Corynebacterium sundsvallense, Facklamia hominis, Actinobaculum schaalii and Helcococcus ovis for the first time in the HPV-negative women. As reported earlier, in healthy women, hydrogen-peroxide producing Lactobacillus spp. are stable in maintaining the vaginal ecosystem 25 . Klebanoff et al. reported that Lactobacillus sp. increase TNF-α and IL-1α production, activate NF-κB in THP-1 cells and increase TNF-α production by human monocytes. This suggests that a higher concentration of Lactobacillus in the vagina can influence the physiology and host defenses therein 25 . In our previous study, we reported that Lactobacillus crispatus and Lactobacillus inners (Cluster III) were the dominant species in HPV-negative women and that high abundance of L. crispatus is associated with a low risk of CIN 26 .
The strength of our present research is that the HPV samplings were longitudinal and that the differences between the cervical microbial compositions (species and genus level) among the HPV-negativity, clearance and persistence samples were identified by 16 S rRNA pyrosequencing. One of the most notable findings is that the bacterial communities in the cervical region of the HPV-persistence and clearance women were distinctive from those in the HPV-negative women, suggesting the presence of certain selective pressures contributing to the shift in the cervical microbiota. However, we also must acknowledge certain study limitations: (1) small sample sizes among the study groups, which could have biased the results; (2) an advanced-age study population (mean age: 44). So, the results should be interpreted with caution, and especially, it should be recognized that they might not apply to the general population; (3) another minor limitation of this study is HPV grouping by the liquid hybridization assay Hybrid Capture 2 (HC2) assays. It has been noted that HC2 assay shows some false positive results of 5%. Despite these limitations, HC2 assay has several advantages over other commercially available PCR-based target amplification method [27][28][29] . Moreover, HC2 assay is technically well designed and can be easily controlled and performed by the laboratory technician. Whereas PCR-based method contains several steps that must be carefully optimized which make it more difficult to standardize PCR in laboratory condition.
In conclusion, this study analyzed various cervical microbial taxa among different HPV groups, and found that increased bacterial diversity with reduced L. crispatus species could be related to HPV persistence and also that the presence and prevalence of a specific cervical microbiome could be relevant to HPV dynamics.

Methods
Subject recruitment. We collected cervical samples from Normal and ASCUS patients (age between 18 and 65 years) who had participated in the Korean human papillomavirus (HPV) cohort study from 2006 to 2013. Detailed information on the inclusion/exclusion criteria for the HPV cohort is provided in our previous paper 30 . The selected participants were interviewed using a structured questionnaire on their socio-demographic characteristics (Table 1), after which they underwent physical and gynecological examinations. Cervical samples were collected using a Cervix brush (Rovers Medical Devices, Oss, The Netherlands), after which the brush was immediately soaked in a vial of PreservCyt solution (Cytyc Corporation, Marlborough, MA, USA) fixed within a Thin Prep processor. The cytological findings were grouped based on the Bethesda system 31 . The samples were examined for the presence of one or more hrHPV types using the Hybrid Capture 2 (HC2) assay (16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, and 68) according to the manufacturer's instructions. The values obtained were recorded in relative light units (RLUs). All relative light units (RLU) measured on a luminometer were divided by the RLU of the respective positive control to provide a ratio. The sample was classified as a positive when the RLU/ PC ratio was 1 pg/mL or greater. After baseline sampling, follow-up visits were scheduled at six-month intervals. The collected samples were stored at −80 °C for further analysis.

Grouping of HPV status for 2-year follow-up.
To investigate the shifts of bacterial communities in the HPV states, we grouped the samples based on the following conditions: (1) HPV negative -persistently negative detection of hrHPV types throughout 24 months (all negative during observation period, 21 samples, 10 subjects); (2) HPV clearance -women who were initially hrHPV (baseline), and then regressed to negative during follow-up (42 samples, 15 subjects); (3) HPV persistence -persistently positive detection of hrHPV types at the baseline and during follow-up (all positive during observation period, 44 samples, 16 subjects) (Fig. 4). Of the 41 subjects, 35 participated at the 2 nd visit, 20 at the 3 rd visit, 10 at the 4 th visit, and 1 at the final visit.
High-risk HPV DNA detection. HPV DNA was detected using the Digene HC2 high-risk DNA test (Qiagen, Gaithersburg, MD, USA) with signal amplification and chemiluminescence for 13 types of HR-HPV scored in RLU/PC. A positive result indicated a concentration of 1 pg/ml or higher than the RLU/cutoff ratio (RLU of the specimen/mean RLU of 2 positive controls).
DNA isolation and pyrosequencing. DNA from cervical samples were isolated using the Fast DNA SPIN kit (MP Biomedicals, Santa Ana, CA, USA). The 16 S universal primers 27 F (5′ GAGTTTGATCMTGGCTCAG 3′) and 518 R (5′ WTTACCGCGGCTGC-TGG 3′) were used for amplification. A 20 ng aliquot of each sample was used for a 50 µl PCR reaction containing 10 × taq buffer, a dNTP mixture (Takara, shiga, Japan), 10 μm of the bar-coded fusion primers, and 2 U of taq Polymerase (extaq, takara). The PCR conditions and pyrosequencing protocols are available elsewhere 32  Next-Generation Sequencing using 454 GS-FLX plus. The raw sequences for the samples were arranged using a unique barcode, and low-quality reads (average quality score <25 or read length <300 bp) were removed 32 . The primer sequences were cut down by employing pairwise sequence alignment, and sequences were gathered to correct for sequencing errors. Taxonomic identification was performed using the EzTaxon-e public database 33 according to the highest pairwise similarity among the BLASTN search results. Possible chimera sequences were removed by the UCHIME algorithm 34 , and the diversity indices were calculated in Mothur after normalization of the read number in each sample. The potential biomarkers linked to HPV negative, HPV clearance and HPV persistence were analyzed by LDA effect size (LEfSe). Finally, the effect relevance was predicted by LDA 35 . Statistical analysis. The Chi-square test was used for testing of the categorical variables, and the Kruskal Wallis test was employed for comparison of the categorical continuous variables. When the number of expected frequencies was less than 5 and the number of cells was more than 25%, Fisher's exact test was performed. To find a significant difference in alpha diversity, the Kruskal Wallis rank sum test was used to evaluate the differences in diversity among the three groups, followed by Dunn's test of multiple comparisons. Multivariate logistic analysis was performed after adjusting for age, menopausal status, oral contraceptive use, and smoking habit. Risk estimates are presented as OR with 95% CI. LEfSe analysis was conducted to find significant differences among the relative abundance taxa 35 . Beta diversity was calculated with principal coordinates analysis (PCoA) according to the Bray-Curtis distances. A permutational multivariate analysis of variance (PERMANOVA) was implemented to determine significance in distance. Diversity and the PERMANOVA results were analyzed using the R packages "vegan" 36 . The statistical analysis was performed with SAS 9.4, and R version 3.3.1 with ggplot2 packages was used for visualization 37 .
Ethical Statement. The Ethics Committee of the National Cancer Center approved this study, and all of the experiments were performed in accordance with the approved guidelines and regulations (IRB No. NCC2016-0147). Written informed consent was obtained from all of the study participants in accordance with good clinical practice.

Data Availability
The datasets generated during the current study are available from corresponding author on reasonable request.