In April 2025, we searched Pubmed, Embase, Science Direct, and Web of Science databases using keyword free search mode, including keywords such as “Cervical neoplasms” or “Cervical cancer” or “Hormone replacement therapy” or “HRT”. We also searched for related topics on Google Scholar.

Research type: all studies were observational studies including cohort studies which can be either retrospective or prospective, or case-control studies, or cross-sectional observational studies.

Research subjects: the research subjects were women who undergo HRT for premenopausal or postmenopausal women (including natural menopause, surgical or pharmacological menopause).

Object: the correlation between exposure (HRT) and outcome (risk of CC), should be studies in literature.

Outcomes: studies should provide the OR of HRT treatment relative to the risk of CC.

Duplicate records were identified by using Endnote 20 (Build 16480) initially. Two researchers independently screened titles and abstracts based on predefined inclusion criteria. Discrepancies were resolved through structured discussion following objective criteria (1): re-examination against specific inclusion/exclusion criteria (2); consultation of the study protocol for borderline cases (3); involvement of a third reviewer when consensus could not be reached. All decisions and rationales were documented with justifications. The full texts of the potential studies (usually a PDF text) were retrieved one by one either from the online databases or by directly contacting the authors. Studies without available full texts were excluded after systematic attempts to obtain them through multiple strategies, including direct contact with corresponding authors, institutional library services, and inter-library loan requests. This exclusion of 11 studies (16.4% of studies sought for retrieval) represents a potential source of selection bias that may affect the generalizability of our findings. Then, the full texts were reviewed for further eligibility. Studies without outcomes or data would be excluded.

The Newcastle Ottawa Scale (NOS) (12) was used to assess the quality of the included studies. The scale evaluated three domains of observational studies including the object selection, comparability, and outcome indicators, with a maximum score of 9 points. Studies with a total score of 5 or above were of acceptable quality, score of 5–7 were of moderate quality, and 8–9 of high quality.

Main outcome: the combined OR of HRT (current and persist use) for cervical cancer incidence.

Secondary outcomes: the combined OR of hormone replacement therapy (current and persist use) for any cytogenic abnormality associated with cervical cancer; The combined OR of past ever use hormone replacement therapy for cervical cancer incidence; The OR of hormone replacement therapy (Duration of use>1 year) for the incidence of cervical cancer; The OR of hormone replacement therapy (age>50) for cervical cancer incidence.

Two researchers independently extracted literature data, including research type, publication date, research location, patient age, main indicators, etc., and then conducted data verification. Discrepancies were resolved through structured discussion with third reviewer consultation when necessary. All extracted data were cross-verified and documented.

The effect of HRT on the incidence of cervical cancer was evaluated using pooled OR and corresponding 95% CI. OR values reported in the literature were first log-transformed, then synthesized using the “metagen” function of the “meta” package, and the results were displayed in a forest plot. Model selection was based on heterogeneity assessment: fixed-effects model (Mantel-Haenszel method) was used when I²≤30% and p≥0.1 for Cochrane’s Q-test, indicating minimal heterogeneity and supporting the assumption of a common true effect size across studies. Random-effects model (DerSimonian-Laird inverse variance method) was used when I²>30% or p<0.1, indicating substantial heterogeneity and allowing for variation in true effect sizes between studies. The choice of I²=30% as a threshold was based on established guidelines that consider I²≤30% as representing low heterogeneity that may not substantially impact the pooled estimate.

I² analysis and Cochrane Q-test were used to detect heterogeneity between literature, with I²>30% or P<0.1 indicating heterogeneity in the results. The interpretation of heterogeneity levels followed established guidelines: I²=0-30% representing low heterogeneity, 31-50% moderate heterogeneity, 51-75% substantial heterogeneity, and 76-100% considerable heterogeneity (13). These thresholds guided our model selection process to ensure appropriate statistical methodology. Subgroup analysis was used to investigate the sources of heterogeneity between literature. Cluster analysis of literature based on heterogeneity and effect size was conducted using GOSH function (14).

The ‘metareg’ function from the ‘meta’ package was called to perform regression analysis, to detect whether the effect size is dominated by factors such as age, publication year, literature quality, and sample size. We simultaneously called the ‘influence’ function of the ‘metafor’ package to conduct influence analysis on the pooled results.

The trim-and-fill method was used to predict the number of documents required to make the funnel fully symmetrical. Egger’s test was used to detect publication bias, and the results were presented in a funnel plot.

All statistical analyses were performed using R software (version 4.4.1). Meta-analyses were conducted using the “meta” package (version 6.5-0) for effect size pooling and forest plot generation. Sensitivity analyses were performed using the “metafor” package (version 4.4-0). The GOSH (Graphical Display of Study Heterogeneity) function was implemented to explore potential clustering patterns among studies. Statistical significance was set at p<0.05 for all analyses.

Figure 1 shows the literature selection process. Initially, 529 articles were retrieved, and after preliminary duplication and screening, the PDF full texts of 56 articles were finally obtained. After detailed screening, 9 articles (10, 11, 14–20) were finally included.

The basic characteristics, grouping information, and outcome indicators of all literature and patients are listed in Table 1 .

Among all 9 articles included, 8 articles (89.9%) had a quality score of 8 or 9, with minimal bias and high quality; One article (1.11%) received a score of 7, with slight bias and moderate quality. The overall quality is good, as shown in Table 2 .

In the pooling of ORs, there was statistical heterogeneity among 14 studies (I² = 49%, p=0.02). We conducted subgroup analysis of the literature according to “Cancer type” and “HRT formula” to investigate the sources of heterogeneity. By dividing the studies into three subgroups according to ‘Cancer type’, the heterogeneity test between groups showed p<0.01, indicating ‘Cancer type’ is a statistically significant source of heterogeneity. The effect size of HRT for the ‘Any type’ cervical cancer subgroup was OR=0.74, 95% CI (0.57, 0.96), for adenocarcinoma subtype was OR=1.82, 95% CI (0.91, 3.65), while the effect size for squamous cell carcinoma subtype was OR=0.51, 95% CI (0.36, 0.71). Mixed type cervical cancer was mentioned in one study but had insufficient sample size for separate statistical analysis. These results suggest that HRT has a protective effect against cervical cancer overall and squamous cell carcinoma specifically but may increase the risk of adenocarcinoma subtype. According to the HRT formula, the studies can be divided into three subgroups with the following effect sizes: E alone OR=0.69, 95% CI (0.55, 0.88), E+P combined OR=0.69, 95% CI (0.48, 0.99), and P alone OR=1.66, 95% CI (0.39, 7.03). The heterogeneity test between subgroups showed p=0.50, indicating no significant heterogeneity between subgroups. These results demonstrate that both estrogen alone and estrogen plus progestin formulations have protective effects against cervical cancer, while progesterone alone shows a non-significant trend toward increased risk ( Figure 3 ).

We also conducted subgroup analysis on the literature according to “Study design” and “Reported indicators” to investigate the sources of heterogeneity. According to the “Study design”, the two subgroups were identified, and the heterogeneity test between groups showed p=0.08, indicating that there was no significant heterogeneity between subgroups, ‘Study design’ is not a statistically significant source of heterogeneity. According to the “Reported Indicator”, the literature can be divided into three subgroups, and the heterogeneity test between subgroups was p=0.23, indicating the “Reported Indicator” is not a statistically significant source of heterogeneity ( Table 3 ).

We conducted additional analyses to examine the effect of specific HRT characteristics on cervical cancer risk ( Figure 4 ). For past use of HRT, analysis of 5 studies revealed a pooled OR of 0.80 (95% CI [0.60, 1.06]), suggesting a protective trend that did not reach statistical significance (p>0.05). This finding indicates that even previous HRT exposure may confer lasting protective effects against cervical cancer development. Duration-specific analysis of 5 studies examining HRT use >1 year showed an OR of 0.81 (95% CI [0.65, 1.01]), demonstrating consistent protective effects regardless of treatment duration, with the upper confidence limit approaching but not exceeding unity. Age-stratified analysis of 6 studies focusing on women >50 years yielded an OR of 0.82 (95% CI [0.54, 1.26]), indicating that HRT’s protective effect persists across different age groups, though with wider confidence intervals reflecting increased heterogeneity in older populations. Notably, all three analyses showed remarkably consistent effect sizes (OR = 0.80-0.82), reinforcing the robustness of HRT’s protective association against cervical cancer across different patient characteristics and exposure patterns.

To investigate the factors that may affect the pooled ES results, we conducted correlation analysis using three sets of data: “publication year”, “sample size”, and “literature quality score” and found a strong correlation of -0.63 between sample size and literature quality ( Figure 5A ). We conducted a multivariate regression analysis using four conventional variables: ‘publication year’, ‘sample size’, ‘literature quality score’, and ‘Continent’ and found that none of these factors had a statistically significant impact on the results (p>0.05). Based on our subgroup analyses, cancer subtype appears to be the primary dominant variable, showing significant between-group heterogeneity (p<0.01). While HRT formulation did not show significant heterogeneity between subgroups (p=0.50), other potentially dominant variables that warrant future investigation include HRT duration, mean participant age, baseline HPV status, and menopausal status, as these clinical and biological factors may have stronger influence on cervical cancer risk than the conventional methodological variables we initially examined. Cancer subtype is currently the only confirmed dominant factor based on our available data, highlighting the need for more comprehensive variable collection in future meta-analyses ( Figure 5B ).

The GOSH diagnosis clearly clustered the studies into 3 groups. Based on the results of previous subgroup analysis, we speculate that these clusters were related to different types of cervical cancer ( Figure 6A ). Influence analysis using multiple diagnostic statistics identified literature (15) [study (3) in Figure 6B ] as an influential outlier across several measures including standardized residuals (rstudent), DFFITS values, Cook’s distance (cook.d), and covariance ratios (cov.r). The study exceeded conventional threshold criteria across these diagnostics, as indicated by the red triangular markers in Figure 6B . However, sensitivity analysis excluding this study showed minimal change in the pooled effect size, confirming the robustness and reliability of our meta-analysis results.

In the analysis of the effect size of HRT (current and persist use) for cervical cancer incidence, 14 studies were assessed using the trim-and-fill method. The result suggested an additional 5 studies need to be added to fully compensate for the asymmetry on both sides of the funnel. However, Egger’s test yielded p=0.386, indicating that the asymmetry on both sides of the funnel is not statistically significant. The unadjusted funnel plots of 14 studies were shown in Figure 7 .