A comprehensive literature search was conducted across Web of Science, PubMed, Embase, Scopus, Cochrane Library, and ClinicalTrials.gov, covering studies from January 1, 2010, to January 31, 2025. The temporal restriction from 2010 onward was applied to prioritize contemporary randomized evidence reflecting current vaccine formulations, standardized outcome reporting, and populations relevant to present vaccination policies. The search utilized Medical Subject Headings (MeSH) related to papillomavirus infections, HPV vaccines (such as Gardasil and Cervarix), immunization, and randomized controlled trials. Additionally, bibliographic references of the selected studies were examined to identify further relevant articles.

The selection criteria were established following the PICO (Population, Intervention, Comparison, Outcome) framework. The study included both males and females across different age groups with documented HPV serostatus (either seropositive or seronegative) and HIV infection status. The intervention was limited to prophylactic HPV vaccination using one of the following formulations: bivalent (HPV16 and HPV18), quadrivalent (HPV6, HPV11, HPV16, and HPV18), or nonavalent (HPV6, HPV11, HPV16, HPV18, HPV31, HPV33, HPV45, HPV52, and HPV58) vaccines. Vaccination regimens included single and two-dose schedules (administered at varying intervals, such as months 0 and 2, 0 and 6, or 0 and 12) and three-dose schedules, either at standard intervals (e.g., months 0, 1, and 6 or 0, 2, and 6) or extended intervals (e.g., months 0, 2, and 12). The comparison group received either a placebo (e.g., saline solution or an adjuvant-containing formulation) or an active control vaccine (e.g., hepatitis B virus [HBV] or hepatitis A virus [HAV] vaccines), administered at the same dosage and schedule as the intervention group. The full list of main outcomes is provided in Supplementary Table S1.

Studies were excluded if they did not specifically evaluate HPV vaccination, involved inappropriate interventions in either the experimental or control groups, included therapeutic HPV vaccines or prophylactic vaccines administered with a therapeutic intent, were unpublished, did not adhere to a randomized controlled trial (RCT) design, or failed to report at least one predefined main outcome measure.

Titles and abstracts were independently screened by two reviewers Zhiliang Wang and Hui Zhang), with discrepancies resolved by a third investigator (Danhe Yang). Extracted data included the first author’s surname, year of publication, patient characteristics, sample size, and intervention details (dose, type, duration). For studies with multiple experimental arms, data from each arm were extracted separately. When multiple publications reported outcomes from the same trial cohort, data were extracted from the most recent and comprehensive report to avoid duplicate inclusion and to capture longer-term follow-up.

Two independent evaluators (Zhiliang Wang and Hui Zhang) conducted a review of study reliability and bias risk, with a third reviewer (Yuhang Zhu) resolving any discrepancies. The Cochrane RoB 2 tool was employed to assess the quality of randomized controlled trials (RCTs). To examine publication bias, funnel plots were analyzed visually, and Egger’s regression test was performed when the meta-analysis included at least 10 studies. If any asymmetry was detected, the trim-and-fill approach was implemented to correct for potential publication bias. Moreover, a leave-one-out sensitivity analysis was conducted to determine the effect of individual studies on the overall findings.

Subgroup analyses and meta-regression evaluations were performed following predefined selection criteria for meta-analyses that included a minimum of 10 studies.

The strength of the evidence was assessed using the GRADE (Grading of Recommendations, Assessment, Development, and Evaluation) framework, which classifies certainty into four levels: high, moderate, low, and very low. Because all included studies were randomized controlled trials, evidence for each outcome was initially rated as high certainty and subsequently downgraded at the outcome level based on prespecified criteria. Downgrading was applied for risk of bias (informed by RoB 2 assessments), inconsistency (substantial heterogeneity, I² > 50%), imprecision (wide confidence intervals or limited sample sizes), and publication bias (funnel plot asymmetry or significant Egger’s test results when applicable). Indirectness was evaluated by examining the applicability of populations, interventions, comparators, and outcomes to the review question. All certainty judgments followed consistent decision rules to enhance transparency and reproducibility.

Data was analyzed statistically with RevMan 5.4 (Nordic Cochran Centre, Copenhagen, Denmark) and Stata 14 (College Station, TX, USA). This study calculated the pooled risk ratio (RR) for dichotomous data as well as the pooled standardized mean difference (SMD) for continuous data, accompanied by 95% confidence intervals (CIs). Fold changes in geometric mean titers (GMTs) and their corresponding 95% CIs were derived from log-transformed GMT ratios to improve precision and comparability across studies. Inter-study heterogeneity was assessed using I² statistics. A fixed-effects model was applied, whenever there was no significant heterogeneity (I² < 50% or P > 0.05); otherwise, a random-effects model was applied. Reference management and data organization were handled using EndNote 21. A p-value < 0.05 was deemed to be statistically meaningful.

Of 8,453 identified records, 3,045 duplicates were removed, leaving 5,229 for screening. After excluding 3,917 records due to review articles (3,063), animal studies (75), non-English language (129) and unrelated topics (650), 1,312 full-text articles were assessed. Among these, 88 reports were not retrieved due to irrelevance (592), unavailable full-text (12), unpublished data (7), or other reasons (56). After assessing 557 articles for eligibility, 412 were excluded due to undesirable intervention (68), comparison (312), outcome (14), or insufficient data (18). Ultimately, 145 (13–58) studies were included (Figure 1).

This meta-analysis included 145 RCTs published between 2010 and 2024, evaluating the efficacy, safety, and immunogenicity of bivalent, quadrivalent, and nonavalent prophylactic HPV vaccines. Studies assessed single-dose, two-dose, and three-dose regimens across diverse populations. The majority were conducted in China (50 studies) and multinational settings (40 studies). Treatment duration ranged from 1 to 16 months, with follow-up periods extending from 6 to 136 months. While most trials focused on female participants, some included male cohorts. Immunogenicity assessments included both humoral responses, measured by HPV-specific antibody titers, and cellular immunity, evaluating CD4⁺ and CD8⁺ T-cell responses. Safety outcomes encompassed systemic and injection-site adverse events (Supplementary Table S2). Risk of bias assessment indicated that 95% of studies had “Some Concerns,” with six classified as “High Risk,” primarily due to selective outcome reporting (Figure 2).

The meta-analysis of 21 studies demonstrated a notable reduction in CIN I risk with the HPV vaccination compared to the control (RR: 0.15; 95% CI: 0.09 to 0.24; p < 0.00001), though heterogeneity was substantial (I² = 91%) (Figure 3A). Removal of high–risk-of-bias studies yielded effect estimates and heterogeneity comparable to the primary analysis (RR ≈ 0.15; I² ≈ 91%), suggesting that study quality did not account for the substantial heterogeneity observed. The funnel plot indicated a mild asymmetry, suggesting potential publication bias (Supplementary Figure S1A). Egger’s test confirmed bias (p = 0.040), reinforcing the asymmetry detected in the funnel plot (Supplementary Table 3). Post trim-and-fill adjustment, the pooled estimate remained unchanged (Supplementary Table S4), indicating robust results despite bias concerns. The most substantial CIN I reduction was observed in individuals ≤ 20 years old (88% decrease). The three-dose (0/1/6) regimen showed the highest efficacy (87% reduction), while among vaccine types, the quadrivalent vaccine provided the greatest protection (95% reduction) (Supplementary Figure S2). Meta-regression revealed no statistically significant associations across age, dose regimen, or vaccine type, suggesting that none independently accounted for variability in CIN I reduction (Supplementary Table S5).

Meta-analysis of 24 studies showed a significant reduction in CIN II risk following HPV vaccination (80% reduction; 95% CI: 70–87%), with moderate heterogeneity (I² = 74%) (Figure 3B). The funnel plot appeared symmetric, suggesting no substantial publication bias (Supplementary Figure S1B), which was further supported by a non-significant Egger’s test (p = 0.057) (Supplementary Table S3). Subgroup analysis identified a significant effect of dose regimen, with the three-dose (0/1/6) schedule yielding the greatest risk reduction (Supplementary Figure S3). Among vaccine types, the bivalent vaccine demonstrated the highest efficacy in reducing CIN II associated with HPV 16/18. Meta-regression analysis revealed no statistically significant predictors, indicating that none of the examined factors independently accounted for variability in CIN II reduction (Supplementary Table S5).

Analysis of five studies revealed a significant 52% reduction in CIN III risk following HPV vaccination (95% CI: 2–77%). Heterogeneity was moderate (I² = 69%), and the overall effect reached statistical significance (p = 0.04), confirming the protective impact of HPV vaccination against CIN III (Figure 3C).

A comprehensive meta-analysis of 51 studies identified a notable reduction in persistent HPV 16/18 infections following HPV vaccination (84% reduction; 95% CI: 79–88%), with high heterogeneity (I² = 95%) (Figure 4A). The funnel plot indicated asymmetry, suggesting potential publication bias (Supplementary Figure S1C), which was confirmed by Egger’s test (p = 0.008) (Supplementary Table S3). Trim-and-fill analysis did not adjust the pooled estimate, indicating that publication bias had minimal impact on the overall findings (Supplementary Table S4). Subgroup analysis showed no statistically significant differences across evaluated factors, and meta-regression identified no significant predictors, suggesting that none independently accounted for the observed variation in persistent HPV 16/18 infections (Supplementary Table S5 and Supplementary Figure S4).

A detailed meta-analysis involving 24 studies identified a notable decline in incident HPV 16/18 infections following HPV vaccination (75% reduction; 95% CI: 66–81%), with substantial heterogeneity (I² = 92%) (Figure 4B). Funnel plot analysis indicated asymmetry, pointing to the potential presence of publication bias (Supplementary Figure S1D). This finding was further verified through Egger’s test, which yielded a statistically significant result (p = 0.001), indicating the presence of small-study effects (Supplementary Table S3). Trim-and-fill analysis did not adjust the pooled estimate, suggesting that while publication bias was detected, it had minimal impact on the overall findings (Supplementary Table S4). Subgroup analysis identified HIV status as a significant factor, with the HIV-negative group experiencing the greatest reduction in incident infections (94% reduction) (Supplementary Figure S5). No statistically significant differences were observed across other subgroup categories, and meta-regression revealed no significant predictors, indicating that none independently accounted for the variation in incident HPV 16/18 infections (Supplementary Table S5).

Meta-analysis of 10 studies demonstrated a significant reduction in ASC-US risk following HPV vaccination (69% reduction; 95% CI: 58–81%), with high heterogeneity (I² = 94%) (Figure 5A). The funnel plot exhibited slight asymmetry, suggesting potential publication bias (Supplementary Figure S1E). Also, Egger’s test indicated the presence of bias (p = 0.004) (Supplementary Table S3). Trim-and-fill analysis did not alter the pooled estimate, suggesting that although publication bias was detected, its impact on the overall findings was minimal (Supplementary Table S4). Subgroup analysis revealed that the greatest reduction in ASC-US risk was observed in individuals aged 20–30 years (87% reduction), while among HIV status categories, the most pronounced effect was seen in the HIV-negative group (96% reduction) (Supplementary Figure S6). Meta-regression analysis found no significant associations, indicating that none of the examined variables independently influenced ASC-US risk reduction (Supplementary Table S5).

Meta-analysis of 14 studies demonstrated a moderate reduction in HSIL risk following HPV vaccination (18% reduction; 95% CI: 0–32%), with low heterogeneity (I² = 23%) (Figure 5B). The overall effect was borderline significant (p = 0.05). The funnel plot appeared symmetric, suggesting no major publication bias (Supplementary Figure S1F). Egger’s test was not statistically significant (p = 0.573), further confirming the absence of small-study effects (Supplementary Table S3).

A systematic review of 10 studies identified a notable decline in LSIL risk among participants in the experimental group compared to the control group (RR: 0.63; 95% CI: 0.49 to 0.81; p = 0.0003), with substantial heterogeneity (I² = 86%) (Figure 5C). The funnel plot suggested slight asymmetry, raising concerns about potential publication bias (Supplementary Figure S1G). However, Egger’s test (p = 0.178) found no statistically significant indication of publication bias (Supplementary Table S3). Subgroup analysis showed the highest LSIL risk reduction in the HIV-negative group (94%) and the bivalent vaccine group (47%) (Supplementary Figure S7). However, meta-regression analysis across these subgroup criteria showed no statistically significant associations (Supplementary Table S5).

A meta-analysis including six studies found a significant reduction in AIN risk among individuals who received the HPV vaccine, compared to those in the control group (RR: 0.64; 95% CI: 0.51 to 0.80; p < 0.0001), with moderate heterogeneity (I² = 37%) (Figure 5D).

The meta-analysis of three studies indicated a notable decrease in AIS risk among individuals receiving the HPV vaccine in comparison with those in the control group (RR: 0.32; 95% CI: 0.11 to 0.98; p = 0.05), with low heterogeneity (I² = 29%) (Figure 5E).

The meta-analysis of seven studies showed no statistically significant difference in EGL risk between the HPV vaccine group and the control group (RR: 0.66; 95% CI: 0.29 to 1.49; p = 0.31), with moderate heterogeneity (I² = 67%) (Figure 6A).

A systematic review of four studies identified a significant decline in tissue-resident CD4⁺ T cell levels among participants who received the HPV vaccine compared to those in the control group (SMD: -1.27; 95% CI: -2.24 to -0.31; p = 0.010), with substantial heterogeneity (I² = 93%) (Figure 6B).

The meta-analysis of three studies showed no significant difference in tissue-resident CD8⁺ T cell levels between the HPV vaccine group and the control group (SMD: -1.12; 95% CI: -2.75 to 0.51; p = 0.18), with high heterogeneity (I² = 96%) (Figure 6C).

The meta-analysis of four studies demonstrated a 3.48-fold increase in anti-HPV6 antibody levels in the HPV vaccine group compared to the control (Ratio for GMT: 3.48; 95% CI: 2.59 to 4.37; p < 0.00001), with substantial heterogeneity (I² = 100%) (Figure 6D).

The meta-analysis of four studies showed a 3.37-fold increase in anti-HPV11 antibody levels in the HPV vaccine group compared to the control (Ratio for GMT: 3.37; 95% CI: 1.79 to 4.96; p < 0.0001), with substantial heterogeneity (I² = 100%) (Figure 6E).

The meta-analysis of nineteen studies showed a 3.09-fold increase in anti-HPV16 antibody levels in the HPV vaccine group compared to the control (Ratio for GMT: 3.09; 95% CI: 2.16 to 4.03; p < 0.00001), with substantial heterogeneity (I² = 100%) (Figure 7A). The funnel plot did not reveal substantial asymmetry (Supplementary Figure S1H), and Egger’s test confirmed the absence of significant publication bias (p = 0.297) (Supplementary Table S3). Subgroup analysis showed the highest increase in anti-HPV16 antibody levels in the ≥30 years age group (5.13-fold, 95% CI: 3.95–6.31), followed by individuals aged ≤20 years (3.37-fold, 95% CI: 1.58–5.17) (Supplementary Figure S8). Antibody titers were higher in HIV-negative individuals (2.01-fold, 95% CI: 0.64–3.37) compared with HIV-positive individuals (0.88-fold, 95% CI: 0.70–1.05). Moreover, across included studies, higher fold changes were observed in sex-not reported (NR) studies (5.74-fold, 95% CI: 5.73–5.74), followed by female-only studies (2.92-fold, 95% CI: 1.95–3.88). Consistently, anti-HPV16 antibody levels increased across assay types, including enzyme-linked immunosorbent assay (ELISA) (2.65-fold, 95% CI: 1.33–3.97), chemiluminescent Immunoassay (CLIA) (2.50-fold, 95% CI: 2.40–2.61), and neutralization assays (5.13-fold, 95% CI: 3.95–6.31). Although, meta-regression based on these subgroup criteria did not yield significant results (Supplementary Table S5).

The meta-analysis of 20 studies showed a 3.10-fold increase in anti-HPV18 antibody levels in the HPV vaccine group compared to the control (Ratio for GMT: 3.10; 95% CI: 2.28 to 3.92; p < 0.0001), with substantial heterogeneity (I² = 100%) (Figure 7B). Egger’s test and funnel plot suggested potential publication bias (p = 0.005) (Supplementary Table S3 and Supplementary Figure S1I), while trim-and-fill analysis did not adjust the pooled estimate, indicating robustness of the results (Supplementary Table S4). Subgroup analysis indicated that the largest increase was observed in the ≥30 years age group (4.70-fold, 95% CI: 3.99–5.40). Across sex-based categories, studies enrolling both female and male participants demonstrated the greatest statistically significant increase (5.23-fold, 95% CI: 4.01–6.45), followed by studies with sex- NR cohort (5.06-fold, 95% CI: 3.91–6.20), whereas female-only studies showed a comparatively lower increase (2.85-fold, 95% CI: 1.92–3.78). Additionally, the subgroup with unspecified HIV status (NR) exhibited the greatest fold change in antibody titers (5.06-fold, 95% CI: 3.91–6.20). Although these increases were consistently observed across assay types, higher fold changes were detected in neutralization assays (4.70-fold, 95% CI: 3.95–5.40) than in CLIA and ELISA methods (Supplementary Figure S9). Meta-regression based on these subgroup criteria did not yield significant results (Supplementary Table S5).

A meta-analysis encompassing 15 studies found no statistically significant variation in the incidence of Grade 3 adverse events between individuals in the HPV vaccine group and those in the control group (RR: 1.42; 95% CI: 0.76 to 2.64; p = 0.27). The analysis reported high heterogeneity (I² = 94%) (Figure 7C). Funnel plot analysis showed no major asymmetry (Supplementary Figure S1J), while Egger’s test confirmed that significant publication bias was not present (p = 0.565) (Supplementary Table S3). Additionally, subgroup analyses, including stratification by follow-up duration, and meta-regression did not yield statistically meaningful outcomes (Supplementary Table S5 and Supplementary Figure S10).

A meta-analysis involving 27 studies found a 10% decrease in the occurrence of serious adverse events among individuals in the HPV vaccine group compared to those in the control group (RR: 0.90; 95% CI: 0.82 to 0.99; p = 0.03). The analysis reported moderate heterogeneity (I² = 39%) (Figure 8A). The funnel plot exhibited a predominantly symmetrical distribution (Supplementary Figure S1K), and Egger’s test (p = 0.291) confirmed the absence of significant publication bias (Supplementary Table S3).

A meta-analysis of 12 studies indicated a 26% increase in the rate of injection-site adverse events in individuals who received the HPV vaccine compared to those in the control group (RR: 1.26; 95% CI: 1.07 to 1.48; p = 0.007). The analysis found high heterogeneity (I² = 95%) (Figure 8B). The funnel plot revealed a mostly symmetrical pattern (Supplementary Figure S11), while Egger’s test (p = 0.164) confirmed no significant evidence of publication bias (Supplementary Table S3). Subgroup analysis showed that sex-stratified effects were significant (P < 0.00001), with males exhibiting the lowest and only statistically significant risk estimate. Further, among vaccine formulations, the bivalent vaccine was associated with the lowest risk estimate (RR: 1.05) (Supplementary Figure S11). Although no significant variation in local adverse events was observed across follow-up durations. Meta-regression analysis further identified a significant association between vaccine type and injection-site adverse events (β = 0.43; 95% CI: 0.17 to 0.70; p = 0.005) (Supplementary Table S5).

The meta-analysis of 11 studies found no significant difference in the incidence of systemic adverse events between the HPV vaccine group and the control (RR: 0.99; 95% CI: 0.92 to 1.06; p = 0.72), with moderate heterogeneity (I² = 56%) (Figure 8C). The funnel plot appeared symmetric, suggesting minimal risk of publication bias (Supplementary Figure S1M). Egger’s test further confirmed the absence of significant publication bias (p = 0.126) (Supplementary Table S3). Subgroup analysis by sex revealed a significant difference, with studies enrolling both female and male participants showing the lowest risk of systemic adverse events (RR: 0.48; 95% CI: 0.31–0.73), whereas no significant associations were observed in the remaining sex-defined subgroups (Supplementary Figure S12). Although the risk of systemic adverse events did not show a meaningful change with longer follow-up. Meta-regression analysis did not identify any significant associations (Supplementary Table S5).

The meta-analysis showed no significant difference between the HPV vaccine and control groups across seven pregnancy-related outcomes. Elective termination rates were similar (RR: 0.99; 95% CI: 0.93 to 1.05; p = 0.62; I² = 0%). No significant differences were observed for live birth rates (RR: 1.01; 95% CI: 0.95 to 1.06; p = 0.84; I² = 0%) or spontaneous abortion (RR: 1.04; 95% CI: 0.92 to 1.18; p = 0.52; I² = 0%). Similarly, no significant differences were found in ectopic pregnancy (RR: 0.77; 95% CI: 0.44 to 1.36; p = 0.37; I² = 0%), congenital anomalies (RR: 1.05; 95% CI: 0.66 to 1.68; p = 0.83; I² = 1%), or stillbirth (RR: 0.67; 95% CI: 0.35 to 1.31; p = 0.24; I² = 0%). Lastly, preterm birth (PTB) rates showed no significant variation between groups (RR: 0.78; 95% CI: 0.51 to 1.20; p = 0.26; I² = 3%) (Figures 9A–E, 10A, B).

The sensitivity analysis revealed that effect sizes remained largely stable following the exclusion of individual studies, with a few notable exceptions. For ASC-US, omitting Lehtinen 2012_TVC and Lehtinen 2012_TVCN adjusted the pooled RR from 0.31 (95% CI: 0.22–0.42; Z = 7.23; p < 0.00001) to 0.14 (95% CI: 0.06–0.35; Z = 4.25; p < 0.0001). Similarly, for LSIL, removing these studies modified the RR from 0.63 (95% CI: 0.49–0.81; Z = 3.66; p = 0.0003) to 0.14 (95% CI: 0.04–0.47; Z = 3.18; p = 0.001) (Supplementary Figures S13-S14).

According to the GRADE assessment, the quality of evidence varied across outcomes. HSIL, AIN, AIS, serious adverse events, elective termination, live birth, spontaneous abortion, ectopic pregnancy, congenital anomalies, stillbirth, and PTB were rated as high. CIN II, CIN III, LSIL, EGL, TCD4+, GMT6, GMT11, GMT16, injection-site adverse events, and systemic adverse events were classified as moderate. CIN I, persistent infections, incident infections, ASC-US, TCD8+, GMT18, and Grade 3 adverse events were categorized as low-quality (Supplementary Table S6).