The impact of cervicovaginal microbiome on HPV infection in vaccinated young Korean women

Human papillomavirus (HPV) is considered the oldest known sexually transmitted (STI) pathogen, with phylogenetic evidence of infection in Homo sapiens as early as 618,500 years ago [ 1 ]. HPV infection remains the most common STI worldwide and a major risk factor for cervical cancer (CC) [ 2 ], with 604,000 new cases and 342,000 deaths reported in 2020 [ 3 ]. HPV is classified by oncogenicity into high-risk HPV (hrHPV) that causes malignancies and low-risk HPV (lrHPV) that is primarily associated with benign lesions [ 4 ]. Although most HPV infections resolve spontaneously within 12 months, persistence occurs in 10%–20% of cases, creating opportunities for oncogenic mutations [ 5 ]. CC is the only malignancy preventable by vaccination [ 6 ]. The quadrivalent HPV vaccine, targeting types 6, 11, 16, and 18, has been widely implemented through NIP, leading to a reported decrease in CC incidence [ 7 ]. In Sweden, one of the first countries to introduce NIP in 2009, has shown a remarkable 88% reduction [ 8 ]. The U.S. has recommended HPV vaccination since 2005, and data from 2001–2017 show a decline in CC incidence [ 9 ]. In Korea, where CC was once the most common cancer among women, the overall incidence declined from 8.6% in 1999 to 3.7% in 2020 following the introduction of a national cervical screening program in 1999 and HPV vaccination through the NIP in 2016 [ 10 ]. However, the number of hospital visits with CC as the chief complaint among South Korean women aged 20 to 39 years increased 1.2-fold between 2016 and 2020, from 14,572 to 17,806 [ 11 ]. The cervicovaginal tract is colonized by various microorganisms. Together with host epithelial and immune components, these microorganisms play a key role in STI defense [ 12 ]. While their role was previously unclear, recent sequencing technologies have defined microbial diversity as the ‘microbiome’ [ 13 ]. Lactobacillus species are the most abundant bacteria in the cervicovaginal microbiome (CVM), exert protective effects through the production of lactic acid and hydrogen peroxide (H 2 O 2 ) [ 14 ]. However, their composition varies with age, ethnicity, geography, hormonal status, delivery history, sexual lifestyle, and medication use [ 15 ]. Ravel et al. [ 16 ] classified CVM into 5 community state types (CSTs): CST I ( L. crispatus ), II ( L. gasseri ), III ( L. iners ), and V ( L. jenseni i) are dominated by Lactobacillus species whereas CST IV is defined by markedly reduced Lactobacillus predominance and enrichment of anaerobic bacteria such as Gardnerella , Prevotella , Atopobium , and others commonly referred to as bacterial vaginosis-associated bacteria (BVAB) [ 16 17 ]. This categorization has provided a foundation for research on the role of CVM in various obstetric and gynecologic diseases. Although HPV vaccination has been widely implemented, the incidence of CC continues to rise among young Korean women [ 18 ]. Given the dynamic nature of the CVM in this age group, this trend suggests that factors beyond vaccination coverage may contribute to cervical carcinogenesis. Recent studies have reported that vaginal dysbiosis increases the risk of hrHPV persistence and cervical dysplasia, underscoring the role of the CVM in HPV-related disease progression [ 19 ]. Korean women of reproductive age, however, represent the first nationwide HPV-vaccinated generation and have not been investigated in this context despite a rising incidence of CC. We therefore hypothesized that CVM composition differs between hrHPV-positive and hrHPV-negative women in this vaccinated population and aimed to investigate hrHPV prevalence and CVM profiles using 16S rRNA sequencing in Korean women aged 20–35 years.

A total of 79 women were initially recruited, and 72 samples were included in the final analysis after excluding cases with low cellularity or insufficient DNA quality ( Fig. 1 ). Baseline demographic and clinical characteristics according to hrHPV infection status are summarized in Table 1 . The hrHPV-positive and hrHPV-negative groups were comparable in age, BMI, menstrual timing, and parity. However, women with hrHPV infection showed a trend toward higher rates of smoking (46.9% vs. 10.0%), abnormal cytology (43.8% vs. 7.5%), and multiple HPV infections (68.8% vs. 2.5%) (p<0.001). Values are presented as mean ± standard deviation or median (interquartile range) or number (%). ASC-US, atypical squamous cells of undetermined significance; BMI, body mass index; HPV, human papillomavirus; hrHPV, high-risk human papillomavirus; HSIL, high-grade squamous intraepithelial lesion; LMP, last menstrual period; lrHPV, low-risk human papillomavirus; LSIL, low-grade squamous intraepithelial lesion; NILM, negative for intraepithelial lesion or malignancy. * p - values were calculated using the Fisher’s exact test for categorical variables and Welch’s t test for continuous variables (<0.05 was considered statistically significant), but these findings should be interpreted with caution due to the small number of events. † The hrHPV-positive group included women with any hrHPV detection, whereas the hrHPV-negative group included women who were either negative for HPV or positive for only lrHPV types. Of the 72 participants, 36 (50.0%) had HPV infection of any type, including 15 (20.8%) with single-type and 23 (31.9%) with multiple-type infections. Of the 32 women with hrHPV infection, 8 (25.0%) were co-infected with at least one lrHPV type. In contrast, only one participant in the lrHPV group had multiple-type infections. The most frequently detected hrHPV types were HPV 39, 51, 52, 53, 56, and 58, which tended to exhibit higher microbial diversity compared with vaccine-covered types ( Fig. 2 ). Among vaccine-covered types, HPV 16 was detected in 4 samples (5.6%), all in the hrHPV-positive group, whereas HPV 6 was found in one hrHPV-negative case; no infections with HPV 11 or HPV 18 were identified. The most common reasons for hospital visits were routine gynecological checkups (n=23), followed by complaints related to endometriosis or dysmenorrhea (n=15). HPV, human papillomavirus; hrHPV, high-risk human papillomavirus. CVM diversity, as measured by the Shannon index, was significantly higher in women with abnormal cytology results (n=17), including high-grade squamous intraepithelial lesions, low-grade squamous intraepithelial lesions and atypical squamous cells of undetermined significance (p<0.05; Fig. 3A ). At the genus level, Megasphaera , Bulleidia , Turicibacter , Fannyhessea , Sneathia , and Corynebacterium were significantly enriched in the abnormal cytology group based on Wilcoxon rank-sum testing (p<0.05; Table S1 ). At the species level, not only BVAB such as Gardnerella leopoldii but also L. iners and mulieris were enriched in the abnormal cytology group (p<0.05; Table 2 ). hrHPV, high-risk human papillomavirus; N, normal cytology; Y, abnormal cytology. Values are presented as mean ± standard error of the mean. CI, confidence interval. * p<0.05 was considered statistically significant. † Normal cytology corresponds to negative for intraepithelial lesion or malignancy; abnormal cytology includes atypical squamous cells of undetermined significance, low-grade squamous intraepithelial lesion, and high-grade squamous intraepithelial lesion. As shown in Fig. 2 , violin plots revealed that non-vaccine hrHPV types exhibited markedly broader inter-individual variation in Shannon diversity, and several subtypes, particularly HPV 39, 52, and 56, showed significantly higher diversity than vaccine-covered types after FDR correction (q<0.05). At the species level, distinct differences in the CVM were also observed between women positive for hrHPV (n=32) and those negative for hrHPV (n=40). L. iners and BVAB were more frequently detected in the hrHPV-positive group ( Fig. S1A ), whereas L. crispatus was more commonly observed in the hrHPV-negative group ( Fig. S1B ). At the genus level, Gardnerella, Neisseria, Blastomonas, and Pelomonas were significantly enriched in hrHPV-positive women (p<0.05; Table S2 ). Alpha diversity, as measured by the Shannon index, was significantly higher in the hrHPV-positive group compared to the hrHPV-negative group (p=0.035; Fig. 3B ). In the hrHPV-positive group, bacterial genera commonly associated with BVAB, including Gardnerella and Neisseria , were significantly more abundant compared to the hrHPV-negative group (p<0.05; Table S2 ). At the species level, G. leopoldii , G. vaginalis , and M. lornae were significantly enriched in the hrHPV-positive group (p<0.05; Table 3 ), along with L. iners , L. mulieris , and L. jensenii . Among these Lactobacillus species, L. iners showed the highest relative abundance in the hrHPV-positive group (32.2±6.5), whereas L. mulieris and L. jensenii were present at much lower levels (6.1±2.6 and 3.1±1.4, respectively), corresponding to approximately 5-fold and 10-fold lower abundances than L. iners ( Table 3 ). Values are presented as mean ± standard error of the mean. CI, confidence interval; hrHPV, high-risk human papillomavirus. * p<0.05 was considered statistically significant. † The hrHPV-negative group included women negative for HPV or positive for only low-risk HPV types, whereas the hrHPV-positive group included women with any hrHPV detection. A similar trend was observed in heatmap analysis, which showed that Lactobacillus -dominant profiles were associated with lower microbial diversity, whereas CST IV profiles, including those with Gardnerella and Megasphaera , were linked to increased diversity ( Fig. S2A ). Among species, CST IV components such as G. leopoldii , G. vaginalis , and M. lornae positively correlated with diversity, while L. crispatus showed a strong negative correlation (p<0.001). L. iners demonstrated a positive correlation with diversity, while other Lactobacillus species, including L. gasseri , showed weaker or no significant correlations ( Fig. S2B ). Although the genus Lactobacillus was not differentially enriched between groups, species-level analysis showed that L. iner s was consistently associated with abnormal cytology, hrHPV positivity and increased microbial diversity. To further explore the relationship between hrHPV infection, cytology status, and microbiome diversity, a brief descriptive subgroup analysis was performed ( Table S3 ). Although the subgroup size was small and statistical testing was not feasible, women with both hrHPV positivity and abnormal cytology (group A) showed a clearly dysbiotic and highly diverse microbiome, whereas hrHPV-negative women with normal cytology (group D) exhibited relatively uniform and low diversity. These descriptive findings suggest that combined hrHPV infection and cytological abnormality may be associated with greater CVM disruption.

Although prophylactic HPV vaccination has markedly reduced infections by vaccine-covered types such as HPV 16 and 18 worldwide, the persistence of hrHPV infection remains a concern. Several population-based studies have reported increasing detection of non-vaccine hrHPV types, including HPV 39, 51, 53, 56, and 58, suggesting possible ecological shifts in genotype distribution [ 21 22 ]. In Korea, the NIP has included bivalent and quadrivalent vaccines for 12-year-old girls since 2016. As of 2020, the vaccination rate among girls under 15 years was 72% [ 23 ]. Despite this relatively high coverage, our cohort showed an hrHPV prevalence of 44.4%, which is higher than that reported in Denmark (35%) [ 24 ]. This discrepancy suggests regional or ecological differences in post-vaccination HPV type distribution and raises an important question: why does hrHPV infection remain common even after widespread vaccination? HPV infection after vaccination, particularly hrHPV, cannot be explained solely by alterations in the CVM [ 25 ]. Alternative mechanisms such as type displacement—the ecological substitution of vaccine-covered types by non-vaccine types under immune pressure —or mask effects, in which non-vaccine types become detectable once dominant vaccine types are eliminated, may also contribute [ 26 ]. For example, a Swedish study demonstrated that although the quadrivalent vaccine reduced targeted HPV types, the prevalence of non-vaccine hrHPV remained high [ 22 ]. In our cohort ( Table 1 ), however, hrHPV-positive women showed more frequent reinfection with vaccine-covered types, including HPV 16 (n=4, 12.5% vs. n=1, 2.5%), and a markedly higher rate of multiple-type infections and microbial diversity (n=22, 68.8% vs. n=1, 2.5%). Furthermore, the 2 groups differed substantially in both overall diversity and microbiome composition. These findings suggest that post-vaccination hrHPV detection may result not only from mask or type displacement but also from interactions with the CMV that may influence HPV persistence. Consistent with previous studies, our findings support the concept that the CVM plays a critical role in hrHPV persistence beyond viral or host factors. In a large Chinese cohort of 448 women, G. vaginalis was significantly associated with hrHPV infection (72.7%) compared with Lactobacillus -dominant microbiota (44.7%) [ 27 ]. Similarly, in our cohort, hrHPV-positive women exhibited greater microbial diversity and enrichment of BVAB species. These taxa are known to induce local inflammation, disrupt the mucosal barrier, and promote viral persistence by creating a permissive microenvironment [ 28 ]. Taken together, these findings suggest that hrHPV infection after vaccination may not simply reflect viral re-emergence but rather an ecological imbalance within the cervicovaginal niche. The role of L. iners in hrHPV infection remains controversial. Some studies have suggested that L. iners exerts protective effects similar to L. crispatus [ 29 ], whereas others have linked it to higher rates of HPV infection and persistence [ 30 ]. L. iners is unique in that it lacks the ability to produce D-lactic acid and hydrogen peroxide and expresses inerolysin, a pore-forming cytolysin similar to Gardnerella’s vaginolysin, which may compromise epithelial integrity [ 14 31 ]. A U.S. cohort study also reported that L. iners dominance was significantly associated with HPV positivity [ 32 ]. In our study, L. iners was a distinctive feature of hrHPV-positive women, supporting the hypothesis that L. iners represents a transitional or borderline pathogenic species that may facilitate viral persistence under specific ecological conditions. This study is the first in Korea to investigate the relationship between HPV infection and the CVM using 16S rRNA sequencing in vaccinated women aged 20–35 years—a population representing the post-vaccination generation at risk for cervical disease. Women with hrHPV infection after vaccination exhibited a distinctive CVM pattern characterized by increased diversity, enrichment of BVAB, and L. iners dominance. These findings may provide a basis for developing next-generation vaccination strategies that incorporate microbiome modulation and for predicting vaccine responsiveness or identifying candidates who may benefit from booster vaccination, with important clinical implications for individualized HPV prevention. There are several limitations in our study. First, because of its cross-sectional design, this study could not determine causal relationships between CVM composition and hrHPV infection after vaccination. Longitudinal and molecular studies are needed to clarify this temporal relationship. Second, the relatively small sample size and single-center design may limit the generalizability of our findings, while reliance on self-reported vaccination history could have introduced recall bias. Another limitation is the absence of an unvaccinated control group. This study focused on vaccinated women to explore post-vaccination HPV infection and microbiome features; however, including unvaccinated women in future studies would help to better distinguish vaccine-related ecological and microbial effects. In conclusion, hrHPV infection remains common among vaccinated Korean women and is accompanied by distinctive CVM changes. These findings emphasize the need to consider microbial and ecological factors when refining HPV prevention strategies in the post-vaccination era.

Participants were recruited from women who visited the Department of Obstetrics and Gynecology at Soonchunhyang University Seoul Hospital for routine gynecologic checkups or follow-up for benign conditions between July 1, 2022, and October 30, 2023. Inclusion criteria were as follows: 1) Korean women aged 19–35 years residing in Seoul, and 2) those who had completed HPV vaccination at least 1 year before sample collection. Exclusion criteria included: 1) age under 18 or over 36 years; 2) recent use (within 3 months) of antibiotics, probiotics, or vaginal cleansers; 3) immunosuppression due to medications or underlying disease; and 4) a known history of abnormal cervical cytology or positive HPV results, or CC. To minimize potential confounders, participants who met any of the exclusion criteria were not enrolled. All vaginal samples were collected by a single gynecologist using a standardized protocol [ 20 ]. The following clinical variables were recorded: body mass index (BMI), recent sexual intercourse within 6 months, days from the last menstrual period to sample collection, parity, marital status (married or unmarried), and smoking status. Information on contraceptive methods was not obtained due to the potential unreliability of self-reported data. Nucleic acid amplification tests for Chlamydia trachomatis and Neisseria gonorrhoeae were selectively performed only in symptomatic women, and vulvovaginal infections were identified in 7 participants (5 with Gardnerella spp. and 2 with Candida spp.), who were included for analysis. A total of 79 women were initially recruited, and 72 samples were included in the final analysis ( Fig. 1 ). This cross-sectional study aimed to identify trends in CVM composition according to hrHPV infection status within a region- and population-specific cohort. Although a formal sample size calculation was not performed, the final cohort (32 hrHPV-positive and 40 hrHPV-negative) provided sufficient sequencing depth and group balance for comparative analysis of microbial diversity. The study protocol and informed consent form were reviewed and approved by the Institutional Review Board (IRB) of Soonchunhyang University Hospital, Seoul (IRB No. 2022-04-005). HPV genotyping was performed using the Anyplex™ II HPV28 Detection system (Seegene, Seoul, Korea), which identifies 19 hrHPV types (16, 18, 26, 31, 33, 35, 39, 45, 51, 52, 53, 56, 58, 59, 66, 68, 69, 73, and 82) and 9 lrHPV types (6, 11, 40, 42, 43, 44, 54, 61, and 70). HPV DNA was extracted using the STARlet automated system (Hamilton, Reno, NV, USA) in accordance with the manufacturer’s instructions. PCR amplification was carried out in a 20 μL reaction mixture comprising 5 μL of DNA template, 5 μL of 4× HPV28 TOM, 5 μL of RNase-free water, and 5 μL of Anyplex PCR Mastermix. Real-time PCR was performed using the CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The protocol included uracil-DNA glycosylase activation (50°C, 4 minutes), initial denaturation (95°C, 15 minutes), and 50 amplification cycles (95°C for 30 seconds, 60°C for 1 minute, and 72°C for 30 seconds), followed by cyclic catcher melting temperature analysis at cycles 30, 40, and 50. Melt curve analysis was performed from 55°C to 85°C with a ramp rate of 0.5°C/5 seconds. HPV detection was semi-quantitatively categorized as “+++” (detected at 30 cycles), “++” (at 40 cycles), and “+” (at 50 cycles). HPV genotypes were interpreted using SeegeneViewer software (Seegene). The 16S rRNA sequencing library was prepared following the Illumina 16S Metagenomic Sequencing Library Preparation protocol. The V3–V4 hypervariable regions of the bacterial 16S rRNA gene were amplified from 5 ng of genomic DNA using universal primers with Illumina overhang adapter sequences. The PCR reaction included 5× reaction buffer, 1 mM dNTP mix, 500 nM of each primer, and Herculase II fusion DNA polymerase (Agilent Technologies, Santa Clara, CA, USA). The initial PCR (25 cycles) consisted of denaturation at 95°C for 3 minutes, followed by 25 cycles of 95°C for 30 seconds, 55°C for 30 seconds, and 72°C for 30 seconds, with a final extension at 72°C for 5 minutes. PCR products were purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA). A second PCR was performed to add Nextera XT index primers, using the same thermal conditions but with 10 cycles. The indexed PCR products were purified with AMPure XP beads and quantified using quantitative PCR (KAPA Library Quantification Kit; Kapa Biosystems, Wilmington, MA, USA). Quality assessment was performed using the TapeStation D1000 ScreenTape (Agilent Technologies, Waldbronn, Germany). Sequencing was conducted on the Illumina MiSeq™ platform (Illumina, San Diego, CA, USA). Raw sequencing data generated by the Illumina MiSeq platform were demultiplexed using index sequences to produce paired-end FASTQ files for each sample. Adapter sequences and target-specific primers were trimmed using Cutadapt (v3.2), and reads were truncated to 250 bp for forward (Read 1) and 200 bp for reverse (Read 2) sequences. Sequence denoising and amplicon error correction were performed using the DADA2 package (v1.18.0) in R (v4.0.3). Paired-end reads with expected errors ≥2 were excluded, and an error model was constructed for each batch to correct sequencing noise. Forward and reverse reads were merged, and chimeric sequences were removed using DADA2’s consensus method to generate amplicon sequence variants (ASVs). ASVs shorter than 350 bp were filtered out using R. Each ASV sequence was aligned against the NCBI 16S Microbial Reference Database using BLAST+ (v2.9.0), and taxonomy was assigned based on the best hit. Sequences with <85% query coverage or identity were excluded from taxonomic assignment. Multiple sequence alignment of ASVs was performed using MAFFT (v7.475), and a phylogenetic tree was constructed using FastTreeMP (v2.1.10). Alpha and beta diversity analyses were performed using QIIME. For alpha diversity, Shannon and Inverse Simpson indices were calculated, and richness was assessed using rarefaction curves and Chao1 estimates. Beta diversity was evaluated based on weighted and unweighted UniFrac distances and Bray–Curtis dissimilarity. Principal coordinates analysis and unweighted pair group method with arithmetic mean clustering were used to visualize microbial community differences. Group-level significance was assessed using analysis of similarity. The Wilcoxon signed-rank test were used to assess differences in alpha and beta diversity and in taxonomic composition at the phylum and genus levels. Differentially abundant taxa were identified accordingly. Spearman correlation analysis was performed to evaluate associations between alpha diversity indices and taxon-level relative abundances across all samples. To adjust for multiple testing, false discovery rate (FDR) correction was applied using the Benjamini–Hochberg method. Data visualization was conducted using the ggplot2 (v3.2.1) package in R (v4.0.3), including scatter plots and heatmaps.