Research on the relationship between HPV infection and alterations in vaginal microbial ecology

Abstract

Objective

To investigate the changes in vaginal microbiota under different HPV infection statuses and explore the correlation between vaginal microbiota alterations and HPV infection.

Methods

151 cervical samples from gynecological outpatients were grouped into HPV-negative (HN, N = 51), transient infection (HTI, N = 42), and persistent infection (HPI, N = 58). Vaginal secretions were collected to assess microecology (pH, vaginal cleanliness, hydrogen peroxide, leukocyte esterase) via genital secretion analyzer. 16S ribosomal RNA (rRNA) sequencing analyzed vaginal microbiota characteristics, community state types (CST), richness, diversity, and biomarkers.

Results

16S rRNA sequencing identified 5 CST II, 52 CST III, and 94 mixed CST IV samples, showing diverse microbiota. Compared with HN, HTI and HPI had lower vaginal cleanliness, higher sialidase activity, elevated pH, and fewer Lactobacilli (P < 0.05). Lactobacillus iners dominated all groups, while Sneathia amnii was significantly higher in HPI (P < 0.05). HPV infection increased vaginal microbiota richness (HPI > HTI/HN, P < 0.05), with distinct group compositions (P < 0.05). Linear Discriminant Analysis Effect Size identified Lactobacillus gasseri, Atopobium vaginae, and Lactobacillus jensenii as biomarkers.

Conclusion

This study found significant differences in microbial community characteristics under different HPV infection statuses. The identification of biomarkers in vaginal microbiota under different infection statuses could provide new targets for clinical screening and prevention of cervical cancer.

Resumen

Objetivo

Investigar los cambios en la microbiota vaginal en diferentes estados de infección por VPH y explorar la correlación entre las alteraciones de la microbiota vaginal y la infección por VPH.

Métodos

Se agruparon 151 muestras cervicales de pacientes ambulatorias ginecológicas en grupos de VPH negativo (VN, N = 51), infección transitoria (IVT, N = 42) e infección persistente (IVP, N = 58). Se recolectaron secreciones vaginales para evaluar la microecología (pH, limpieza vaginal, peróxido de hidrógeno, esterasa leucocitaria) mediante un analizador de secreciones genitales. Se realizó secuenciación del ARN ribosomal (rARN) 16S para analizar las características de la microbiota vaginal, los tipos de estado comunitario (TEC), la riqueza, la diversidad y los biomarcadores.

Resultados

La secuenciación del 16S rRNA identificó 5 muestras de TEC II, 52 de TEC III y 94 de TEC IV mixto, lo que mostró una microbiota diversa. En comparación con el grupo VN, los grupos IVT y IVP presentaron menor limpieza vaginal, mayor actividad de sialidasa, pH elevado y menos Lactobacillus (p < 0,05). Lactobacillus iners dominó en todos los grupos, mientras que Sneathia amnii fue significativamente más abundante en el grupo IVP (p < 0,05). La infección por VPH aumentó la riqueza de la microbiota vaginal (IVP > IVT/VN, p < 0,05), con composiciones grupales distintas (p < 0,05). El análisis de tamaño del efecto por análisis discriminante lineal identificó a Lactobacillus gasseri, Atopobium vaginae y Lactobacillus jensenii como biomarcadores.

Conclusión

Este estudio encontró diferencias significativas en las características de la comunidad microbiana en diferentes estados de infección por VPH. La identificación de biomarcadores en la microbiota vaginal en diferentes estados de infección podría proporcionar nuevos objetivos para el cribado clínico y la prevención del cáncer cervical.

Introduction

Cervical cancer is a malignant tumor originating in the cervix, with the human papillomavirus (HPV) serving as one of its primary causative agents. HPV belongs to the Papillomaviridae family and Papillomavirus genus,¹ and it is a double-stranded DNA (dsDNA) virus primarily transmitted through sexual activity. HPV strains involved in genital tract infections are categorized into high-risk HPV (HR-HPV) and low-risk HPV (LR-HPV) types. HR-HPV infection is internationally recognized as a carcinogenic factor for cervical cancer.² Studies have indicated that the overall positivity rate for HR-HPV among the Chinese population is 22.7%,³ with the primary infection types being HPV52, 58, 16, 51, and 66. Among these, HPV16 is the primary carcinogenic type, followed by HPV52 and HPV18.⁴ Vaccination against HPV16, HPV18, and other high-risk types can prevent up to 90% of HPV-related infections.^5,6 On the other hand, LR-HPV is associated with mild squamous epithelial lesions and urogenital warts.

Studies have shown that HPV-positive women exhibit an increased diversity in their vaginal microbiota, accompanied by a decrease in the relative abundance of Lactobacillus spp. and an elevation in pH levels.^7,8 The vaginal microbiota constitutes a vital component of the female vaginal microenvironment and can be classified into five community state types (CST) based on differences in species composition and their relative abundances. CST-I, CST-II, and CST-III are predominantly characterized by Lactobacillus crispatus, Lactobacillus gasseri, and Lactobacillus iners, respectively, while CST-IV-A and CST-IV-B are marked by a reduced proportion of Lactobacilli and an increased proportion of anaerobes.⁹ A healthy vaginal environment is typically dominated by Lactobacillus spp., which maintains an acidic vaginal milieu in women¹⁰ and serves as the first line of defense against sexually transmitted infections.¹¹ HPV disrupts the original balance of the vaginal microbiota by inhibiting the growth of Lactobacilli and increasing species diversity, particularly through enzymes and metabolites produced by strict anaerobes.^12,13 This disruption compromises the protective barrier of cervical epithelium, making it easier for HPV to integrate into epithelial cells and promote cervical intraepithelial neoplasia (CIN) lesions.^14,15 Therefore, vaginal microbiota can be considered as a biomarker for screening HPV infection or precancerous lesions of the cervix.

Bacterial vaginosis (BV) arises from the disruption of the vaginal microenvironment and microbial community, primarily associated with the reduction of Lactobacillus species and the increase of anaerobes such as Gardnerella, Atopobium, and Prevotella.^16–18 The microbial community structure of BV corresponds to CST-IV. This study focused on investigating the differences in microbial community characteristics under various infection states by combining clinical samples with 16S rRNA high-throughput sequencing technology. In particular, it explored the evolutionary patterns of microbial communities from transient to persistent infections, laying a crucial foundation for understanding the microenvironmental regulatory mechanisms of HPV infection. These findings not only deepen our comprehension of HPV-associated vaginal microbial dysbiosis but also provide a scientific basis for the development of an early risk warning system based on microbial biomarkers. The proposed microecological regulatory strategies, especially the novel approach of intervening in the HPV infection process by maintaining the dominance of lactobacilli, open up potential avenues for future clinical translational research aimed at preventing persistent HPV infections and cervical lesions. This research holds significant theoretical value and clinical application potential.

Materials and methods

Study subjects

Building on our research group's prior epidemiological findings indicating significantly elevated HPV infection rates in Xinjiang compared to other Chinese regions, this study aimed to further elucidate the epidemiological characteristics and clinical pathogenesis patterns through prospective case collection. Between June 2021 and June 2022, we enrolled 151 eligible female patients at the Gynecology Clinic of the People's Hospital of Xinjiang Uygur Autonomous Region. All clinical data were systematically collected using standardized protocols to ensure data quality and comparability for subsequent analysis of regional HPV infection profiles. HPV infection in cervical exfoliated cell samples was detected using an HPV nucleic acid typing detection kit (flow-through hybridization method, covering 37 genotypes) in combination with the Hybribio HPV-DNA Chip Typing Flow Hybridization Technique. The principle is based on the chemical color development of hybridization signals resulting from the binding of probes to the viral genome. Specifically, if no color appears in the spots, it indicates the absence of HPV infection. Conversely, if the spots exhibit a distinct color, it signifies the presence of HPV infection. The number and color intensity of the spots are used to determine multiple infections. HPV testing was conducted again one year after the initial infection. If no HPV infection was detected at this follow-up, it was classified as a temporary infection; if the infection persisted, it was categorized as persistent infection. Participants were divided into three groups: the HPV negative (HN, N = 51), the HPV transient infection (HTI, N = 42), and the HPV persistent infection (HPI, N = 58). As this study involved biological samples, ethical principles outlined in the Helsinki Declaration were strictly adhered to. The study was approved by the Ethics Committee of People's Hospital of Xinjiang Uygur Autonomous Region, and written informed consent was obtained from patients before single-use sterile flocked cervical sampling swabs were used to collect vaginal secretions from the middle to upper sections of the vaginal wall.

Inclusion and exclusion criteria

Female participants meeting the following criteria are eligible for inclusion in this study: women participants who have abstained of sexual intercuourse for three days and have no prior history of hysterectomy. Additionally, All the following requirements must have been met within the past six months: Participants should not have had cervical or vaginal surgery; they should have no history of cervical lesions or cervical cancer; they must be free from autoimmune diseases; and they should not be taking oral immunosuppressants or antibiotics. Participants must also not be pregnant at the time of sampling. All initially HPV-positive cases demonstrated monotypic infections (only single HPV genotype) at baseline, with systematic follow-up confirming either complete viral clearance or persistent carriage of the identical HPV genotype at the one-year follow-up assessment. All HPV patient samples selected for this study were high-risk samples (one of the 14 high-risk genetypes). In contrast, the negative samples were those that tested negative for all 37 of the aforementioned HPV genotypes. Females who fulfill these criteria will be included in the study, while those who do not will be excluded.

Vaginal microecological evaluation

In accordance with the ‘Expert Consensus on the Clinical Application of Vaginal Microbiological Evaluation’ established by the Chinese Obstetrics and Gynaecology Infection Collaboration Group, vaginal swab samples are analyzed using a reproductive tract secretions comprehensive analyzer supplied by Jiangsu Shuoshi Biotechnology Co., Ltd. This analyzer utilizes dry chemistry enzyme-based and latex immunochromatography-based detection kits for in vitro qualitative testing of female vaginal secretions. The analysis evaluates indicators such as vaginal (or leukorrhea) cleanliness (determined by the number of Lactobacilli and epithelial cells), hydrogen peroxide levels, leukocyte esterase activity, sialidase activity, pH value, and assesses the microbial ecological balance.

16S rRNA sequencing

The gathered samples were dispatched to Novogene for examination. Initially, genomic DNA was extracted using the magnetic bead method, followed by PCR amplification targeting the 16Sv34 region with the primer sequences: CCTAYGGGRBGCASCAG and GGACTACNNGGGTATCTAAT. The PCR amplification protocol is as follows: “Hot lid set at 105 °C; Initial denaturation at 98 °C for 1 min; followed by 14 cycles of [denaturation at 98 °C for 20 s, annealing at 60 °C for 15 s, and extension at 72 °C for 15 s]; a final extension at 72 °C for 2 min; and then hold at 4 °C indefinitely.” The PCR products were then subjected to electrophoresis using a 2% agarose gel for detection. Qualified PCR products were purified using magnetic beads and quantified using an enzyme-linked immunosorbent assay (ELISA). After mixing equal amounts of PCR products based on their concentrations and ensuring thorough mixing, the PCR products were again subjected to electrophoresis using a 2% agarose gel. The target bands were recovered using a universal DNA purification and recovery kit (Tiangen, Catalog No.: DP214). The NEB Next® Ultra™ II FS DNA PCR-free Library Prep Kit (NEB/E7430L) was utilized for library preparation. The constructed libraries were quantified using Qubit and Q-PCR. Upon qualification, the libraries were sequenced on a NovaSeq 6000 platform with PE 250 sequencing for the detection and analysis of vaginal microbial flora species and density.

Bioinformatics analysis

Extracting individual sample data from sequencing output based on Barcode sequences and PCR primer sequences. Subsequent to trimming the Barcode and primer sequences, FLASH¹⁹ (Version 1.2.11) is utilized to concatenate the reads of each sample, yielding raw tags. These raw tags are then rigorously filtered using the fastp software (Version 0.23.1).²⁰ The raw tags undergo rigorous filtering using the fastp software (Version 0.23.1). The filtering criteria are as follows: a minimum quality score (Q20) of 90 or above, a minimum sequence length retention of 200 base pairs (bp), and a GC content maintained within the range of 40–60%. Through this stringent filtering process, high-quality Clean Tags are obtained. Following this process, the tags undergo chimera sequence removal. This involves aligning the tag sequences with species annotation databases (Silva database for 16S/18S at https://www.arb-silva.de/) to detect and subsequently eliminate chimera sequences, resulting in the final Effective Tags.²¹ The effective data is further denoised using DADA2 to yield the ultimate Amplicon Sequence Variants (ASVs).

Alpha Diversity was utilized to analyze the microbial community diversity within each sample. The diversity and evenness of pathogenic microorganisms were evaluated using indices such as Shannon and Simpson. Subsequently, Beta Diversity was employed to compare and analyze the microbial community compositions across different samples, with the Unweighted Unifrac distance metric chosen to measure the dissimilarity coefficients between samples. Non-metric Multidimensional Scaling (NMDS) analysis was then conducted to visualize the inter-group and intra-group differences among samples. Finally, the MetaStat and LEfSe methods were applied to analyze the species abundance data between groups, aiming to identify species with significant differences and thereby discover statistically distinct Biomarkers between the groups.

Sequencing annotation and diversity analysis

For the obtained ASVs, the classify-sklearn algorithm in QIIME2²² was employed to annotate each ASV using a pre-trained Naive Bayes classifier, yielding corresponding species information and abundance distribution based on species. Simultaneously, abundance and Alpha diversity calculations were conducted for the ASVs, with Chao1, Shannon, and Simpson indices used to evaluate species richness and evenness within samples. On the other hand, multiple sequence alignments were performed on the ASVs to construct a phylogenetic tree. Beta diversity was explored to investigate differences in community structure among different samples or groups through dimensionality reduction analyses such as NMDS, as well as the presentation of sample clustering trees. The Kruskal–Wallis rank-sum test and LEfSe statistical analysis methods were selected to test for significant differences in species composition and community structure among grouped samples.

Statistical analysis

All data were processed using GraphPad Prism 8.0. Comparisons between groups were conducted using the chi-square test. A P-value < 0.05 indicated a statistically significant difference, while a P-value < 0.01 indicated a highly statistically significant difference. Additionally, this study utilized the R package treeio for reading phylogenetic trees, ggtree for visualization, phyloseq for integrating microbiome data, and ggplot2 for labeling differential groups. The analysis was conducted using R version 4.3.0 to ensure reproducibility.

Results

Patient demographics and outcomes of vaginal microecological assessments

To investigate the variations in microbial communities across diverse HPV infection statuses, this study categorized samples based on specific criteria: an HN cohort comprising 51 cases, where HPV tests remained negative for a consecutive 12 months; an HTI cohort of 42 cases, where patients tested positive for HPV at enrollment but subsequently tested negative within the 12-month follow-up period; and an HPI cohort encompassing 58 cases, where patients maintained an HPV-positive status throughout the follow-up. In total, 151 patients were enrolled in this study. The mean ages of HN, HTI, and HPI cohorts were (41.72 ± 8.99) years, (39.90 ± 10.05) years, and (42.79 ± 11.83) years, respectively, with no statistically significant age differences observed among the three cohorts (P > 0.05). This suggests that age may not be a confounding factor in the analysis of microbial community variations across different HPV infection statuses.

The results of the patients’ vaginal microecological assessments revealed that (Table 1), in comparison to HN, both HTI and HPI cohorts exhibited a decline in vaginal cleanliness, an increase in sialidase positivity, an elevation in pH levels, and a decrease in Lactobacillus counts. These disparities were statistically significant (P < 0.05). When comparing HTI with HPI, there were no significant differences observed in any of the indicators (P > 0.05) except for Hydrogen peroxide (P < 0.05).

Table 1.

Comparison of biochemical indices of vaginal microecology under different HPV infection states.

Indices	HN (51 cases)	HTI (42 cases)		HPI (58 cases)
		Number	HTI vs HNP value	Number	HPI vs HNP value	HPI vs HTIP value
Vaginal cleanliness (III–IV)	15 (29.4%)	30 (71.4%)	<0.001	42 (72.4%)	<0.001	>0.99
Hydrogen peroxide (+)	37 (72.6%)	34 (80.9%)	0.46	33 (56.9%)	0.11	0.017
Sialidase (+)	16 (31.7%)	28 (66.7%)	0.002	38 (65.5%)	0.0005	>0.99
Leukocyte esterase (+)	28 (54.9%)	28 (66.7%)	0.29	32 (55.2%)	0.34	0.30
pH (>4.5)	15 (29.4%)	28 (66.7%)	0.0004	33 (56.9%)	0.007	0.41
Lactobacillus (few)	20 (39.2%)	30 (71.4%)	0.003	40 (69.0%)	0.002	0.83

Note:

(1) Vaginal cleanliness: Grades I and II fall within the normal range (negative). Grade III indicates moderate inflammation, while Grade IV indicates severe inflammation. Grades III–IV are thus defined as positive.

(2) There is a negative correlation between vaginal cleanliness and the quantity of Lactobacilli. The more severe the vaginal cleanliness issue (i.e., the higher the grade), the fewer the Lactobacilli present. Grades III–IV are defined as positive in terms of cleanliness.

(3) Hydrogen Peroxide (H₂O₂): An abnormal red color indicating a value ≤2 μmol/L is considered positive.

(4) Sialidase: A concentration > 9 U/mL, resulting in a red or purple color, is positive.

(5) Leukocyte Esterase: A value > 7 U/mL, showing a light blue color, is positive.

(6) Lactobacilli play a crucial role in maintaining the normal acidic environment (pH 3.8–4.5) of the vagina. When the pH > 4.5, it indicates that the vaginal microbial balance has been disrupted.

16S rRNA genomic sequencing and microbial community characteristics

After conducting 16S rRNA genomic sequencing on samples from 151 participants, we obtained a total of 21,537,348 raw reads. Following rigorous quality control, which involved the removal of adapters, barcodes, primers, low-quality bases, short reads, and chimeras, we were left with 18,438,201 valid reads, representing an average efficiency rate of 85.61%. Utilizing the DADA2 algorithm for denoising and dereplication, we identified a total of 8352 unique operational taxonomic units (OTUs), with an average of 55.31 OTUs per sample. Within these OTUs, HN harbored 2814 OTUs (an average of 55.18 per sample), HTI had 3076 OTUs (an average of 73.24 per sample), and HPI contained 4815 OTUs (an average of 83.02 per sample). Notably, 528 OTUs were shared among all three groups, accounting for 6.32% of the total OTUs detected (Fig. 1A).

Fig. 1.

Differences in microbial communities between samples and groups. Note: A is a Venn diagram showing the number of OTUs detected for characteristic sequences in samples grouped by different HPV statuses; B is a bar chart displaying the top ten types of pathogenic microorganisms detected in samples grouped by different HPV statuses and their respective proportions; C is a heatmap illustrating the classification of bacterial CSTs and the relationship between CSTs and sample groupings.

At the species level, the dominant taxa included L. iners, Gardnerella vaginalis, Prevotella melaninogenica, Streptococcus anginosus, L. gasseri, Sneathia amnii, Atopobium vaginae, Prevotella timonensis, Lactobacillus delbrueckii, and Prevotella bivia, with three species belonging to the Lactobacillus genus. Across all three groups, L. iners was the predominant species, exhibiting a trend of initial decline followed by an increase. In comparison to HN, L. gasseri levels significantly decreased following HPV infection. No statistically significant differences were observed in the levels of G. vaginalis, A. vaginae, and P. bivia between HTI and HN; however, these species were significantly reduced in HPI. Conversely, S. amnii levels significantly increased in HPI (Fig. 1B).

Characteristics of vaginal microbiota profiles across diverse populations

To comprehensively evaluate the overall architecture of cervical-vaginal microbiota, a statistical analysis was conducted on 151 samples, categorizing them based on their CSTs (Fig. 1C). Notably, neither CST I, which is characterized by L. crispatus as the predominant species, nor CST V, dominated by Lactobacillus jensenii, were identified in this study. Instead, among the 151 samples, only 5 (representing 3.31%) were classified as CST II, with L. gasseri serving as the dominant species. Additionally, 52 samples (representing 34.44%) belonged to CST III, which is characterized by L. iners as the primary species.

CST IV, on the other hand, is characterized by a mixed microbiota lacking a specific dominant species. This type typically encompasses a diverse array of bacterial species, including G. vaginalis, P. bivia, and A. vaginae, reflecting a highly diversified microbial community state. In this study, the vast majority of samples fell into CST IV, with a total of 94 samples (representing 62.25%). Intriguingly, the distribution of samples between HTI and HPI was remarkably similar in both CST III and CST IV.

Comparative analysis of vaginal microbiota richness and α-diversity across diverse HPV infection statuses

At OTU level, we leveraged the Chao1 index and Shannon index to evaluate the richness and diversity of vaginal microbiota, respectively. Our findings revealed no statistically significant difference (P > 0.05) in Chao1 index values between HN (150.90 ± 15.53) and HTI (167.87 ± 14.39). However, a marked and statistically significant difference (P < 0.0001) emerged when comparing to HPI (225.18 ± 14.77). Furthermore, a highly significant difference (P = 0.004) was also noted between HTI and HPI (Fig. 2A).

Fig. 2.

Vaginal microbial richness and diversity. Note: A is a box plot showing the difference in Chao1 index between different HPV groups; B is a box plot demonstrating the difference in Shannon index between different HPV groups.

Conversely, the Shannon index indicated no significant differences (P > 0.05) in microbial diversity among HN (2.72 ± 0.17), HTI (2.34 ± 0.15), and HPI (2.46 ± 0.13) (Fig. 2B). These results collectively suggest that HPV infection is associated with an increase in vaginal microbiota richness, yet the proportional distribution of dominant species remains relatively unchanged across the different infection statuses.

Compositional differences and β-diversity of vaginal microbiota across diverse HPV infection statuses

To delve deeper into the diversity of microbial community compositions among various HPV infection statuses, we constructed a phylogenetic tree utilizing OTU sequences and employed the Unweighted UniFrac metric to compute the UniFrac distance. The computed dissimilarity coefficients were as follows: 0.672 between HN and HTI, 0.728 between HN and HPI, and 0.709 between HTI and HPI. These coefficients underscore substantial variations in microbial community compositions among the three groups (Fig. 3A).

Fig. 3.

Dissimilarity coefficient calculation based on unweighted UniFrac and NMDS analysis. Note: A presents the correlation analysis among samples grouped by different HPV statuses. The numbers in the squares represent the dissimilarity coefficients between pairs of samples. The smaller the dissimilarity coefficient between two samples, the smaller the difference in species diversity. B displays the NMDS analysis of community differences among grouped samples.

Furthermore, NMDS analysis of community differences revealed a clear segregation among the three sample groups, with a Stress value of 0.15 indicating a good fit of the data to the two-dimensional space. ANOSIM test for assessing community dispersion uniformity confirmed that there were statistically significant differences in microbial communities between HN and both HTI and HPI. Notably, the primary source of this discrepancy was the significant difference observed between HN and HPI (Fig. 3B). These findings collectively highlight that the vaginal microbial community composition exhibits notable variations across different HPV infection statuses.

Identification of vaginal microbial biomarkers

To ascertain the taxa exhibiting relatively high abundance across HN, HTI, and HPI, and to identify statistically significant biomarkers, we utilized LEfSe with a significance threshold of α = 0.05 (Fig. 4A). At the species level, our analysis revealed three differential microbial biomarkers: L. gasseri, A. vaginae, and L. jensenii (Fig. 4B). Meanwhile, dot-box plots were employed to visually illustrate the distribution of absolute abundances of the biomarker bacterial species across each sample group (Fig. 4C–E).

Fig. 4.

Microbial differences and abundance analysis between groups. Note: A, B utilize LEfSe to identify biomarkers that distinguish between samples grouped by different HPV statuses; C–E represents the distribution of absolute abundances of corresponding biomarkers for each sample within the three sample groups, namely HN, HTI, and HPI. F–H presents volcano plots of pathogenic microorganisms with significant differences between pairs of sample groups; I displays a heatmap of changes and differences in microbial abundance at the top 10 species level among different HPV groups.

To gain deeper insights into the differences in species abundance among these groups, we employed the MetaStat method and generated a volcano plot based on the P-values of the differential species. When compared to HN, HTI exhibited a total of 92 differential microbial species. Notably, species such as Bacteroides sartorii, Lactobacillus coleohominis, Lachnospiraceae bacterium, and Weissella koreensis displayed a significant increase in abundance (P < 0.01), whereas species like Neisseria perflava, Streptococcus parauberis, Anaerostipes hadrus, and Bifidobacterium pseudolongum showed a significant decrease (P < 0.01) (Fig. 4F). In contrast, when compared to HN, HPI demonstrated 164 differential microbial species. Species such as Anaerococcus rubeinfantis, Paenibacillus taiwanensis, Ruminococcus callidus, Bacillus coagulans, and Chromobacterium aquaticum exhibited a significant increase in abundance (P < 0.01), while Corynebacterium maris and Megasphaera micronuciformis displayed a significant decrease (P < 0.01) (Fig. 4G). Additionally, the abundance of P. melaninogenica was notably increased in both HTI and HPI. When comparing HPI to HTI, 51 differential microbial species were identified. Species such as Bacillus infernus, Mageeibacillus indolicus, Rhabdanaerobium thermarum, and Chromobacterium aquaticum showed higher abundance in HPI, whereas L. coleohominis, Corynebacterium sp., Campylobacter ureolyticus, and Acinetobacter guillouiae exhibited significantly higher abundance in HTI (Fig. 4H).

To investigate the variations and disparities in microbial abundance across diverse HPV infection groups, focusing on the top 10 species, a biclustering heatmap was constructed (Fig. 4I). Among these species, Five microorganisms exhibited significant differences between groups, namely P. melaninogenica, L. delbrueckii, L. iners, A. vaginae, and L. gasseri. Notably, L. delbrueckii and A. vaginae exhibited statistically significant variations exclusively in HPI (P < 0.01). Specifically, L. delbrueckii abundance increased compared to HN, whereas A. vaginae abundance decreased markedly. P. melaninogenica and L. gasseri displayed significant associations with HPV infection (P < 0.05). An upward trend in P. melaninogenica abundance was observed with HPV infection progression, whereas L. gasseri abundance decreased post-HPV infection. However, the differences in abundance between HTI and HPI were relatively minor and substantially less pronounced than those observed in HN, aligning with anticipated findings. Furthermore, a significant difference in L. iners abundance was solely identified between HTI and HPI, with HTI showing substantially lower abundance than HPI. The abundances of the remaining five microorganisms underwent notable fluctuations before and after HPV infection but showed no statistically significant differences between groups (P > 0.05).

Discussion

The female urogenital microecosystem stands as one of the four cardinal ecological systems in the human body, critically contributing to human reproduction and immune defense mechanisms. HPV infections are often accompanied by vaginal inflammation. Our study reveals that, in comparison to those without HPV infection, individuals with the virus demonstrate notable differences: decreased vaginal cleanliness, an elevated positive rate of sialidase, a higher pH level, and a diminished presence of Lactobacilli. These findings align with previous research reports.^23,24 In a healthy vaginal microecosystem, the vaginal flora is primarily dominated by Lactobacilli species, including L. crispatus, L. gasseri, L. iners, and L. jensenii. These Lactobacilli species secrete unique proteins such as lactic acid, hydrogen peroxide, and defensins, effectively inhibiting the growth of other microorganisms. Studies have indicated that cervical-vaginal dysbiosis can weaken the cervical-vaginal barrier function²⁵ and modify metabolic characteristics,²⁶ which may, in turn, facilitate HPV infection and the progression of CIN or cancer.

According to reports, in HN group, L. iners and L. gasseri emerge as the predominant species, with L. gasseri showing a significantly higher proportion compared to women experiencing HTI group or HPI group, potentially linked to a lower risk of CIN.²⁷ However, in our study, we found no samples belonging to CST I or CST V types. The absence of CST I and V in our cohort might partially reflect the limited species-level resolution of 16S rRNA sequencing. Future studies using shotgun metagenomics may improve taxonomic precision. Instead, we observed five CST II samples (5/151, 3.31%) with L. gasseri as the dominant species and 52 CST III samples (52/151, 34.44%) with L. iners predominating. The remaining 94 samples (94/151, 62.25%) fell into CST IV, characterized by the absence of any specific dominant species. The distribution proportions among the groups were largely consistent within CST III and CST IV. This contrasts with the findings of Ravel et al.²⁸ study on Asian populations, which reported CST III as the most prevalent type in the normal group and CST IV as the dominant type in the HPV-positive group. We hypothesize that the observed differences in microbial community distribution may stem from variations between our study population and that of Ravel et al.,²⁸ leading to distinct microbial community compositions.

Research has established a link between an elevated vaginal microbial diversity and both HPV acquisition and persistent infection.²⁹ Nevertheless, it is noteworthy that variations in community composition or a high degree of diversity do not unequivocally signify microbial dysbiosis.³⁰ In our study, after conducting α-diversity analysis, the Chao1 index revealed that women with HPI exhibited significantly greater diversity in vaginal pathogens compared to those with HTI or HN, thereby reinforcing the notion that HPI augments the diversity of vaginal pathogenic microorganisms.^18,31 Moreover, an increase in vaginal microbial diversity may also pose as a risk factor for HPI. Conversely, the Shannon index indicated no discernible difference among HN women, those with HTI, and those with HPI, implying that while HPV infection boosts the richness of vaginal microorganisms, it does not significantly alter the proportion of dominant species across various infection types. Nevertheless, some studies have reported that HPV-positive women exhibit lower richness of vaginal pathogenic microorganisms compared to HN women,³² and metagenomic sequencing has shown no clear distinction in α-diversity between the two groups.³¹ These discrepancies in sequencing outcomes may be attributed to variations in sequencing methodologies.^33,34

β-Diversity analysis further underscored significant differences in microbial communities among HN group, HTI group, and HPI group. Both UniFrac distance calculations and community difference analyses (NMDS and Anosim) highlighted distinct microbial diversity among these three groups. This underscores that the HPV infection status may, to some extent, influence the composition of vaginal microbial communities, with HPI potentially exerting a more enduring impact on the vaginal microbial ecosystem, thereby exacerbating the differences observed in comparison to HN group. Notably, in women with HPI, a minority of strictly anaerobic microorganisms, such as those belonging to the genera Gardnerella, Prevotella, Atopobium, and Sneathia, which dominate in CST IV, can elicit potent inflammation within the vaginal milieu.³⁵

LEfSe analysis pinpointed L. gasseri, A. vaginae, and L. jensenii as the distinctive biomarkers for HN cohort, HTI cohort, and HPI cohort, respectively. When compared to the HN group, the abundance of these three biomarkers significantly declined. Notably, L. gasseri exhibited the highest abundance in the normal group. Although it was also present in the vaginas of patients infected with HPV, its abundance level was markedly reduced. This observation implies that L. gasseri might play a beneficial role in reducing the HPV load. Wang³⁶ conducted an analysis of the vaginal microbiome in women with varying degrees of HPV infection in the South China region. The study revealed a significant reduction in the abundance of L. jensenii among HPV-positive patients with cervical lesions (P < 0.001). This result aligns with the findings of our study, which designated L. jensenii as a biomarker. However, in a study by Laura L,³⁷ which analyzed the vaginal microbiome of African American women in Chicago, no significant association was observed between HPV and the overall relative abundance of the Lactobacillus genus. Similarly, there was no significant correlation with the relative abundance of L. gasseri, L. iners, and L. jensenii individually. We hypothesize that the discrepancy in results may stem from differences in the distribution of the study populations. These differences could be significantly influenced by variations in environmental factors and lifestyle habits. It is noteworthy that A. vaginae, a strictly anaerobic microorganism, demonstrates a substantial predictive capacity for BV.³⁸ To delve deeper, a differential abundance analysis was performed across the detected pathogenic microorganisms in each group, culminating in the creation of a heatmap. Among the top ten species, five microorganisms exhibited significant differences among the groups: P. melaninogenica, L. delbrueckii, L. iners, A. vaginae, and L. gasseri. Specifically, these microorganisms exhibited a higher diversity in microbial abundance within HPI group compared to HN group. Of particular interest, three out of these five species belong to the Lactobacillus genus, indicating a close correlation between the decrement of Lactobacilli and both the incidence and persistence of HPV infections, which concurs with previous research findings.³⁹ Furthermore, Prevotella and Atopobium have been implicated in biofilm formation, and their heightened abundance may facilitate the persistence of the virus.⁴⁰

This study investigated the alterations and correlations of the vaginal microbiota under different HPV infection statuses. It unveiled the potential link between vaginal microbial dysbiosis and HPV infection, and identified three key biomarkers: L. gasseri, A. vaginae, and L. jensenii. Moreover, the findings suggested that a microenvironment dominated by Lactobacilli might offer protective effects against persistent HPV infection. Future interventions targeting the vaginal microbiota could potentially modulate the progression of HPV infection. This research lays a theoretical foundation for the development of HPV risk prediction tools based on microbial biomarkers and microbiota-modulating therapies.

CRediT authorship contribution statement

Min Guo: Writing-original draft, Supervision, Resources, Investigation, Funding acquisition, Data curation. Chaoyang Chen and Wenling Wang: Visualization, Validation, Software, Formal analysis, Conceptualization. Cairong Zhang and Jie Ma: Visualization, Software, Resources, Investigation, Conceptualization. Maidinamu Sadike and Xiaoli Feng: Writing-review & editing, Validation, Supervision, Conceptualization. Mayinuer Niyazi and Kaichun Zhu: Writing-review & editing, Writing-original draft, Validation, Data curation, Conceptualization, Project administration.

Ethics approval

We confirm that our study complies with all relevant national and international regulations governing research involving human subjects. We have taken all necessary steps to protect the privacy and confidentiality of the participants’ information, ensuring that their identities remain anonymous in all published materials.

Funding

Thanks for the support by Special Programme for Key Research and Development Project of Xinjiang Uygur Autonomous Region (Grant Number: 2022B03018-3).

Conflict of interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Kaichun Zhu reports financial support was provided by Special Programme for Key Research and Development Project of Xinjiang Uygur Autonomous Region, China. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.