Abstract
Objective
This study aimed to investigate the relationship between cervical intraepithelial neoplasia associated with human papillomavirus (HPV) infection and the vaginal microbiome.
Materials and Methods
In this study, the vaginal microbiota profile was compared among three groups of women: those with HPV infection and no cervical intraepithelial neoplasia (NILM, n=35), those with low-grade squamous intraepithelial lesions (LSIL, n=28), and those with high-grade squamous intraepithelial lesions (HSIL, n=24). Vaginal bacterial diversity was analyzed by deep sequencing of the V3-4 region of the barcoded 16S rRNA gene using the Illumina MiSeq platform, considering alpha diversity, beta diversity, and taxon classifications. Statistical significance was set at p<0.05.
Results
In the analyses performed using Chao1, Inverse Simpson, Shannon, and Observed indices, statistically significant differences were found among the groups in terms of all indices (p<0.05). Among groups, beta diversity did not show any notable differences. According to the “Linear Discriminant Analysis Effect Size” analysis, the taxa enriched in the HSIL group were Roseburia inulinivorans (p=0.0308), Micromonosporaceae family (p=0.0331), and Pirellula genus and species, (Planctomycetes), (p=0.0165); the taxa enriched in the LSIL group were Actinobaculum genus and species (p=0.0183). Lactobacillus helveticus and Faecalibacterium prausnitzii were more abundant in the NILM group, while Prevotella copri, Akkermansia muciniphila, and Fusobacterium species were more abundant in the LSIL and HSIL groups.
Conclusion
Our findings indicate that variations in the severity of cervical lesions are associated with notable alterations in vaginal microbiota composition. Further research is required to conclude which contribute to the formation of the cervical lesion and which are a consequence, among those that cause changes in the vaginal microbiota, of the lesion.
PRECIS: Comparison of human papillomavirus genotype distribution and microbiota profiles in women with varying grades of cervical epithelial abnormalities.
Introduction
Human papilloma virus (HPV) is the most common infection among sexually transmitted infections, and the lifetime infection rate in women who continue their sexual life can reach 70% and above(1). Epidemiological data show that most of these infections disappear spontaneously in the short term thanks to the body’s own immune mechanisms. However, in 10% to 20% of cases, the infection becomes persistent, and a long-term clinical course is observed(1, 2). Persistent infections predispose the development of cervical squamous intraepithelial lesions (SIL), which can develop into low-grade squamous intraepithelial lesions (LSIL) or high-grade intraepithelial lesions (HSIL). HSIL carries a high risk of carcinogenesis and requires early intervention against the risk of viral DNA integration into the cell genome. In contrast, LSIL cases are usually associated with transient HPV infections, and approximately 60% of cases tend to resolve naturally within one year(3).
The vaginal microbiome has long been identified as playing a critical role in maintaining vaginal health. Next-generation sequencing technologies provide more detailed and in-depth information about the structure and functions of bacterial communities in the vagina(4, 5). Dysbiosis in the vagina causes damage at the cellular and DNA level by producing bacterial genotoxins, thus triggering the formation of a chronic inflammatory environment favorable for cancer development(6). Chronic inflammation triggers the production of nitrosamines, paving the way for metabolomic profiles considered to be associated with cancer to play a decisive role in the interaction between cervical inflammation, HPV infection, and the local microbiome(7, 8).
Previous studies have revealed that Lactobacillus dominance decreased and bacterial diversity increased in the vaginal microbiome of women diagnosed with HPV persistence, LSIL, and HSIL(9, 10). In addition, although specific bacterial species associated with HPV and SIL have been identified, no consensus has been reached on this issue(11). Therefore, in this study, only the differences between vaginal dysbiosis and the normal microbiome of women diagnosed with LSIL and HSIL were analysed.
Materials and Methods
Ethical Approval
This is a prospective multi-center study of adult women admitted to the colposcopy units of Bozok University Research and Application Hospital and Etlik Zübeyde Hanım Gynecology Training and Research Hospital. Participants were enrolled between January 2023 and January 2024. The research protocol received ethical clearance from the Clinical and Interventional Research Ethics Committee of Yozgat Bozok University Research and Application Hospital, with official approvals granted on 25 August 2022 and 22 September 2022 (decision no: 2017-KAEK-189_2022.08.25_10, date: 28.08.2022 and decision no: 2017-KAEK-189_2022.09.22_06, date: 22.09.2022). All authors fully complied with the ethical principles of the Declaration of Helsinki.
Study Design and Sample Collection
Pre-menopausal, non-pregnant women over the age of 30 years who presented to the Clinic of Colposcopy, Yozgat Bozok University Research and Application Hospital and Etlik Zübeyde Hanım Training and Research Hospital, with a diagnosis of HPV positivity, were included in the study. In this study, we excluded individuals who met the following criteria: 1) those who had received antibiotic treatment in the last 7 days, 2) those who had used suppositories, vaginal medication or vaginal cleaning products in the 48 hours before the visit, 3) those who had a history of coitus within the last 48 hours, 4) those currently diagnosed with a sexually transmitted infection or vaginal infection, 5) those with diseases affecting the immune system, 6) those who had undergone cervical treatment in the past, 7) those using hormone replacement therapy or oral contraceptives, and 8) active smokers. All participants provided written informed consent.
A total of 87 women who participated in the study were divided into three groups according to their cervical pathological findings based on the Bethesda 2001 system(12) from the histopathological evaluation of the biopsy samples taken under colposcopy: No intraepithelial lesion or malignancy (NILM) (n=35), LSIL (n=28), and HSIL (n=24). The presence of HPV was detected by Linear Array HPV Genotyping Assays developed by Roche (Indianapolis, IN) and the subjects were divided into three categories as follows: 1) HPV 16 and 18 positive, 2) Other high-risk subtypes positive, and 3) Low-risk subtypes positive (Table S1). Other high-risk HPV types include HPV26, 31, 33, 35, 39, 45, 51, 52, 53, 56, 58, 59, 66, 68, 69, 73, and 82. The low-risk HPV comprises HPV6, 11, 40, 42, 43, 44, 54, 61, and 70. Some demographic characteristics of all included patients were also recorded. It is shown in Table 1.
Two gynecologists from each of the two centers used eNATTM kits (606CS01L, Copan Group, Copan Italia), to collect vaginal swabs without lubrication after inserting the vaginal speculum. Each vaginal sample was collected using a different swab and collecting tube, and it was promptly stored at -80 °C until DNA extraction.
DNA Isolation
After the samples were collected, they were all sent to the Diagen lab in Ankara, Türkiye, on dry ice so that DNA could be extracted. As directed by the manufacturer, vaginal swab samples were examined using the Kurabo QuickGene DNA tissue kit S (Japan).
DNA Amplification and Sequencing of 16S RNA Gene Fragments
Total genomic DNA was extracted from vaginal samples with the CTAB/SDS technique. DNA concentration and purity were evaluated via electrophoresis on a 1% agarose gel, after which the samples were adjusted to a final concentration of 1 µg/µL using sterile water. The V3-V4 hypervariable portions of the bacterial 16S rRNA gene were subsequently amplified according to the designated polymerase chain reaction (PCR) conditions. In this procedure, the primers 515F: 5’-GTGCCAGCMGCCGCGCGGTAA-3’ and 806R: 5’-GGACTACHVGGGTWTCTAAT-3’ were utilized. DNA amplification was conducted following the Illumina technique for 16S Metagenomic Sequencing Library Preparation (Part No. 15044223 Rev. A, Illumina, CA, USA). PCR processes commenced with an initial denaturation at 98 °C for 1 minute, succeeded by 30 cycles including denaturation at 98 °C for 10 seconds, annealing at 50 °C for 30 seconds, and elongation at 72 °C for 30 seconds. PCR products were assessed using 2% agarose gel electrophoresis and subsequently purified using the Qiagen Gel Extraction Kit (Qiagen, Germany). The PCR products were indexed using the Nextera XT Index Kit (Illumina), followed by the multiplexing step. Purified DNA libraries were constructed utilizing TruSeq DNA PCR-free sample preparation kits (Illumina, USA) and sequenced to acquire 250 bp paired-end reads using the Illumina NovaSeq platform. Following sequencing, barcodes were eliminated from the arrays, and the data were analyzed using the QIIME V.2.0 analytical pipeline. The amplicon sequence variations (ASVs) for each sample were limited to 10,000 reads through random selection(13, 14). The ASV taxonomy was established using the Silva 138 SSURef NR99 16S rRNA gene reference database(15). Reliability was maintained throughout the experimental process by employing negative controls in DNA extraction, PCR, and sequencing. The precision of the studies was enhanced by employing an internal standard, which included an internal control for 16S and 1007 reads, during sequencing.
Bioinformatics Analysis of Sequences
First, a quality evaluation was conducted using the Prinseq-lite (v0.20.4) software (http://prinseq.sourceforge.net) to begin bioinformatic analysis of the raw read data. The quality control step was conducted using the following parameters: 50 for min length, 30 for trim qual right, 20 for trim qual window, and mean for trim qual type. The FLASH software was used with the default parameters to merge the quality-controlled forward and reverse read sequences. After merging, human-derived sequences were identified by aligning individual reads to the human genome (GRCh38.p11, December 2013) using Bowtie2 software (very-sensitive mode: -D 20 -R 3 -N 0 -L 20 -i S,1,0.50) and subsequently removed.
Using the Bayesian-based ribosomal database project (RDP) Classifier (v2.12) algorithm, reads that did not match the human genome were taxonomically categorized after they were aligned to the RDP. A contingency table with taxonomic counts was created once each sample’s taxonomic distribution was established. QIIME (Quantitative Insights Into Microbial Ecology) v1.9.0 software was used to convert this table to Biom format, and the ecological diversity, abundance, and composition of microorganisms were analyzed(16).
The diversity within samples was measured using the alpha diversity assessment. With 9,500 random readings, 1,000 dilutions were made for every sample. The Shannon diversity index, Chao1 index, observed index, and inverse Simpson index were used per sample to determine diversity. The statistical software R (v3.1.0) was used to create box plots to visualize the results of the alpha diversity analysis.
The microbiological diversity of the samples was examined using beta diversity assessments. To do this, Bray-Curtis distance index matrices were used in principal coordinate analysis. Utilizing Adonis nonparametric analysis of variance, statistical differences across groups were assessed. Additionally, the non-parametric Wilcoxon test was used in R scripts to analyze groups and sample types.
Statistical Analysis
Data analysis was performed using IBM SPSS Statistics for Windows (version 27.0; Armonk, NY, USA). For continuous variables, descriptive statistics such as mean, standard deviation, and median were calculated. Categorical variables were presented as frequencies and percentages. The Student’s t-test was applied to compare numerical data between two groups. For comparisons involving more than two groups, either the Kruskal-Wallis H test or one-way ANOVA was utilized, depending on the normality of the data distribution. A p-value of less than 0.05 was considered statistically significant.
Results
There were no statistically significant differences among the groups in their demographic profiles, including age, height, weight, body mass index, gravida, parity, duration of marriage, or the number of prior abortions (p>0.05) (Table 1).
In the study, analyses of vaginal microbiota alpha diversity between NILM, LSIL, and HSIL groups using Chao1, Inverse Simpson, Shannon, and Observed indices, showed statistically significant differences between the groups in terms of all indices (p<0.05). The Chao1 index (Figure 1A) showed that the LSIL group had significantly higher species richness than the other groups (p<0.001). The Inverse Simpson index (Figure 1B), a measure that emphasizes the distribution of dominant species within the community, indicated higher diversity in the LSIL group; a statistically significant difference was found between the groups (p<0.001). The Shannon index (Figure 1C) was evaluated considering both species richness and even distribution. Alpha diversity analysis revealed significantly increased richness in the LSIL and HSIL groups relative to the NILM group (p=0.009). As shown in Figure 1D, the LSIL group exhibited the highest observed species count, whereas the NILM group demonstrated the lowest (p<0.001).
Bray-Curtis non-metric multidimensional scaling (NMDS) analyses (Figure 2) revealed that the groups did not differ significantly in terms of beta diversity. No significant difference was observed between the LSIL, HSIL, and NILM (normal) groups. Unweighted UniFrac, and weighted UniFrac NMDS analyses (Figure 2) also showed similar results. There was no significant difference in microbiota composition between the groups. Stress plot analyses showed that the NMDS model was reliable (R2>0.96) (Figure 3). The high model fit indicates that the analyses were technically correct, but there was no significant segregation between the groups (Figures 2 and 3).
According to linear discriminant analysis effect size (LEfSe) analysis used to determine the enriched microbiota species specific to the groups, the enriched taxa in the HSIL group were Roseburia inulinivorans (p=0.0308, LDA score=2.50), Micromonosporaceae family (p=0.0331, LDA score=2.07), Pirellula genus and species (Planctomycetes) (p=0.0165, LDA score=2.07); in the LSIL group, the enriched taxa were Actinobaculum genus and species (p=0.0183, LDA score=3.16); and in the NILM group, the enriched taxa were Deinocococci class (p=0.0301, LDA score=2.34), Thermus genus and Thermaceae family (p=0.0301, LDA score=2.34) (Figure 4 and S1).
In this study, the distribution of vaginal microbiota at phylum, class, order, family, genus, and species levels among NILM, LSIL, and HSIL groups was analyzed.
When analyzed at the phylum level, Firmicutes emerged as the predominant phylum across all groups. In contrast, the relative abundances of Bacteroidetes and Actinobacteria were reduced in both LSIL and HSIL groups. Meanwhile, phyla like Fusobacteria and Proteobacteria showed a marked increase, particularly in the HSIL group (Figure 5A).
At the class level, Bacilli and Clostridia were observed as dominant in all groups. Gammaproteobacteria and Actinobacteria had higher proportions in LSIL and HSIL groups; and especially Verrucomicrobiae and Synergistia increased in the lesion group (Figure S2). Lactobacillales, the most common bacterial group at the order level, had the highest proportion in all groups. Bacteroidales and Clostridiales were more prominent in the HSIL and LSIL groups, while Fusobacteriales and Campylobacterales were significantly increased, especially in the HSIL group (Figure S2). At the family level, Lactobacillaceae was the most common family in all groups. Prevotellaceae and Leptotrichiaceae were more prevalent in the HSIL group, while Bifidobacteriaceae and Ruminococcaceae were more prevalent in the LSIL and NILM groups. Fusobacteriaceae and Veillonellaceae were among the other important families showing an increase in the HSIL group (Figure S2). At the genus level, Lactobacillus was the dominant genus in all groups. Prevotella and Gardnerella were found at higher rates in both the LSIL and HSIL groups, while Atopobium and Fusobacterium appeared more frequently in the HSIL group. Streptococcus and Bifidobacterium were more common in the NILM group (Figure S2).
At the species level, Lactobacillus iners was the most dominant species in all groups. Lactobacillus helveticus and Faecalibacterium prausnitzii were more abundant in the NILM group, while Prevotella copri, Akkermansia muciniphila, and Fusobacterium species were more abundant in the LSIL and HSIL groups. Species such as Roseburia faecis and Dorea formicigenerans also showed higher abundance in the NILM group compared to the other groups (Figure 5B).
Discussion
Maintaining the balance of the cervicovaginal microbiota is essential for modulating the acquisition, persistence, and eventual elimination of HPV from the female reproductive tract. An imbalance in this microbial ecosystem, referred to as dysbiosis, can contribute to the initiation of pathological processes that may lead to carcinogenesis. These processes are mediated through disruptions such as impairment of epithelial integrity, metabolic dysfunction, irregular cellular proliferation, genomic instability, chronic inflammatory responses, and enhanced angiogenesis(17-19).
The vaginal microbiome is a constantly changing ecosystem that can be influenced by several factors(20). Ravel et al.(5) conducted a study where they categorized the vaginal microbiome of healthy females into four distinct community status types (CST), based on the dominating species of Lactobacillus and the corresponding pH values. The four CST types are CST-I, CST-II, CST-III, and CST-V, which correspond to the presence of Lactobacillus crispatus (26.2%), L. gasseri (6.3%), L. iners (34.1%), and L. jensenii (5.3%), respectively. In subsequent years, by assembling the largest dataset of human vaginal microbiota profiles, CST-IV was divided into seven subgroups according to the dominant non-lactobacillus species: CST IV-A; Candidatus lachnocurva vagina high rate G. vaginalis a medium rate, CST IV-B; G. vaginalis at a high rate, L. vaginae at a low rate, CST IV-C0; an equal community with moderate amounts of Prevotella, CST IV-C1; Streptococcus dominant, CST IV-C2; Enterococcus dominant, CST IV-C3; Bifidobacterium dominant, CST IV-C4; and Staphylococcus dominant(21). These members of the flora defend against harmful bacteria by generating hydrogen peroxide (H2O2) and bacteriocins, specifically generated by Lactobacillus species to maintain an acidic vaginal pH(22, 23). They produce hemolysin to enhance cytotoxicity (Prevotella)(24) and facilitate adhesion for epithelial colonization (Atopobium, Bifidobacteria, Gardnerella)(25).
CST III and IV have frequently been linked to the emergence of premalignant cervical lesions and the progression to invasive cervical cancer, particularly in the presence of HPV infection. Furthermore, Lactobacillus gasseri, the predominant strain of CST II, has been suggested to be associated with the fastest clearance of acute HPV infection in HPV-positive women(26, 27). Some studies report that HPV alters the composition of the vaginal microbiota by shifting it from CST III to CST IV, and that CST IV contributes to the progression and increased severity of cervical neoplasms(28). In the present study, Lactobacillus iners increased more in the LSIL group than in the HSIL and in the normal group. Lactobacillus gasseri was not detected in HSIL and LSIL groups. In accordance with the report by Lee et al.(29), which stated that Prevotella species are associated with HPV, Prevotella species were more prevalent in HSIL groups in our study. Fusobacteria were the predominant phylum, especially in the HSIL group. This is consistent with other studies emphasizing the association between Fusobacteria and cervical cancer so far(9, 29, 30).
LEfSe analysis conducted in our study revealed several candidate biomarkers that were linked to distinct stages of cervical pathology. Roseburia inulinivorans and members of the Micromonosporaceae family were found to be significantly more abundant in the HSIL group, while Actinobaculum was predominantly detected in LSIL samples. Roseburia inulinivorans is a motile bacterial species classified under the phylum Firmicutes, typically residing in the human colon, where it aids in butyrate synthesis through the fermentation of dietary polysaccharides(31). On the other hand, Actinobaculum is known to be a common pathogen in urogenital tract infections(32). These two agents are of interest for HSIL and LSIL, considering the exclusion criteria and meticulous sampling.
Akkermansia muciniphila was expressed at higher levels in HSIL and LSIL groups compared to other groups. Considering the reports in previous studies that Akkermansia muciniphila has the potential to prevent obesity, diabetes, and atherosclerosis and increase the efficacy of cancer immunotherapy, this dysbiosis member raises the question “Could it be a repair function for HSIL and LSIL?”(33, 34).
Our findings regarding alpha diversity mirrored the previously reported pattern of changes in species richness and community composition of the vaginal microbiota along the progression of cervical lesions(35). In contrast, the absence of a statistically significant difference in beta diversity across the groups could be attributed to the application of the NMDS model in our analysis. To assess the robustness of this model, we also performed stress plot evaluations. We also included stress plot analyses in our research to understand the reliability of this model.
Study Limitations
Given the relatively small sample size, there is a possibility that the statistical power may have been insufficient. Therefore, the results should be approached and interpreted cautiously. Among the patients included in the study, only specific HPV types -HPV6, 16, 18, 31, 35, 42, 45, 51, 52, 53, 54, 56, 58, 59, 61, 66, and 68- were detected (Table S1), while other HPV types were not observed. This restricted HPV spectrum may limit the generalizability of our results and could hinder a comprehensive understanding of the relationship between vaginal microbiota composition and cervical intraepithelial neoplasia. Future research with larger, more diverse cohorts, encompassing a broader range of HPV genotypes, is warranted to clarify whether the observed associations persist across different viral subtypes and to provide a more complete insight into this interaction.
Conclusion
The results of this study demonstrate that alterations in vaginal microbiota composition are associated with the severity of cervical lesions. However, further investigations are necessary to determine whether these microbial shifts contribute to the development of cervical lesions or arise as a consequence of the disease process.
