Non-invasive host transcriptome and HPV oncogene expression map the molecular landscape of HPV-driven cervical lesions
Mohamad Ammar Ayass, Naila Zaman, Natalya Griko, Victor Pashkov, Kevin Zhu, Ghulam Abbas, Melesse Ghelan, Ramya Ramankutty Nair, Tutku Okyay, Jin Zhang, Lina Abi-Mosleh
TL;DR
This study uses transcriptomic profiling to identify molecular changes in cervical lesions caused by HPV, offering new biomarkers for better risk assessment and treatment.
Contribution
The study introduces stage-specific transcriptomic signatures and integrates HPV oncogene expression to improve cervical cancer screening.
Findings
NILM showed upregulated growth signaling and immune suppression, indicating stealth infection.
HSIL was marked by cell cycle hyperactivation and epithelial dedifferentiation, showing molecular discontinuity from earlier stages.
Higher E6/E7 mRNA levels correlated with more severe lesions, suggesting their use as biomarkers.
Abstract
Cervical cancer remains a significant global health burden. Current screening methods are not yet capable of detecting molecular alterations preceding cytological abnormalities. In this study, we performed integrative transcriptomic profiling of 132 HPV-positive cervical Pap Smear specimens (NILM, ASCUS, LSIL, HSIL), combining HPV genotyping, and E6/E7 mRNA quantification to map molecular progression. Our analysis revealed stage-specific signatures: NILM displayed a “stealth infection” profile marked by upregulated protein synthesis and growth signaling (EGFR/ERBB2) alongside immune suppression. ASCUS presented a critical tipping point with introduction of early oncogenic drivers (CCND1, SHH), while LSIL prioritized viral productivity with suppressed antimicrobial defenses (MPO, DEFA1). HSIL was distinct from earlier stages and defined by cell cycle hyperactivation (CDK1, PLK1),…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig 1
Fig 2
Fig 3
Fig 4
Fig 5Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsCervical Cancer and HPV Research · Endometrial and Cervical Cancer Treatments · Molecular Biology Techniques and Applications
1 Introduction
Cervical cancer remains a significant global health burden, particularly among women aged 20–39 [1], with an estimated 13,820 new U.S. cases and 4,360 deaths expected in 2024 [2]. Persistent infection with high-risk human papillomavirus (HPV), particularly types 16 and 18, drives nearly 70% of cases by expressing the viral oncoproteins E6 and E7, which disrupt host tumor suppressor pathways such as p53 and Rb [3]. While the majority of HPV infections are transient (67% resolving within 12 months and 90% within two years [4]), persistent infections can progress to precancerous lesions and ultimately, invasive cancer.
Cervical cancer screening guidelines have evolved significantly in recent years to improve early detection while minimizing both overtreatment and missed diagnoses. The United States Preventive Services Task Force (USPSTF) now recommends primary HPV testing every 5 years as the preferred screening strategy for women 30 years and older, with alternatives including HPV/Pap co-testing every 5 years or Pap testing alone every 3 years [5]. The American Cancer Society has taken an even stronger position, recommending primary HPV testing starting at age 25 for all individuals with a cervix. These guideline shifts reflect mounting evidence that HPV testing provides superior sensitivity for detecting precancerous lesions compared to cytology alone. Despite cervical cancer being largely preventable through HPV vaccination and routine screening, cervical cancer detection still faces dual challenges of underdiagnosis and overtreatment. Research has demonstrated that cytology is subjective and offers low sensitivity, missing up to 60% of high-risk infections with progression potential, especially in early stages [6]. Most cervical cancer cases in the U.S. occur in women who have never been screened or who haven’t been screened in over 5 years, highlighting the serious consequences of underdiagnosis. Conversely, while HPV DNA testing offers excellent sensitivity (*>*95%) for detecting high-grade lesions, its positive predictive value remains relatively low (approximately 20–40% for cervical intraepithelial neoplasia grade (CIN2^+^)), resulting in substantial overdiagnosis and unnecessary intervention for many women with transient infections that would resolve spontaneously [7].
The central clinical dilemma in cervical cancer prevention lies in understanding the molecular mechanisms that distinguish HPV infections destined to progress from those that will resolve spontaneously. Current screening methods detect viral presence but provide no functional insight into the critical host-virus interactions at the transcriptome level that determine disease trajectory. This limitation stems from the fact that HPV DNA detection merely confirms viral presence without information about viral activity or host cellular response, critical determinants of progression risk. The consequences of this diagnostic gap are significant: approximately 70% of women referred for colposcopy after positive HPV tests do not ultimately have significant lesions requiring treatment [8], leading to unnecessary procedures, psychological distress, and healthcare costs.
Transcriptomics, the study of the complete set of RNA transcripts produced by the genome, has emerged as a powerful approach for addressing this challenge. By examining the entire transcriptome, researchers can gain holistic insights into cellular processes, regulatory networks, and molecular pathways involved in disease development and progression [9–11]. This approach offers a promising solution by capturing both viral oncogene expression and host cellular changes, potentially enabling more precise risk stratification than current methods allow. The clinical application of transcriptomics offers transformative potential through the development of novel biomarkers with enhanced specificity. For example, quantification of viral oncogene expression (E6/E7) combined with host gene signatures involved in cell cycle dysregulation (p16/Ki-67) or immune evasion (PD-L1/CXCL10 ratios) has demonstrated improved specificity for detecting Cervical Intraepithelial Neoplasia Grade 2 (CIN2^+^) compared to conventional methods. Identification of such biomarkers could potentially reduce unnecessary colposcopies by 30–50% while maintaining high sensitivity for detecting significant lesions [12].
Recent advances in transcriptome profiling have significantly enhanced our understanding of cervical carcinogenesis by mapping dynamic host-virus interactions at the mRNA level. Several studies have identified key gene signatures associated with cervical cancer progression. From cervical biopsies, Yi et al. identified 107 hub genes associated with different stages of cervical intraepithelial neoplasia (CIN1–CIN4), with FN1, ITGB1, and MMP9 linked to tumor cell metastasis. Additionally, STAT1 was found to play a dominant role in phase IV via cancer-associated signaling pathways, while CDK1 and CCNB1 were implicated in regulating proliferation and differentiation through the cell cycle and viral tumorigenesis [13]. Furthermore, IFI6, SLC39A9, and ZNF185 were strongly correlated with tumor progression and patient survival in the OncoLnc database, whereas AKAP12 and DUSP5, previously linked to poor prognosis in other cancers, were identified as novel HPV16 E7-regulated genes with potential roles in early tumorigenesis [14]. Li et al. highlighted the prognostic potential of the histone gene family, demonstrating that genes such as HIST1H2BD, HIST1H2BJ, HIST1H2BH, HIST1H2AM, and HIST1H4K could predict survival outcomes in cervical cancer patients [15]. Other candidate genes, including SPP1, FOXM1, LYN, COL6A3, CCL21, TTK, and MELK have emerged as potential prognostic markers and therapeutic targets for cervical cancer [16].
These findings underscore the molecular complexity of HPV-driven cervical cancer, highlighting dysregulated pathways such as cell cycle control, immune evasion, and epithelial-mesenchymal transition (EMT) that precede cytological abnormalities [17]. Integrating viral mRNA data with host transcriptome profiles enhances risk stratification by capturing both viral activity and the host’s molecular response. For example, combined analysis of E6/E7 mRNA and immune-related transcripts (e.g., CXCL9, PD-L1) has revealed immune suppression mechanisms that contribute to viral persistence and disease progression [18,19].
Immune evasion represents a critical yet under-addressed mechanism through which HPV establishes persistence, the fundamental prerequisite for malignant transformation. High-risk HPV types have evolved sophisticated strategies to subvert host immune surveillance at multiple levels. The virus maintains a completely intraepithelial lifecycle without viremia, minimizing exposure to systemic immune responses. At the molecular level, HPV E6/E7 oncoproteins interfere with interferon signaling pathways by disrupting STAT activation and downregulating interferon-responsive genes. Additionally, viral proteins reduce the expression of pattern recognition receptors (including TLR9) that would otherwise detect viral DNA. E5 protein decreases surface MHC-I expression, hampering antigen presentation and cytotoxic T-cell recognition of infected cells. These evasion mechanisms create a localized immune suppressive microenvironment characterized by reduced pro-inflammatory cytokine production and impaired recruitment of effector immune cells. Transcriptomic profiling across progressive cytological categories offers unique insights into the evolution of this immune evasion from initial subversion of innate immunity in NILM (Negative for Intraepithelial Lesion or Malignancy) samples to profound immunosuppression in higher-grade lesions.
The molecular evolution from initial HPV infection to invasive cervical cancer follows a discernible stepwise progression that occurs prior to and independent of cytological abnormalities. This progression reflects a carefully orchestrated series of viral strategies that unfold sequentially: first, the virus must establish persistence by evading host immune detection [20,21]; second, it drives cellular proliferation to expand the infected cell population [22]; third, it induces genomic instability to facilitate host cell transformation [23,24]; and finally, it promotes invasive potential through disruption of cellular adhesion and differentiation programs. Each step corresponds to specific molecular alterations detectable through transcriptomic analysis, potentially enabling more precise disease classification than morphology-based approaches. Importantly, this molecular progression is not always synchronized with cytological changes, significant transcriptional alterations often precede observable cellular abnormalities by months or years.
Despite these advancements in HPV research, the molecular mechanisms underlying the transition from normal cytology to precancer, a critical window for early intervention, remain not well understood. While previous transcriptomic cancer studies have largely relied on invasive tissue biopsies or cell-lines, our study leverages on non-invasive Pap Smear specimen to characterize early molecular events. In this study, we conducted a comprehensive omics-analysis of host-virus interplay in routine, non-invasive cervical Pap Smear specimens, integrating HPV DNA genotyping, E6/E7 mRNA quantification, and cervical transcriptome profiling of non-invasive cervical Pap Smear samples. Our study systematically maps these molecular transitions across the spectrum from HPV-positive NILM (Negative for Intraepithelial Lesion or Malignancy) through increasingly severe abnormalities (ASCUS (Atypical Squamous Cells of Undetermined Significance), LSIL(Low-Grade Squamous Intraepithelial Lesion), HSIL(High-Grade Squamous Intraepithelial Lesion)), revealing distinct gene expression signatures at each stage. Our findings revealed stage-specific molecular dynamics: in HPV-positive NILM, HPV adopts a stealth infection profile marked by heightened host protein synthesis and growth signaling, alongside suppression of immune detection mechanisms; progression to ASCUS represents a critical tipping point with lowered immune signaling, sustained protein synthesis, and early proliferative activation; LSIL exhibited a virus-centric stage dominated by high productive activity with blunted host immune defense; at the HSIL stage, the host-cell transcriptional landscape becomes overtly oncogenic, with viral hijacking of cell cycle regulation. This structured progression model offers a framework for developing stage-specific biomarkers that could identify high-risk infections before cytological abnormalities emerge, potentially enabling earlier and more targeted intervention.
2 Materials and methods
2.1 Clinical specimens
Cervical cytology specimens preserved in PreservCyt^®^ solution (Hologic, Inc., Marlborough, MA, USA) were procured from multiple commercial biorepositories, including BIOFLUIDS.com (Los Osos, CA), Discovery Life Sciences (DLS, Huntsville, AL), and iSpecimen Inc. (Lexington, MA). All specimens had undergone prior cytological and histopathological assessment to confirm diagnostic status before their inclusion in this study. Each supplier attested that the biospecimens were collected under appropriate institutional protocols with documented informed consent from donors, ensuring ethical compliance for downstream research use. These Pap Smear derived specimens represent a non-invasive sampling method commonly used in routine cervical screening, making them highly suitable for clinical transcriptomic analysis.
All samples used in this study were de-identified Pap smear specimens purchased between April 17, 2024, and April 28, 2025. The research protocol was reviewed and approved by the Salus IRB with a waiver of informed consent, as the study posed no more than minimal risk (Protocol #ABS014 06 01 2025). Data collection and analysis were conducted only after IRB approval on June 27, 2025. At no point during or after the study did the researchers have access to information that could identify individual participants.
2.2 DNA and RNA extraction from clinical specimens
Cellular material from PreservCyt^®^ preserved specimens (Hologic, Inc., Marlborough, MA, USA) was recovered by pelleting. Approximately 3.5 mL of each specimen was used for nucleic acid extraction, with 1 mL allocated for genomic DNA isolation and the remaining 2 mL for total RNA purification. For DNA extraction, samples were centrifuged at 4,000 × g for 10 minutes at 4°C to pellet the cells. After discarding the supernatant, the resulting cell pellet was resuspended in 200 µL of phosphate-buffered saline (PBS). Genomic DNA was then extracted using the GenFind^TM^ V3 Blood, Cell 161 and Serum Genomic DNA Isolation Kit (Beckman Coulter, Brea, CA, USA; Cat. No. 162 C34881), following the manufacturer’s protocol. Purified DNA was stored at −20°C until further use.
RNA extraction followed a similar initial centrifugation step (4,000 × g for 10 minutes at 4°C). After removal of the supernatant, the pellet was immediately resuspended in 1 mL of TRIzol^TM^ Reagent (Thermo Fisher Scientific, Cat. No. 15596018), and total RNA was subsequently isolated using the PureLink^TM^ RNA Mini Kit (Thermo Fisher Scientific, Cat. No. 12183025), following the manufacturer’s “TRIzol Plus Total Transcriptome Isolation” protocol. RNA yield and purity were determined spectrophotometrically using the NanoDrop^TM^ Lite Plus (Thermo Fisher Scientific, Madison, WI, USA), while integrity was assessed via the Qubit^TM^ RNA HS Assay Kit (Invitrogen, Cat. No. Q32855) on a Qubit^TM^ 3.0 Fluorometer. In samples exhibiting visible red blood cell contamination, additional RNA processing was conducted using the GlobinClear^TM^ Human Globin mRNA Removal Kit (Thermo Fisher Scientific, Cat. No. AM1980) to eliminate abundant globin mRNA, which could interfere with downstream transcriptome profiling. All purified RNA samples were stored at −80°C until further downstream application.
2.3 HPV DNA genotyping
Genomic DNA from each sample was evaluated for HPV genotype using 20 parallel quantitative PCR (qPCR) reactions, each incorporating a custom-designed primer–probe set targeting a specific HPV genotype. The panel included high-risk types (HPV16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, and 68) and low-risk types (HPV6, 11, 26, 53, 66, 73, and 82). All primer–probe sets were acquired as 20 × stocks from Thermo Fisher Scientific, with the following specific Assay IDs: HPV16 (AIX02E2), HPV18 (AIY90LA), HPV31 (AI0IYRI), HPV33 (AI1RWXQ), HPV35 (AI20U3Y), HPV39 (AI39S96), HPV45 (AI5IRGE), HPV51 (AI6RPMM), HPV52 (AI70NSU), HPV56 (AI89LY2), HPV58 (AIAA04V), HPV59 (AIBJZA3), HPV68 (APZTGM6), HPV6 (AIGJRZZ), HPV11 (AIFATTR), HPV26 (AIHSP57), HPV53 (AII1OCF), HPV66 (AID1VNJ), HPV73 (AIKAMIN), and HPV82 (AILJKOV) (Supplementary S1 Table). Assays were conducted with appropriate positive controls (TrueMark^TM^ Comprehensive Microbiota Control; Thermo Fisher Scientific, Cat. No. A50382) and internal reference controls (TaqMan^TM^ Copy Number Reference Assay targeting the RNase P gene; Cat. No. 4403328). Human genomic DNA confirmed to be HPV-negative was included as a negative control for each qPCR run.
Reactions were prepared using the TaqMan^TM^ Fast Advanced Master Mix (Thermo Fisher Scientific, Cat. No. 4444557) and performed on a QuantStudio^TM^ 7 Pro Real-Time PCR System. Each 10 µL reaction mixture contained 5 µL of master mix, 0.5 µL of the genotype-specific primer–probe set, 2 µL of DNA template (10 ng/µL), and 2.5 µL of nuclease-free water. Thermal cycling conditions were standardized across all reactions and consisted of an initial incubation at 50°C for 2 minutes, followed by enzyme activation at 95°C for 2 minutes, and then 40 amplification cycles comprising denaturation at 95°C for 5 seconds and annealing/extension at 60°C for 20 seconds.
2.4 Assessment of E6/E7 mRNA expression by RT-qPCR
Quantification of E6 and E7 mRNA expression was performed on RNA samples confirmed to be HPV-positive by prior genotyping. The RT-qPCR assays were conducted using the iTaq^TM^ Universal SYBR^®^ Green One-Step Kit (Bio-Rad, Hercules, CA, USA; Cat. No. 1725151). Primer sequences targeting the E6/E7 regions were adopted from the study by Pan et al., which developed and validated a multiplex assay for high-risk HPV mRNA detection [25]. The specific primers used are listed in Supplementary S1 Table. Each 10 µL RT-qPCR reaction consisted of 5 µL of SYBR Green master mix, 0.125 µL of iScript^TM^ reverse transcriptase, 0.6 µL of a forward/reverse primer mix (0.3µM each final concentration), 2 µL of RNA template (20ng/µL), and 2.3 µL of nuclease-free water. Reactions were run using the following thermal protocol: reverse transcription at 50°C for 10 minutes, polymerase activation at 95°C for 4 minutes, followed by 40 amplification cycles consisting of denaturation at 95°C for 15 seconds and annealing/extension at 60°C for 30 seconds.
2.5 Transcriptome sequencing and gene expression quantification
Transcriptomic profiling was performed on extracted RNA samples following dilution to a concentration ≤ 30ng/µL, with a maximum input of 100 ng per reaction. Library preparation was carried out using the Ion AmpliSeq^TM^ Transcriptome Human Gene Expression Kit (Thermo Fisher Scientific, USA; Cat. No. A26327), according to the manufacturer’s protocol. Amplified libraries were purified using Agencourt^®^ AMPure^®^ XP magnetic beads (Beckman Coulter, USA; Cat. No. A63881) and quantified using the Ion Library TaqMan^TM^ Quantitation Kit (Thermo Fisher Scientific, USA; Cat. No. 4468802) on a StepOnePlus^TM^ Real-Time PCR System.
Eight RNA libraries were pooled per sequencing run, generating a 300 pM pooled library, which was loaded onto the Ion Torrent Genexus^TM^ Integrated Sequencer. Sequencing was performed under the instrument’s “Library-to-Result” automated workflow, enabling integrated template preparation, chip loading, sequencing, and initial data processing. Sequence reads were aligned to the hg19 human reference genome using Torrent Suite^TM^ Software (Thermo Fisher Scientific; version 6.8.1.1), and raw gene expression was quantified using the AmpliSeq RNA plug-in. Gene counts were exported for downstream differential expression and pathway analysis.
2.6 Bioinformatics and computational analyses
Normalization and differential gene expression (DGE) analyses were performed using the DESeq2 R package (v1.36.0). Raw count data were normalized using the median-of-ratios method to control for library size and composition bias. To address the inherent intra-group heterogeneity of ASCUS lesions, a group-wise differential design was implemented, comparing each cytological category (ASCUS, LSIL, HSIL) directly to the HPV-negative NILM reference (n = 51). This approach mitigated lesion-level variability by aggregating expression trends across individuals rather than relying on pairwise or hierarchical contrasts. Differentially expressed genes (DEGs) were defined by |log_2_ fold-change|≥2 and adjusted p-value*<*0.05 (Benjamini–Hochberg correction). Empirical Bayes shrinkage of *|*log_2_ fold-changes| and Wald tests were used to control for overdispersion and sample imbalance.
2.7 Functional enrichment analyses
Gene Ontology (GO) and pathway enrichment analyses were performed using WebGestalt (2023 release) and the STRING database enrichment tool. GO enrichment analysis encompassed three domains: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), with an FDR threshold of *<*0.05. For pathway analysis, KEGG, Reactome, and WikiPathways databases were queried to identify statistically overrepresented pathways (FDR < 0.05), using the full human genome as the background set.
To control for intergroup variability, enrichment scores were aggregated across all DEGs within each cytological class, allowing the identification of consensus pathways rather than individual sample–specific events. This pathway aggregation minimized the influence of outlier expression particularly in heterogeneous ASCUS lesions and highlighted reproducible biological themes such as early oncogenic signaling and immune modulation.
2.8 Hub gene identification via network centrality
The hub genes were identified through integrative network analysis using multiple centrality measures. Protein–protein interaction (PPI) networks were generated in STRING v11.5 with a minimum interaction confidence score of 0.7. All transcriptomic hub genes were mapped to their corresponding proteins to ensure that network topology reflected biologically relevant protein-level associations. Networks were filtered by node degree to retain high-confidence connections and visualized in Cytoscape (v3.10.3) and hub genes were subsequently identified using CytoHubba, a Cytoscape plugin that ranks nodes based on 12 centrality algorithms: Degree, Edge Percolated Component (EPC), Maximum Neighborhood Component (MNC), Density of MNC (DMNC), Maximum Clique Centrality (MCC), Bottleneck, Eccentricity, Closeness, Radiality, Betweenness, Clustering Coefficient, and Stress. Genes ranked in the top 20 across at least two algorithms were considered robust hub genes. To further reduce cytology-related heterogeneity and false positives, only connected network components were analyzed, excluding isolated nodes.
2.9 Driver-passenger analysis of transcriptomic alterations
To distinguish driver from passenger gene expression changes, an integrative framework was employed incorporating temporal expression trends across cytological progression (NILM → HSIL), network centrality metrics (betweenness, degree, closeness) from PPI analysis, functional enrichment context, and cross-referencing with curated cancer driver databases, including COSMIC and IntOGen. A gene was classified as a driver if (i) up-regulated across *≥*2 progressive stages, (ii) hub centrality *>*0.8 in ≥ 1 cytology, and inclusion in FDR < 0.01 pathways. Genes meeting none or only one criterion were classified as passengers.
2.10 Pathway crosstalk and network integration
To evaluate inter-pathway relationships across cytological stages, a pathway crosstalk network was constructed using Cytoscape, where each node represented a pathway and edges reflected gene overlap. Only pathways with FDR < 0.05 and ≥3 constituent genes were retained. Pairwise similarity between pathway gene sets was calculated using both Jaccard Coefficient (JC) and Overlap Coefficient (OC) based on set sizes and overlap. The JC was calculated as JC =(| A ∩ B |)/(| A ∪ B*|) and OC as OC = (|* A ∩ B |)/(min (|A|,|B|), where A and B depict gene sets of two distinct pathways [26,27]. Both coefficients range from 0 (no similarity) to 1 (complete similarity), but they emphasize different aspects of the relationship between sets. The Jaccard coefficient provides a broader, more balanced measure, while the overlap coefficient highlights the degree of containment relative to the smaller set. Pathway pairs with fewer than two shared genes were excluded. Topological analysis of the crosstalk network was conducted using the CentiScape plugin to compute nodal degree and identify highly connected “hub” pathways.
2.11 PPI network analysis
The STRING-derived PPI network for all hub genes (n = 151 nodes) was refined to include only edges with interaction confidence *≥*0.7 and node degree *>*5. The final network included 2,459 edges and was stratified into functional modules. Cluster analysis using STRING’s functional categorization highlighted three major modules: (i) chemokine response cluster, (ii) mitochondrial translation cluster, and (iii) mitotic spindle and replication cluster, corresponding to early immune evasion, mid-stage viral productivity, and high-grade dysplastic transformation, respectively.
2.12 Statistical analysis
Statistical analysis was performed using SPSS software (version 29; SPSS Inc, Chicago, IL, USA). The association between cytology categories and HPV type and E6/E7 mRNA Ct levels was evaluated using the Chi-square test. The correlation between cytology categories and E6/E7 mRNA Ct value was evaluated using Spearman’s Rank correlation. One-Way analysis of variance (ANOVA) was used to compare mean Ct value across different cytological categories, followed by Tukey’s post-hoc test for pairwise comparison. A p-value *<*0.05 was considered statistically significant.
3 Results
We conducted an integrative analysis of HPV-positive cervical Pap smear samples across four cytological stages (NILM, ASCUS, LSIL, HSIL), comparing them to HPV-negative NILM controls. Below, we report clinical specimen characteristics, E6/E7 expression patterns, differentially expressed genes (DEGs), hub gene profiles, enriched functional categories, and pathway crosstalk associated with cervical lesion progression.
3.1 Clinical specimen characterization
A total of 183 Pap Smear specimens were included in this study and analyzed, stratified into four cytological categories and one control category: Negative for Intraepithelial Lesion or Malignancy but HPV infected (HPV-positive NILM, n = 37), Atypical Squamous Cells of Undetermined Significance (ASCUS, n = 44), Low-grade Squamous Intraepithelial Lesion (LSIL, n = 35), High-grade Squamous Intraepithelial Lesion (HSIL, n = 16), and HPV-negative NILM (base control) (n = 51) (Table 1). HPV genotyping identified: high-risk types (HPV 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, and 68) and low-risk types (HPV 6, 11, 26, 53, 66, 73, and 82). Of the 132 specimens, 73 (55.3%) exhibited single HPV strain infections, while 59 (44.7%) showed multiple HPV strains infections. HPV E6 and E7 mRNA expression levels measured by reverse transcription PCR cycle threshold (Ct) values, were categorized as: strong (Ct < 25, n = 74), moderate (25 ≤ Ct < 28, n = 19), weak (28 ≤ Ct < 35, n = 30), and undetectable (Ct ≥ 35, n = 9).
Table 1: Characteristics of 132 HPV-Positive Pap Smear Specimens.
3.2 E6/E7 mRNA Ct with HPV infection and cytology categories
Statistical analyses were conducted to explore the relationship between cytological categories, HPV infection patterns (single vs. multiple), and E6/E7 mRNA expression measured by cycle threshold (Ct) values. A chi-square test showed a significant association between cytology categories and Ct level categories (X^2^ = 24.407, df = 9, p = 0.003). Spearman Rank Correlation showed a significant negative correlation between cytology categories and Ct value (ρ = −0.321, p < 0.001), suggesting that Ct values (indicating higher viral mRNA expression) tend to decrease as cytological abnormalities more severe. No significant relationship was found between cytology categories and HPV infection pattern (single vs. multiple) (X^2^ = 0.249, df = 3, p = 0.983). One-Way ANOVA indicated a significant difference in mean Ct values across the cytological categories (F = 3.17, df = 3, p = 0.027). Tukey’s post-hoc comparison revealed that the mean Ct value for LSIL cases (22.3 ± 4.2) was significantly lower than for NILM cases (25.83 ± 6.3, p = 0.019). No significant differences in mean Ct values were observed among other cytological category comparisons (ASCUS: 24.68 ± 4.7; HSIL: 22.87 ± 4.74).
3.3 Identification of DEGs
Transcriptomic analysis identified 9,194 differentially expressed genes (DEGs) across all cytological categories (HPV-positive NILM, ASCUS, LSIL, HSIL) when compared to HPV-negative NILM controls (|log_2_FC| ≥ 1, p < 0.05). Of these, 891 genes (9.7%) met the stringent criteria of |log_2_FC| ≥ 2 and p < 0.05 across the cytology groups when compared to normal controls. The distribution of significant DEGs per cytological category was as follows: NILM (498 upregulated, 12 downregulated); ASCUS (85 upregulated, 3 downregulated); LSIL (165 upregulated, 6 downregulated); and HSIL (70 upregulated, 52 downregulated). In total, 818 genes were upregulated and 73 were downregulated across all categories. DEG identification was performed using DESeq2 with false discovery rate correction for multiple testing.
3.4 Stage-specific hub gene signatures reveal molecular progression
Hub gene analysis identified 308 highly connected genes across the four cytological categories, revealing distinct functional signatures at each stage of HPV-driven progression. These hub genes represent the most influential molecular players driving stage-specific biological processes. Hub genes were subsequently identified using CytoHubba, a plugin for cytoscape.
3.4.1 HPV-positive NILM: establishing viral persistence through metabolic hijacking.
HPV-positive NILM samples (n = 91 hub genes, 56 upregulated and 35 downregulated genes) displayed a characteristic ”stealth infection” profile with two dominant themes: enhanced cellular productivity and immune suppression. The upregulated hub genes in NILM primarily functioned in protein synthesis and mitochondrial translation (37 genes, 66%; including RPL5, RPS6, MRP19, EPRS, CCT2, BYSL, NCL, PDCD11, etc.) and growth factor signaling and epithelial proliferation (EGFR, ERBB2, MYC, PLK1, PTK2, PRKCA, PRKACB). It also involved tumor suppression and DNA repair (TP53, PARP1, PALB2); hormonal regulation (ESR1), cell adhesion and cytoskeletal organization (CDH1 [28]); cellular stress response and protein homeostasis (XBP1, CUL7, GOSR2); nucleotide metabolism (TYMS), and cytokine signaling (FLT3LG). These NILM hub genes highlight their diverse roles in maintaining epithelial homeostasis.
The downregulated hub genes in NILM were involved in: immune response and inflammation (IL1β, IL8, CXCL1,CXCL2, CCL20,CCL3, PTGS2, NFKB1 [29],IKBKG, NFKBIA); immune cell signaling and activation(TYROBP, TREM1 [30]; FCER1G, CLEC77A, LCP2, HCK), apoptosis and cell death regulation(CASP5 [31], CARD16, AIM2, PLAU); cell adhesion and migration(LAPTM5, CXCR4 [32], PLEK, TNFAIP6); metabolic and stress response(HIF1A [33], FPR1, FGF13, END1); and innate immune sensors(NLRC4, FCGR1A, FCAR). These hub genes highlight their diverse roles in maintaining epithelial homeostasis and suppressing inflammation and immunity in non-dysplastic cervical epithelium (Supplementary S2 Table).
3.4.2 ASCUS: critical tipping point with early oncogenic activation.
ASCUS samples (n = 79 hub genes, 47 upregulated and 32 downregulated genes) represented a pivotal transition state, maintaining the metabolic hijacking and immunosuppresive characteristic of NILM while introducing the first oncogenic drivers that signal cellular transformation. Similar to NILM, upregulated ASCUS hub genes were involved in protein synthesis and mitochondrial translation (including RPL8, MRPL19 [34], PTCD3, GADD45GIP1 [35], etc.); growth factor signaling and epithelial proliferation (EGFR [36], CD9, PLCG1, CCND1 [37], SHH [36], NT5E, IGF1); cell adhesion and cytoskeletal organization (CDH1 [28], OCLN); and nucleotide metabolism (GMPS, CTPS2, NME3). Additionally, ASCUS upregulated hub genes functioned in hormonal and nuclear receptor signaling (ESR1, AR, PPARGC1A); DNA repair and genomic stability (PARP1, RAD51); and protein folding and glycosylation (STT3A). The downregulated hub genes in ASCUS were involved in pro-inflammatory cytokine and chemokine signaling (IL8 [38], CCL20 [39], CCL1, CXCL2, TNFAIP6 [40], IL22); chemokine receptor and GPCR signaling (CCRL2, CXCR4, C3AR1), inflammatory signaling regulation (IKBKG, TANK, MAP3K8); and inflammatory protease and effector pathways (PLAU). Additionally, downregulated hub genes in ASCUS also involved neuronal and structural proteins (MAPT, NEFH, CRYGC) and other immune and signaling pathways. Because ASCUS lesions exhibit considerable biological and transcriptomic heterogeneity, arising from mixed cellular composition, variable HPV transcriptional activity, and patient-specific immune contexts; a group-wise DESeq2 framework combined with pathway aggregation and network centrality analyses (Methods 2.6–2.8) was used to minimize intra-class variance. This approach stabilized differential expression estimates and allowed identification of consistent ASCUS-associated programs, ensuring that the hub-gene profile described above represents shared molecular features rather than stochastic inter-sample variation. These hub genes reflect ASCUS’s intermediate state, balancing growth with homeostasis to resist full dysplasia while preventing chronic inflammation and allowing subtle oncogenic changes (Supplementary S3 Table).
3.4.3 LSIL: viral productivity maximization with suppression of antimicrobial defense.
LSIL samples (n = 63 hub genes, 44 upregulated and 19 downregulated genes) exhibited a distinctive “viral factory” profile, characterized by intensified protein synthesis machinery and strategic suppression of antimicrobial defenses to optimize conditions for viral replication and assembly. Unlike the hub genes identified in NILM and ASCUS, LSIL upregulated hub genes began to involve epithelial differentiation and secretion (AGR2, DAG1) and RNA processing and centrosome function (DROSHA [41], CEP290 [42], POLR1B, ZNF512B). The downregulated hub genes started to involve neutrophil and antimicrobial defense (MPO [43], DEFA4, CEACAM8) and leukocyte adhesion and migration (SELL, AIF1 [44]). These hub genes indicate a more advanced dysplastic state than ASCUS with enhanced proliferation and immune evasion, marking progression toward malignancy (Supplementary S4 Table).
3.4.4 HSIL: genomic turbulence and epithelial dedifferentiation.
HSIL samples (n = 75 hub genes, 47 upregulated and 28 downregulated genes) demonstrated the most dramatic molecular reprogramming, with hub genes reflecting uncontrolled proliferative machinery, DNA replication stress, and complete breakdown of epithelial differentiation programs characteristic of pre-malignant transformation. HSIL hub genes exhibited distinct functional categories as compared to other cytological grades. The upregulated hub genes were involved in cell cycle regulation and proliferation (MYC [45], CDK1 [46], PLK1 [47], EZH2 [48], FOXM1 [49]); DNA replication and repair (TYMS, CTPS1, TK1, RAD51, CDC45, MCM4, MCM3, MCM5, MCM6 [50], RFC4, RFC3, FEN1, PCNA [51], POLA2, CHTF18); mitotic spindle organization and chromosome segregation (KIF11 [52], NUF2, TPX2, WDHD1, DLGAP5, CENPF, KIF2C, PTTG1 [53], KIF14, KIF15, KIF23, KIF4A, ZWILCH, MAD1L1); growth factor and tyrosine kinase signaling(KDR [54], PLCG1, ABL1), nucleotide metabolism (CMPK2); epigenetic and post-translational modification(PRMT1, MELK), RNA processing and nucleolar function(FBL, NONO, FTSJ2); and cytoskeletal and extracellular matrix remodeling (TLN2, LOX, UBE2C, KPNA2). The downregulated hub genes in HSIL were involved in cornified envelope and epidermal differences(SPRR1A, FLG [55], CDSN, LCE3D, RPTN, LCE3E, LCE2D, LCE2C, LCE3A, TGM3 [56]); keratinocyte structure protein (KRT6A, KRT1 [57], KRT6C, KRT6B, DSG1, PKP1); epidermal protease inhibition and defense(SPINK5, WFDC12, KLK14); oxidative stress (HMOX1, GPX3); lipid metabolism(MGLL, ABHD5, RBP7); and metallothionein and ion transport (MT1E, MT1G, MT1X, SLC9A7). HSIL hub genes reflect a severe dysplastic state more aggressive than LSIL, characterized by intense cell cycle dysregulation and differentiation suppression, approaching an invasive malignant state (Supplementary S5 Table).
3.5 Summary of progression
Transcriptomic analysis reveals a stepwise molecular progression from HPV infection to precancerous lesions. In NILM, despite normal cytology, we observed a coordinated upregulation of protein synthesis machinery (RPL5, RPS6) and growth signaling (EGFR, MYC) alongside strategic downregulation of inflammatory pathways (IL1B, IL8, CXCL1/2), establishing viral persistence while maintaining epithelial homeostasis. These dual strategies-upregulating host translational and proliferative machinery while simultaneously suppressing immune signaling-represent HPV’s initial foothold in establishing persistent infection despite normal cytology.
ASCUS maintains this translational emphasis while introducing early oncogenic drivers (CCND1, SHH) and continued immune suppression. These oncogenic factors signal the beginning of cellular transformation without full dysplastic commitment. As lesions progress to LSIL, hub genes shift toward epithelial differentiation alteration (AGR2, DAG1) and intensified viral productivity, with suppression of antimicrobial defenses (MPO, DEFA4). This stage represents active viral replication and assembly, reflected in the robust translational signature and emerging viral escape from epithelial defense mechanisms.
The transition to HSIL marks the most dramatic molecular reprogramming, characterized by intense cell cycle hyperactivation (CDK1, PLK1), DNA replication licensing (MCM3–6, RFC3/4), and chromosomal instability (KIF family, CENPF), concurrent with significant downregulation of epithelial differentiation markers (SPRR1A, KRT family) and barrier function genes. This represents HPV’s transition from productive infection to transformative integration, with viral oncoproteins driving uncontrolled proliferation, genomic instability, and loss of differentiation.
Preliminary trends, further supported by enrichment analysis in subsequent sections, revealed a clear progression from homeostatic maintenance in NILM to increasingly dysregulated cell cycle and suppressed differentiation in HSIL (p < 0.01). These findings demonstrate HPV’s molecular evolution from immune evasion to host genome subversion, with distinct hub gene signatures marking critical transitions in dysplastic severity.
3.6 GO enrichment analysis of hub genes
Gene ontology (GO) analysis was performed across cytological groups to identify the functional roles of hub genes associated with cervical lesion progression. This analysis included three primary aspects: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The functional characterization of genes associated with four cytological stages was performed using WebGestalt. An over-representation analysis was performed, selecting “Biological Process” as the GO functional category, with gene symbols as identifiers (ID type) and entire genome as the enrichment background.
3.6.1 Biological Process (BP) analysis reveals functional reprogramming.
Gene Ontology analysis revealed a systematic shift in cellular priorities across disease progression. Early-stage lesions (NILM/ASCUS) prioritized metabolic processes (51 and 42 genes respectively), reflecting the virus’s hijacking of host biosynthetic machinery. However, HSIL demonstrated a fundamental reprogramming toward regulatory and stress-response functions, with significant enrichment in biological regulation (37 genes, p < 0.001), stimulus response (32 genes), and cellular reorganization (35 genes) (Fig 1A). This transition from metabolic support to dysregulated control mechanisms marks the shift from viral persistence to malignant transformation. Downregulated processes showed consistent immune suppression across all stages, with response to stimuli (34 genes) and cell communication (33 genes) prominently suppressed in NILM, establishing the immunotolerant foundation for progression (Fig 1B).
GO Biological Process Enrichment of Differentially Expressed Hub Genes Across Cervical Lesion Stages.Bar plots illustrating Gene Ontology (GO) biological process (BP) enrichment of upregulated (Fig 1A), downregulated (Fig 1B) and hub genes across four cervical cytological stages: NILM (Negative for Intraepithelial Lesion or Malignancy), ASCUS (Atypical Squamous Cells of Undetermined Significance), LSIL (Low-grade Squamous Intraepithelial Lesion), and HSIL (High-grade Squamous Intraepithelial Lesion). Enrichment analyses were performed using WebGestalt (2023 release), applying a false discovery rate (FDR) threshold of <0.05. Cytological groups are color-coded as follows: blue (NILM), green (ASCUS), red (LSIL), and orange (HSIL). Bar length indicates the number of genes associated with each enriched term, and the x-axis shows –log10 transformed FDR values indicating statistical significance.
3.6.2 Cellular component (CC) analysis reveals subcellular reorganization.
Cellular component enrichment demonstrated a progressive shift from cytoplasmic metabolic compartments in early lesions to nuclear and chromosomal structures in advanced disease. Early-stage lesions (NILM/ASCUS) showed dominant enrichment in protein-containing complexes (49 and 35 genes respectively) and membrane-enclosed lumens, reflecting enhanced biosynthetic activity and organellar function. HSIL exhibited a dramatic transition toward nuclear compartments (40 genes), cytoskeletal apparatus (19 genes), and chromosome-associated structures (25 genes), indicating the cellular reorganization accompanying malignant transformation (Fig 1A). Downregulated components consistently involved membrane systems and extracellular structures, suggesting progressive breakdown of epithelial barrier function (Fig 1B).
3.6.3 Molecular function (MF) analysis reveals enzymatic reprogramming.
Molecular function analysis revealed systematic changes in cellular enzymatic priorities throughout disease progression. While protein binding remained consistently enriched across all stages (41–45 genes), early lesions (NILM/ASCUS) emphasized nucleic acid binding (36 and 28 genes respectively) and structural molecular activities, supporting the enhanced transcriptional and translational programs. HSIL demonstrated a striking shift toward catalytic functions, with elevated hydrolase activity (15 genes), transferase activity (12 genes), and ion binding (27 genes), reflecting the intensive enzymatic remodeling required for uncontrolled proliferation and chromatin reorganization. This functional reprogramming marks the transition from metabolic support to dysregulated cellular machinery (Fig 1B).
3.7 Functional progression across cytological categories
Gene ontology analysis demonstrates distinct functional profiles from NILM to HSIL across multiple GO domains. In Biological Process terms, upregulated genes shift from metabolic activity to regulatory functions. Metabolic processes are enriched in earlier stages: NILM (24.4%), ASCUS (24.1%), and LSIL (25.6%) but decrease in HSIL (15.4%). Conversely, genes involved in biological regulation increase from NILM (11.96%) to HSIL (15.81%), becoming the most enriched category in HSIL. Additionally, we observed stress response elevation and tissue architecture disruption in HSIL: response to stimulus peaks in HSIL (13.7% vs ≤ 9.8% in other stages), while cellular component organization doubles in HSIL(15% vs ≤ 8.6% in other stages). Meanwhile, cell communication is significantly reduced in ASCUS (6.3%) and LSIL (7.5%) compared to NILM (8.6%), showing mid-grade signaling suppression. Downregulated genes also exhibit distinct stage-specific patterns: HSIL shows severe suppression of cell communication (4% vs 13–15% in other stages) and cell population proliferation (2% vs 5% in NILM/LSIL), indicating progressive loss of growth control mechanisms. Furthermore, HSIL doubled disruption of developmental processes (18.2%) compared to early stages (*≤*8.7%).
In the Cellular Component terms, upregulated genes in HSIL reflect marked architectural changes. Nuclear components are significantly elevated in HSIL (16.9%) compared to NILM (8.6%), and chromosomal elements rise dramatically to 10.5%, representing a five-fold increase over other stages (*≤*2.8%). In contrast, components related to mitochondria and protein synthesis machinery, such as ribosomes and envelopes are notably reduced or absent in HSIL (2.5% vs. 7.9–9.1% in earlier stages). Cytoskeletal reorganization is distinctly enriched in HSIL (8.0% vs. ≤ 2.2% in other grades), while membrane and protein complexes remain relatively stable across all stages. Among downregulated genes, membrane components remain consistent across stages. However, HSIL shows pronounced suppression of cytosolic (15.1%) and nuclear (11.8%) components. LSIL, on the other hand, is marked by disruption in the extracellular space (17.7%) and vesicular transport systems, including the endomembrane system and vesicles (16.1% each). Notably, protein complexes appear to be selectively retained in HSIL downregulation (2.2%), contrasting with higher suppression rates in earlier stages (8–10%).
In Molecular Function terms, upregulated genes in HSIL reveal a distinct pattern of functional reprogramming. Ion binding becomes the dominant activity in HSIL (16.6% vs. ≤ 12.9% in earlier stages), while structural molecular activity is nearly eliminated (0.6% compared to 12–15% in other stages). Nucleic acid binding is most prominent in normal cytology (22.9%) but decreases to 14.7% in HSIL. Catalytic activity shows stage-dependent changes, with hydrolase activity rising fourfold from NILM (2.5%) to HSIL (9.2%), and chromatin binding increasing during dysplastic progression. Protein binding remains stable across all stages (26–28%), suggesting preservation of core molecular interactions. Downregulated genes show pronounced suppression of structural molecule activity (16.1%), and a near-complete loss of molecular transducer activity (1.8% vs. 5–13% in other stages) in HSIL. Protein binding downregulation increases progressively, peaking in LSIL (50%), while catalytic hydrolase activity is entirely absent at this stage. ASCUS displays a specific decline in molecular adaptor functions (7.1% downregulation), which are completely lost in HSIL.
This GO term progression provides a functional framework for understanding the systematic reprogramming of cellular machinery during HPV-driven neoplastic transformation.
3.8 Pathway enrichment analysis of hub genes
Pathway enrichment analysis was performed to highlight the key molecular mechanisms and their implications for cervical lesion progression, focus on statistical significance (FDR), biological relevance, and gene involvement. STRING enrichment tools, Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and WikiPathways databases were used. To ensure that the enrichment analysis carried out across these databases provides comprehensive pathway coverage, we considered pathways with a false discovery rate (FDR), p < 0.05. The analysis revealed a distinct stage-specific pathway progression across cytological categories. In NILM group, the upregulated hub genes were most significantly enriched in mitochondrial translation (elongation, initiation, termination), metabolism of proteins, ribosome/translation (Supplementary S6.1 Table), while downregulated genes were enriched in immune system, innate immunity, NOD-like receptor signaling, chemokine signaling, and IL-17 signaling pathways(Supplementary S6.2 Table), reflecting a balance between active proliferation and suppressed inflammation that maintains epithelial homeostasis.
ASCUS showed a shift towards mitochondrial translation and cancer-related pathways (e.g.,pathways in cancer, proteoglycans in cancer, breast cancer, prostate cancer, bladder cancer) in upregulated genes, with notable EGFR tyrosine kinase inhibitor resistance (Supplementary S7.1 Table), while downregulated genes remained enriched in inflammatory signaling (IL-17, TNF, NF-κ B) and cytokine-cytokine receptor interactions, indicating cells under proliferative stress with emerging immune evasion (Supplementary S7.2 Table).
LSIL maintained upregulation in mitochondrial translation and protein metabolism while adding hormonal and receptor related pathways (mammary gland development and nuclear receptors) (Supplementary S8.1 Table), with downregulation in cytokines/chemokines pathways representing a transition toward dysplasia with metabolic reprogramming in a tolerogenic microenvironment (Supplementary S8.2 Table).
HSIL demonstrated significant upregulation in cell cycle/DNA replication, DNA repair/replication stress response, and oncogenic signaling/transcriptional dysregulation pathways (Supplementary S9.1 Table), with downregulation in cornified envelope formation, oxidative stress response, and mineral homeostasis pathways, reflecting uncontrolled proliferation, genomic instability, and epithelial barrier disruption (Supplementary S9.2 Table). These findings illustrate the progression from metabolic equilibrium in NILM to immune evasion in ASCUS, metabolic/hormonal rewiring in LSIL, and ultimately cell-cycle hyperactivation, replication stress, and chromosomal instability characteristic of high-grade dysplasia and epithelial disintegration in HSIL, providing a comprehensive molecular framework for understanding cervical carcinogenesis.
3.9 HPV-related pathways by cytology
After excluding pathways unrelated to HPV, the most biologically relevant pathways associated with HPV infection and cervical lesion progression were presented as follows (Supplementary S10 Table):
3.9.1 NILM pathway analysis.
Analysis of HPV-positive NILM samples revealed enrichment of multiple cellular and immune pathways. Enriched downregulated pathways included immune-related processes such as Innate Immune System (KEGG pathway ID(hsa)168249, FDR = 1.25e-13), Adaptive Immune System (hsa1280218, FDR = 0.0147), Toll-like receptor signaling (hsa04620, FDR = 7.47e-07), NOD-like receptor signaling (hsa04621, FDR = 5.65e-12), and RIG-I-like receptor signaling (hsa04622, FDR = 0.00012). The presence of these immune signatures potentially reflects baseline immune surveillance in morphologically normal cervical epithelium exposed to various microbes. Additionally, the Apoptosis pathway showed significant enrichment in NILM (WikiPathways identifier(WP254, FDR = 0.0059). This finding is noteworthy as functional apoptotic mechanisms are essential for eliminating damaged cells in normal tissue, whereas HPV inhibits apoptosis through E6-mediated targeting of p53 and various pro-apoptotic factors. The presence of intact apoptotic machinery in NILM, contrasted with its subsequent downregulation in progressive lesions, highlights another critical mechanism through which HPV promotes carcinogenesis by enabling damaged cells to evade programmed cell death (Supplementary S10 Table).
Cell adhesion and communication regulated pathways were also significantly enriched in NILM samples, including Adherens Junction (hsa04520, FDR = 0.0164), Gap Junction (hsa04540, FDR = 0.0247), and Focal Adhesion pathways (hsa04151, FDR = 0.00029). These findings establish a reference baseline for epithelial integrity, as normal epithelium maintains intact E-cadherin-mediated adherens junctions and connexin gap junctions. In subsequent HPV-induced lesion stages, these adhesion structures are frequently compromised, as HPV E6 targets PDZ-domain cell polarity/adhesion proteins (including hDlg, Scribble, and Magi), leading to junction disruption. The presence of these intact adhesion pathways in NILM underscores that their progressive downregulation represents a hallmark of neoplastic transformation (Supplementary S10 Table).
Additionally, several signaling cascades were enriched in NILM, including PI3K-Akt signaling (WP4172, FDR = 0.0024), MAPK signaling (hsa04010, FDR = 0.00068), Wnt(hsa04310, FDR = 0.016), and Hippo signaling pathways (WP4540, FDR = 0.0065). These represent fundamental cell regulatory mechanisms that maintain normal tissue homeostasis but become progressively dysregulated during HPV-driven transformation (Supplementary S10 Table).
3.9.2 ASCUS pathway analysis.
In ASCUS samples, pathway enrichment analysis revealed a notable pattern of immune signaling pathways. Significantly enriched downregulated pathways included Chemokine signaling (hsa04062, FDR = 0.00019), Cytokine–cytokine receptor interactions (hsa04060, FDR = 0.00068), IL-17 (hsa04657, FDR = 0.0014) and IL-18 signaling (WP4754, FDR = 0.0119), Toll-like and NOD-like receptor Signaling (hsa04620, FDR = 0.0172), and NF-κB signaling (hsa04064, FDR = 0.0016). This pattern reflects the ongoing interaction between host defense mechanisms and viral immune evasion strategies at this early stage of abnormality (Supplementary S10 Table).
While our hub gene analysis identified downregulation of specific immune mediators in ASCUS, pathway analysis reveals that the broader immune signaling architecture remains intact and detectable, suggesting complex regulation.
In early HPV infection, keratinocytes produce pro-inflammatory chemokines/cytokines (e.g., IL-8, CCL20) to recruit immune cells, creating a local inflammatory milieu to combat the virus [58]. Indeed, in ASCUS we observe chemokine pathways with genes like CXCL8 (IL-8) and CCL20, which are known to be induced in HPV-infected epithelial cells [58]. These factors recruit neutrophils and Th17 cells to the site of infection. The concurrent presence of IL-10 signaling (an immuno suppresion cytokine pathway) suggests HPV may simultaneously trigger regulatory feedback to dampen inflammation. HPV employs multiple immune evasion strategies, for example, high-risk HPV E6/E7 oncoproteins can downregulate CCL20 to evade immune detection [59]. This balance between immune activation and suppression aligns with our hub gene findings, where key immune mediators showed downregulation despite pathway-level detection.
Additionally, ASCUS samples showed upregulated enrichment of the ErbB signaling pathway (WP673, FDR = 0.0445), which involves EGFR and related growth factor signals. This finding is consistent with our hub gene analysis that identified upregulation of growth factor signaling components in ASCUS. Early HPV infection can enhance growth factor signaling; for instance, the HPV16 E5 protein amplifies EGFR signaling, promoting cell proliferation. The ”Proteoglycans in cancer” pathway was also enriched (hsa05205, FDR = 0.00085), reflecting changes in the extracellular matrix and cell signaling in the tumor microenvironment. Even at the ASCUS stage, HPV-infected cells begin altering the local environment (e.g., via proteoglycans, integrins) as a prelude to transformation.
These pathway findings complement our hub gene analysis, demonstrating that while specific immune hub genes are strategically downregulated by HPV, the immune signaling architecture remains detectable at the pathway level, highlighting the complex molecular dynamics in early HPV-induced abnormalities.
3.9.3 LSIL pathway analysis.
In LSIL samples, pathway enrichment analysis demonstrated a continued predominance of immune and inflammatory signaling pathways, similar to the pattern observed in ASCUS. Significantly enriched downregulated pathways included Chemokine (hsa04062, FDR = 7.14e-12) and Cytokine signaling (hsa04060, FDR = 6.00e-12), TNF signaling(hsa04668, FDR = 0.0038), Toll-like (hsa04620, FDR = 0.003) and NOD-like receptor pathways (hsa04621, FDR = 0.00092), IL-17 (hsa04657, FDR = 0.00012)/IL-18 signaling (WP4754, FDR = 8.74e-07), as well as broader “Immune System” and “Cytokines and inflammatory response” pathways (Supplementary S10 Table). These pathway findings complement our hub gene analysis, which identified downregulation of specific antimicrobial defense genes (MPO, DEFA4) in LSIL, suggesting that while certain immune effectors are suppressed, broader immune signaling networks remain activated.
The persistence of these immune pathways indicates ongoing host immune engagement even in LSIL lesions. Pro-inflammatory cytokines like TNFα and IL-17 identified in these pathways contribute to immune surveillance and tissue remodeling in lesions. HPV’s continued presence maintains activation of these pathways, a pro-inflammatory microenvironment commonly accompanies HPV lesions. However, the virus successfully evades complete clearance, as evidenced by persistent infection. The continued enrichment of IL-10 signaling in LSIL (also observed in ASCUS) suggests HPV-induced immune suppressive mechanisms, as IL-10 promotes an anti-inflammatory microenvironment favorable to HPV persistence [58].
Notably, LSIL samples retained enrichment of NF-κB signaling (hsa04064, FDR = 0.00015) and Toll-like receptor signaling pathways (hsa04620, FDR = 0.003), both crucial for innate immune responses. Keratinocytes express Toll-like receptors (particularly TLR9) that detect HPV DNA and activate NF-κB and type I interferon responses [60]. HPV can modulate these signals through mechanisms such as repressing TLR9 expression and downstream IL-6/IL-8 production to evade effective immune responses [60]. The enrichment of these pathways underscores their importance in HPV containment. The pathway “Viral protein interaction with cytokine and cytokine receptor”(hsa04061, FDR = 3.89e-16) was also enriched, reflecting the ongoing molecular interaction between viral immune evasion strategies and host cytokine responses.
These findings align with our hub gene analysis showing that LSIL represents a stage characterized by enhanced viral productivity while simultaneously suppressing specific antimicrobial defenses. The persistence of immune pathway signatures despite viral evasion mechanisms suggests a dynamic equilibrium that allows for viral replication while preventing complete immune collapse, creating conditions conducive to lesion persistence.
3.9.4 HSIL pathway analysis.
HSIL represents a precancerous state characterized by high-risk HPV driving significant cellular transformation. Upregulated pathway enrichment analysis revealed a dramatic shift in the molecular landscape of HSIL samples, dominated by cell cycle dysregulation (hsa04110, FDR = 1.71e-12), DNA replication (WP466, FDR = 7.99e-13), and DNA damage repair pathways (WP4946, FDR = 0.001), hallmark processes of HPV-driven oncogenesis (Supplementary S10 Table).
These pathway findings strongly align with our hub gene analysis, which identified upregulation of cell cycle regulators (CDK1, PLK1, FOXM1), DNA replication factors (MCM family, RFC3/4), and mitotic mediators (KIF family, CENPF) as key features of HSIL. High-risk HPV16/18 E6 and E7 oncoproteins, persistently expressed in HSIL, drive uncontrolled cell proliferation and genomic instability [61]. The significantly enriched pathways in HSIL encompass multiple aspects of the Cell Cycle, including Cell Cycle checkpoints (hsa-69620, FDR = 5.35e-14), Mitotic cell cycle (hsa-69278, FDR = 1.88e-25), and APC/C complex regulation for mitosis (hsa-176409, FDR = 0.019). Additionally, DNA damage response pathways were prominently enriched, including ATR activation (hsa-176187, FDR = 2.22e-09), G2/M DNA damage checkpoint (hsa-69473, FDR = 0.0202), DNA double-strand break repair (WP3959, FDR = 0.0154), Homologous recombination repair (hsa5693538, FDR = 3.17e-05), Mismatch repair (hsa03430, FDR = 0.0035), Translesion DNA synthesis (hsa5655862, FDR = 0.00059), and Telomere maintenance pathways (hsa-174417, FDR = 8.28e-08; hsa5656121, FDR = 0.00059; hsa110320, FDR = 0.00075)(Supplementary S10 Table).
This comprehensive dysregulation of cell-cycle checkpoints and repair processes reflects the fundamental mechanisms of HPV oncogenesis: HPV E7 inactivates Rb, forcing cells into S-phase, while E6 inactivates p53, disabling the G1/S checkpoint and DNA damage checkpoints. The result is uncontrolled DNA replication accompanied by accumulating genomic damage that evades proper repair. HPV oncoproteins impair multiple DNA repair pathways (particularly p53- and pRb-regulated processes), promoting cell cycle progression at the expense of genomic integrity [62]. The enrichment of repair pathways in HSIL suggests cellular attempts at DNA damage management that are simultaneously undermined by viral oncoproteins, leading to progressive mutations and instability.
Activation of ATR in response to replication stress emerged as a key enriched pathway in HSIL (hsa176187, FDR = 2.22e-09), HPV-induced unscheduled replication triggers replication stress that activates ATR. Notably, HPV co-opts the ATR/ATM DNA damage response to facilitate its own viral DNA replication in differentiating cells [63]. However, in the context of neoplastic progression, chronic activation of ATR/ATM signifies persistent DNA damage. Related HSIL-enriched pathways such as DNA damage response (WP707, FDR=0.0021), G2/M checkpoint (HSA-69473, FDR=0.0202), and APC/C-mediated degradation of cell cycle proteins (HSA-174048, FDR=0.0351) further demonstrate how E6/E7 drive cells through division despite accumulating damage. Research confirms that E7 causes replication stress and activates ATM/ATR pathways [64], while the combination of E6/E7 constitutively activating the cell cycle while disabling p53 leads to progressive genomic instability [62].
HSIL samples also showed enrichment of select immune-related pathways, though fewer than observed in lower-grade lesions. Notably, the Adaptive Immune System and MHC class II antigen presentation pathways remained enriched (hsa2132295, FDR = 0.00054). At the HSIL stage, adaptive immune responses (including T cell recognition of HPV antigens) should be active. However, many HSIL lesions develop an immune suppressive microenvironment (characterized by PD-L1 upregulation and regulatory T-cell accumulation), facilitating lesion progression despite immune surveillance. The persistence of antigen presentation pathway enrichment indicates ongoing presentation of HPV antigens via MHC II to helper T-cells, though HPV employs various mechanisms to evade immune clearance.
These pathway findings reveal a clear progression across cervical lesion stages: pathways in NILM were primarily associated with “immune signaling” and “viral oncogenesis”; pathways in ASCUS were mainly categorized as “immune signaling” or “viral oncogenesis”; pathways in LSIL were predominantly classified as “immune signaling”; while pathways in HSIL were predominantly categorized as “viral oncogenesis.” This molecular progression aligns precisely with our hub gene analysis, demonstrating the transition from immune evasion and viral persistence in early stages to genomic instability and cellular transformation in HSIL.
3.10 Pathway crosstalk analysis
Pathway crosstalk analysis was performed to identify shared molecular signatures and interaction networks between the cytological categories. We constructed a network graph where nodes represent unique pathways and edges represent gene overlap between pathways. Pathways with FDR adjusted p-value >0.05 or those represented by fewer than 3 genes were excluded to maintain analytical rigor. Furthermore, to calculate the overlap between pathways, the Jaccard Coefficient (JC) and Overlap Coefficient (OC) was calculated. Specifically, the JC was calculated as JC =(| A ∩ B |)/(| A ∪ B*|) and OC as OC = (|* A ∩ B |)/(min (|A|,|B|), where A and B depict gene sets of two distinct pathways. Cytoscape was then used to visualize the interrelationships among these pathways. Network analysis showed a fair interconnection among the cytological categories representing 81 pathways shared across *>*2 cytological categories (1513 edges) and none of them are isolated (Fig 2).
Pathway Crosstalk Network Across Cervical Cytological Categories (Full Network).This network diagram illustrates the inter-pathway crosstalk among statistically enriched biological pathways identified across four cervical cytological categories: NILM (blue), ASCUS (green), LSIL (red), and HSIL (orange). Each node represents a unique pathway enriched in one or more cytological stages (FDR < 0.05, ≥ 3 constituent genes) and edges represent gene overlap between pathways. The overlap between pathways is quantified using the Jaccard Coefficient (JC) and Overlap Coefficient (OC), with a minimum of two overlapping genes required for inclusion. Pathways exhibiting the highest centrality are visualized as hub nodes and reflect biologically critical signaling axes during neoplastic progression. Notably, NILM pathways clustered toward immune regulation and epithelial adhesion, ASCUS and LSIL retained significant overlap with NILM (OC > 0.9), while HSIL pathways formed a distinct cluster dominated by cell cycle, replication, and chromatin remodeling processes. Network topology was visualized in Cytoscape and pathway identities are annotated according to KEGG, Reactome, and WikiPathways databases.
To identify the most interconnected and influential pathways, we performed nodal degree analysis using the CentiScape plugin to identify key pathway nodes with nodal degree *>*5, which resulted in 668 edges depicted in Fig 3. NILM showed the largest distinct enriched pathways (n = 234), followed by ASCUS (n = 46), LSIL(n = 41), and HSIL(n = 89). JC and OC demonstrated high overlap between ASCUS and LSIL with NILM (OC = 0.978, 0.927, respectively), indicating that ASCUS shared most of its pathway profile with NILM (45 out of 46 ASCUS pathways were found in NILM), and LSIL sharing a large portion of its pathways with NILM (38 out of 41 LSIL pathways were found in NILM). Additionally, a moderate similarity (OC = 0.512, 21 out of 41 LSIL pathways were found in ASCUS) was observed between ASCUS and LSIL. In contrast, HSIL showed minimal overlaps with other cytological categories with NILM (OC = 0.067), ASCUS (OC = 0.043), and LSIL (OC = 0), indicating the set of pathways in HSIL was largely unique to the high-grade lesion category.
Pathway Crosstalk Network Across Cervical Cytological Categories (Node degree ≥5).This network visualization depicts the subnetwork with nodal degree ≥ 5 representing crosstalk among significantly enriched molecular pathways (FDR < 0.05; ≥ 3 genes) across four cytological stages of HPV-positive cervical lesions: NILM (blue), ASCUS (green), LSIL (red), and HSIL (orange). Each node represents a unique biological pathway identified through enrichment analysis using KEGG, Reactome, and WikiPathways databases. Node color indicates the cytological stage(s) in which the pathway was enriched. Edges between nodes reflect shared gene content, with a minimum threshold of two overlapping genes for edge retention. The network illustrates how pathway architecture evolves with disease progression: early-stage NILM and ASCUS pathways cluster around immune signaling and homeostatic regulation, while LSIL pathways show transitional overlap with metabolic and growth factor signaling. In contrast, HSIL forms a distinct, tightly connected module centered on cell cycle regulation, DNA replication, mitotic spindle assembly, and chromatin remodeling. High-connectivity nodes such as “Cell Cycle,” “DNA Replication,” and “ATR Activation in Response to Replication Stress” emerge as central hubs in HSIL, reflecting oncogenic pathway consolidation. The increasing modularity and divergence of HSIL-associated pathways highlight a molecular shift from viral persistence to host genome destabilization and dysplastic transformation.
ASCUS shares the immuno suppression and hyper-translational state with NILM and LSIL, indicating these are fundamental to HPV infection. What differentiates ASCUS is the incipient activation of oncogenic pathways: the first signs of the host cell responding to viral oncoproteins with changes characteristic of cancer development (growth signals, cell cycle dysregulation). These changes signal a molecular infection point, where immune evasion gives way to nascent oncogenic transformation.
The LSIL molecular signature shares the two main themes of early HPV infection (immune evasion and enhanced protein synthesis) with NILM and ASCUS. What becomes most pronounced in LSIL is the protein synthesis aspect, corresponding to actual virion production. LSIL can thus be distinguished by extremely high levels of ribosomal and mitochondrial protein genes. In contrast, LSIL does not yet share the aggressive cell-cycle signature of HSIL. LSIL is considered as a state where the virus is reproducing efficiently, but the host cell genome is largely intact (no integration in most cases) and the cell’s malignant transformation is not fully realized. Thus, LSIL is often regarded as a transient infection state and can regress if the immune system eventually responds. The HSIL signature is quite specific and distinct from NILM/ASCUS/LSIL. The massive upregulation of cell cycle and replication pathways is unique to high-grade lesions. Early-stage signatures (immune down, translation up) are less conspicuous in HSIL, largely because these aspects are already established (immune pathways long suppressed, and protein synthesis is generally high in all dividing cells). In practice, markers like MCMs, TOP2A, KI-67, and p16 are used to identify HSIL cells, all reflecting the uncontrolled cell cycle that our pathway analysis also captures [65]. HSIL does not share many of the unique “productive infection” features seen in LSIL. For example, HSIL exhibited diminished viral late gene expression (the virus often integrates and stops making L1/L2). Therefore, LSIL and HSIL have almost opposite dominant signatures; LSIL is defined by viral gene expression, HSIL by host cell gene dysregulation.
3.11 PPI network analysis
The PPI was constructed using STRING database and visualized in Cytoscape containing 151 nodes and 2459 interconnected edges. The nodal degree *>*5 was maintained for PPI network analysis depicted in Fig 4. The STRING analysis categorized a list of genes into 3 distinct categories which comprise cellular response to chemokine, mitochondrial cluster, and mitotic cell cycle processes (Fig 5).
Global Protein–Protein Interaction (PPI) Network of Hub Genes Across Cytological Categories.This figure illustrates the comprehensive PPI network of all hub genes identified across the four cytological categories of HPV-positive cervical specimens: NILM, ASCUS, LSIL, and HSIL. The network was constructed using the STRING database (v11.5), with a high-confidence interaction score threshold of ≥ 0.7 to ensure biologically meaningful associations. Only nodes (genes) with a degree > 5 were retained for visualization, resulting in a network of 151 nodes and 2,459 edges. Each node represents a unique hub gene, and edges indicate protein–protein interactions based on experimental data, curated databases, and predictive algorithms. The spatial organization of the network was generated using force-directed layout in Cytoscape (v3.10.3), allowing clusters and highly interconnected regions to emerge visually without imposed functional categorization. This unlabeled network provides a global overview of hub gene connectivity and highlights the dense interaction patterns that underpin the transcriptomic reprogramming during HPV-driven cervical lesion progression. The visualization serves as a reference for understanding network topology prior to functional modular classification.
Protein–Protein Interaction (PPI) Network of Hub Genes with Functional Cluster Labeling.This figure displays the high-confidence PPI network constructed from all identified hub genes across the four cytological stages of HPV-positive cervical lesions. Network construction was performed using the STRING database (v11.5), with a minimum interaction confidence score threshold of 0.7 to retain high-confidence associations. The resulting network includes 151 nodes (hub genes) and 2,459 edges, visualized in Cytoscape (v3.10.3). Only genes with a node degree >5 were retained to minimize noise and highlight core interactors. Nodes represent individual hub genes, and edges denote experimentally validated or predicted protein–protein associations. Functional modules were identified through STRING functional enrichment clustering and labeled based on biological relevance. Three major functional clusters are highlighted: 1. Mitochondrial Translation Cluster: Encompassing genes involved in ribosomal biogenesis and protein synthesis, reflecting enhanced metabolic and translational activity characteristic of early lesion stages (e.g., NILM and ASCUS). 2. Chemokine Response Cluster: Including immune modulators and signaling components involved in cytokine–cytokine receptor interactions and inflammatory signaling pathways, largely downregulated in early lesions but indicative of host immune surveillance. 3. Mitotic and DNA Replication Cluster: Comprising cell cycle regulators, DNA replication licensing factors, and mitotic spindle components, predominantly enriched and upregulated in high-grade lesions (HSIL), representing proliferative dysregulation and genomic instability. This cluster-labeled network underscores the temporal transition from immune evasion and metabolic support (early lesions) to proliferative transformation (HSIL), providing insight into molecular reprogramming along the cervical neoplastic continuum.
4 Discussion
Although previous research had enhanced our understanding of transcriptome profiles in cervical cancer, few studies have systematically tracked the activation of key molecules and pathways across the full cytological spectrum (NILM, ASCUS, LSIL, HSIL). This gap hinders the identification of targetable molecular processes for cancer prevention and therapy. A deeper understanding of the molecular events driving cervical cancer initiation and progression is essential to refine screening protocols, enhance predictive diagnostics, and guide the development of effective treatment strategies.
In this study, we conducted a comprehensive omics-analysis of 132 HPV-positive Pap smear samples by integrating HPV genotyping, E6/E7 mRNA expression, and transcriptomic profiling. By using HPV-negative NILM as our transcriptomic baseline, we were able to distinguish infection-associated alterations from normal homeostatic gene expression, revealing the early molecular foot print of HPV, even in cytologically normal samples. Our findings revealed that HPV infection induces distinct molecular changes in cervical epithelial cells that vary depending on the cytological stage.
4.1 Molecular prodrome in HPV-positive NILM: gateway to early prevention
Our transcriptomic analysis of HPV-positive cervical samples reveals that even cytologically normal (NILM) specimens harbor significant molecular alterations that precede morphological abnormalities. This ”molecular prodrome” offers a unique opportunity for early risk prediction and preventive intervention.
In HPV-positive NILM samples, we identified significant molecular alterations. A persistent upregulation of growth factor receptor signaling (EGFR, ERBB2, PTK2, PRKCA; log_2_FC > 1.5, p < 0.001) with high betweenness centrality (*>*0.85) was observed. These genes promote cell survival and proliferation and may be influenced by the activity of HPV oncoproteins. Although our study did not directly assess E5 expression, prior studies have shown that HPV16 E5 enhances EGFR signaling by promoting receptor recycling and sustained activation [66]. In addition, HPV E6 and E7 broadly disrupt multiple growth signaling cascades and reduce the expression of cell-surface adhesion molecules to evade immune surveillance [66].
In parallel, the upregulation of ribosomal and protein translation pathways (RPL5, RPL8, RPS6, MRPL14, MRPS11, MRPS2, MRPL48, GADD45G1P1) indicate a hypermetabolic state, consistent with reports of HPV-driven metabolic reprogramming via E6-mediated activation of the mTOR signaling pathway [66,67].
Concurrently, strategic downregulation of immune mediators (IL1B, IL8, CXCL1/2, CCL20, NFKB1, IKBKG) contribute to an immune suppressive microenvironment conducive to viral persistence. Prior studies have shown that HPV E6/E7 can inhibit NF-kB signaling, suppressing antiviral cytokine expression [65], and downregulate CCL20 through NF-κB pathway suppression [68].
These findings reveal a “stealth infection” profile of the cervical epithelium during HPV persistence, characterized by proliferative signaling, metabolic activation, and immune evasion.
4.2 Transcriptome-defined stage-specific signatures: beyond current paradigms
ASCUS represents a crucial early stage in the development of cervical cancer. At this point, ASCUS maintains the HPV-driven immune evasion and metabolic activation seen in NILM. Simultaneously, it introduces oncogenic signaling pathways, such as “Pathways in Cancer” and “Proteoglycan in Cancer”. Key genes involved include CCND1, SHH, IGF1, CDH1, and EGFR. Early overexpression of CCND1 is a known event in cervical dysplasia, contributing to unscheduled cell cycle entry [65]. This, alongside persistent immune suppression and enhanced protein synthesis establishes a molecular environment primed for progression. Notably, ESR1 appears in the ASCUS “Pathway in Cancer” gene list, suggesting that estrogen signaling may be active in ASCUS lesions, potentially promoting growth of infected cells.
LSIL shares two key features of early HPV infection: immune evasion and enhanced protein synthesis with NILM and ASCUS, but it also exhibits a distinct “viral productivity” signature. This is characterized by heightened translational activity and downregulation of antimicrobial defenses (MPO, DEFA4, CEACAM8), creating an environment optimized for viral replication and assembly. The robust translational signature (MRPS9, MRPL19, RPL8, RPS6, GADD45GIP1) reflects the virus’s metabolic hijacking of host machinery for virion production, though LSIL has not yet acquired the aggressive cell-cycle dysregulation seen in HSIL. At this stage, the virus reproduces efficiently, but the host genome remains largely intact, and malignant transformation is not yet fully realized. This may explain why LSIL is often a transient infection state that can regress if the immune system mounts an effective response. This stage-specific signature offers potential targets for antiviral therapies aimed at disrupting viral assembly before genomic integration occurs.
HSIL, by contrast, displays a distinct molecular profile from NILM, ASCUS and LSIL. It is characterized by dramatic molecular reprogramming dominated by cell cycle hyperactivation (CDK1, cyclin B1, CENPF, KIF2C, PLK1, FOXM1) and DNA replication licensing (MCM family, CDC45, PCNA, POLA2, FEN1, RFC3/4), along with epithelial dedifferentiation (SPRR1A, KRT family). The upregulation of cell cycle genes indicates HSIL cells are rapidly proliferating without normal regulatory control. In infections with HPV16, E7 was shown to be primarily responsible for this, as it binds and inactivates Rb, freezing E2F to transactivate S-phase genes and override cell cycle checkpoints [65].
This molecular transformation reflects HPV’s transition from productive infection to transformative integration, with viral oncoproteins driving uncontrolled proliferation and genomic instability. Notably, HSIL and LSIL do not share any significant pathways (OC = 0). In fact, they exhibit almost opposite dominant signatures: LSIL is defined by viral gene expression, whereas HSIL is characterized by host cell gene dysregulation. This underscores a fundamental molecular discontinuity between these stages, marking a critical threshold in disease progression.
Our transcriptome analysis reveals dynamic interplay between molecular drivers and passengers across lesion progression: from foundational driver activation (EGFR/ERBB2) in early stages, to mitotic deregulation (CDK1/PLK1) in HSIL, culminating in replication stress and licensing deregulation (MCMs) that characterize high-grade lesions. This progression suggests distinct therapeutic windows for intervention at each stage.
While genomic mutation analysis can identify structural alterations, only transcriptomic profiling can capture the dynamic functional consequences of these mutations and viral integration events. Our approach reveals not just which genes are altered, but how their expression patterns change across disease progression, providing a functional readout of disease activity that surpasses the capabilities of DNA-based or cytology testing alone.
4.3 Candidate biomarkers
We identified a novel class of transcriptome-based biomarkers repressed in a stage-specific manner in response to HPV E6/E7 oncoproteins. Unlike genotyping or morphology, these biomarkers report the actual dynamic biological state of the cervical epithelial cells. Here, we propose specific gene sets aligned to key stages of cervical lesion progression as candidate biomarkers: For NILM, indicative of viral persistence with immune evasion, candidates include RPL5, RPS6, and EGFR, with downregulated IL1β and CXCL1. The ASCUS stage, marking oncogenic initiation with suppressed inflammation, candidates include CCND1, SHH, PARP1, and IGF1, with downregulated IL8. LSIL, characterized by active viral replication and antimicrobial suppression, candidates include DROSHA and RPL8, plus downregulated MPO and DEFA4. Finally, HSIL, associated with genomic instability and loss of epithelial identity, candidates include CDK1, PLK1, MCM4−6, and EZH2, alongside downregulated KRT1 and SPRR1A. Critically, these signatures not only delineate the current lesion stage but also hold potential for predicting future trajectory (Supplementary S11 Table). For instance, IL1β is one central key player in the immune surveillance interactome by driving inflammation and connecting innate and adaptive immunity. In HPV-infected NILM samples, HPV E6/E7 oncoproteins interfere with the normal processing of pro-IL-1, resulting in reducing the secretion of IL-1β in virus-infected keratinocyte, a key mechanism that helps HPV evade immune evasion and initiate malignancy [69]. The SHH ligand enhances the proliferative and migratory capacity of cervical cancer cells and contributes to the development of radioresistance during cancer treatment [70,71]. Increased expression of EZH2 has been linked with the carcinogenesis, proliferation, biology behavior, and the clinical outcome of cervical cancer. High EZH2 level was correlated with more advanced clinical stage, histologic differentiation, infiltration depth, and lymph node metastasis. Patients exhibiting EZH2-positive expression tend to have a decreased overall survival compared with those with EZH2-negative expression [72].
4.4 E6/E7 mRNA expression: the critical link between viral activity and disease progression
Our analysis of E6/E7 mRNA expression levels (measured by PCR Ct values) provides compelling evidence for the superiority of transcriptomic approaches over conventional HPV DNA testing. We demonstrated a significant association between viral oncogene expression and cytological abnormalities (X^2^ = 24.407, p = 0.003). Spearman Rank Correlation showed a significant negative correlation between cytology categories and Ct value (ρ = −0.321, p < 0.001), suggesting that Ct values (indicating higher viral mRNA expression) tend to decrease as cytological abnormalities more severe. These findings directly validate the central role of active viral oncogene expression rather than mere viral presence in driving disease progression.
The transformative power of transcriptomics is particularly evident when comparing our E6/E7 expression data with conventional HPV genotyping results. While 44.7% of our specimens showed multiple HPV strain infections, we found no significant relationship between infection pattern (single vs. multiple) and cytological categories (X^2^ = 0.249, p = 0.983). This striking contrast reveals the fundamental limitation of DNA-based testing: the presence or number of HPV genotypes tells us little about disease activity or progression risk, whereas oncogene expression levels directly correlate with pathological changes.
The current paradigm of HPV DNA testing fails on a crucial front—it cannot distinguish between transient infections that will resolve spontaneously and transformative infections with high oncogene expression that will progress. Our transcriptomic approach demonstrates that E6/E7 mRNA quantification provides critical information about disease activity that DNA testing fundamentally cannot capture. The presence of HPV DNA without significant E6/E7 expression likely represents a biologically distinct and clinically benign disease state compared to infections with active viral oncogene transcription, explaining why approximately 90% of HPV infections resolve without progression while a minority advance to cancer.
Our comprehensive transcriptomic analysis further revealed that the molecular cascade initiated by E6/E7 expression—including growth factor receptor activation, immune evasion, and metabolic reprogramming—provides far more clinically relevant information than simply detecting viral DNA. This functional readout of viral oncogene activity and host response offers the precision needed for truly personalized risk assessment and intervention strategies that current DNA-based testing methods cannot achieve.
4.5 Viral DNA integration
The integration of high-risk human papillomavirus (HPV) DNA into the host genome is a pivotal event in the oncogenic transformation of cervical epithelial cells, marking the transition from persistent infection to precancerous and ultimately invasive carcinoma. Integration typically disrupts the viral E2 regulatory gene, leading to deregulated, persistent overexpression of E6 and E7, inactivation of p53 and Rb, genomic instability, and unchecked proliferation. In our study, although integration events were not directly sequenced, we experimentally verified active viral transcription by RT-qPCR quantification of E6/E7 mRNA in HPV-positive samples, supporting the functional consequence of E2 dysregulation and aligning with the integration-associated transcriptomic changes observed in HSIL. Consistent with E6/E7-driven effects, HSIL showed hyperactivation of cell-cycle/replication programs (e.g., CDK1, PLK1, MCM3–6, RFC3/4; log2FC*>*2, p < 0.001) together with induction of DNA damage response pathways, including ATR/ATM signaling (hsa176187, FDR = 2.22e-9), G2/M checkpoint control (hsa69473), homologous recombination, and mismatch repair, reflecting a cellular environment attempting to manage oncogene-mediated replication stress. These alterations indicate chronic DNA damage and stalled forks that HPV exploits to sustain replication, ultimately fostering host genomic rearrangements (oncogene amplifications, tumor-suppressor disruptions). This integration-linked HSIL profile contrasts with the largely productive, episomal viral state in LSIL, where minimal pathway overlap (OC = 0) suggests limited integration and instability. HSIL enrichment in DNA repair/replication stress responses (WP707, FDR = 0.0021) and mitotic dysregulation underscores integration as a key trigger of malignant progression, potentially facilitated by inflammation-derived reactive oxygen species and chromatin state changes. Our findings, in line with prior reports [23,24,62–64] highlight integration as a hallmark of transformative infection in HSIL and point to biomarker opportunities targeting integration sites or downstream effects (e.g., oncogenic fusions) for risk stratification and early intervention in HPV-driven cervical pathogenesis.
4.6 Transcriptomics: revolutionary advancement beyond current ”stone age” diagnostic tools
Current cervical cancer screening relies primarily on cytology and HPV DNA testing—approaches that detect abnormal cell morphology and viral presence but fail to capture the complex molecular dynamics driving disease progression. These methods represent ”stone age” tools in the era of precision medicine, unable to distinguish between transient infections and those on trajectory toward malignancy.
Our transcriptomic analysis represents a quantum leap in diagnostic capability by revealing the functional alterations that precede and drive morphological changes. Unlike conventional methods, transcriptomics captures the dynamic interplay between viral oncogene expression, host immune response, and cellular transformation pathways.
This approach provides not just diagnostic information but offers prognostic insights by identifying active driver events that promote progression.
The remarkable molecular distinctiveness of HSIL, with minimal pathway overlap with earlier stages (OC = 0.067 with NILM, OC = 0.043 with ASCUS, OC = 0 with LSIL), demonstrates how transcriptomics can delineate disease states with unprecedented precision. Such molecular stratification could revolutionize clinical management by identifying patients requiring aggressive intervention versus those who can be safely monitored.
4.7 Clinical implications: new era of molecular risk stratification
The identification of stage-specific molecular signatures has profound implications for cervical cancer prevention and management. Current approaches treat all HPV-positive patients within a cytological category as having equivalent risk, leading to both over-treatment and under-treatment.
Our transcriptomic profiling enables molecular risk stratification within each cytological category. For example, HPV-positive NILM patients with activated EGFR/ERBB2 signaling and downregulated immune mediators might warrant more aggressive surveillance than those without these signatures. Similarly, ASCUS cases with oncogenic driver activation (CCND1, SHH) might benefit from early intervention rather than surveillance.
The early driver events we identified represent potential biomarkers for progression risk assessment and targets for preventive intervention. Small molecule inhibitors targeting EGFR/ERBB2 signaling, immunomodulatory agents that restore immune detection, or metabolic modulators that target enhanced protein synthesis could potentially disrupt the carcinogenic cascade before irreversible genomic destabilization occurs.
4.8 Future directions: pap smear transcriptomics as the new frontier
Our successful transcriptomic profiling of routine Pap smear specimens demonstrates the feasibility of incorporating molecular analysis into standard clinical practice. Future developments should focus on translating these comprehensive signatures into clinically applicable biomarker panels that can be readily implemented in screening programs.
Prospective studies are needed to validate the predictive value of these molecular signatures for progression across the spectrum of HPV-induced lesions. Integration of transcriptomic data with other molecular features (methylation, microRNA profiles) could further refine risk prediction and therapeutic targeting. The dramatic molecular discontinuity between LSIL and HSIL warrants further investigation to identify the critical molecular switches that mediate this transition. Understanding these pivotal events could reveal novel targets for preventing progression to high-grade lesions.
Our findings call for a paradigm shift in how we conceptualize HPV-induced carcinogenesis, not as a continuous progression but as a series of distinct molecular states with critical thresholds and driver events that can be targeted for prevention and early intervention.
4.9 Study limitation
Although variations in group sizes and the absence of experimental validation may introduce certain constraint, the integrated transcriptomic approach applied in this study provides valuable insights into the molecular mechanisms underlying cervical lesion progression. Future work will extend these findings through immunohistochemical, proteomic, and in vivo validations of representative hub proteins to substantiate their biological significance and clarify their roles in disease pathogenesis. While our network analyses integrated transcript-to-protein associations inferred from high-confidence STRING interaction scores (*≥*0.7), future work will include experimental protein-level confirmation. Key HSIL hub genes (CDK1, PLK1, EZH2, MCMs) and early-lesion immune markers will be validated through immunohistochemistry (IHC) or multiplex immunoassays in a prospective cohort. These studies will establish RNA-to-protein concordance and strengthen the translational relevance of the identified molecular networks. Future integration of targeted proteomic validation will further corroborate the transcriptome-derived network signatures of cervical lesion progression.
5 Conclusion
In conclusion, transcriptome analysis of cervical specimens represents a revolutionary advancement beyond current screening methods, providing unprecedented insight into the molecular evolution of HPV-induced neoplasia and opening new avenues for precision prevention and treatment of cervical cancer.
Supporting information
S1 TablePrimers for E6/E7 mRNA detection of high-risk human papillomavirus.(XLSX)
S2 TableHub genes in NILM.(XLSX)
S3 TableHub genes in ASCUS.(XLSX)
S4 TableHub genes in LSIL.(XLSX)
S5 TableHub genes in HSIL.(XLSX)
S6.1 TablePathway enrichment analysis of upregulated hub genes in NILM.(XLSX)
S6.2 TablePathway enrichment analysis of downregulated hub genes in NILM.(XLSX)
S7.1 TablePathway enrichment analysis of upregulated hub genes in ASCUS.(XLSX)
S7.2 TablePathway enrichment analysis of downregulated hub genes in ASCUS.(XLSX)
S8.1 TablePathway enrichment analysis of upregulated hub genes in LSIL.(XLSX)
S8.2 TablePathway enrichment analysis of downregulated hub genes in LSIL.(XLSX)
S9.1 TablePathway enrichment analysis of upregulated hub genes in HSIL.(XLSX)
S9.2 TablePathway enrichment analysis of downregulated hub genes in HSIL.(XLSX)
S10 TableHPV-Related pathways enriched in four cytology categories.(XLSX)
S11 TableBiomarker candidates in cytological categories.(XLSX)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Society AC. cancer.org. [cited 29-05-2025]. https://www.cancer.org/content/dam/cancer-org/research/cancer-facts-and-statistics/annual-cancer-facts-and-figures/2020/special-section-cancer-in-adolescents-and-young-adults-2020.pdf
- 2Society AC. cancer.org;. [cited 29-05-2025]. https://www.cancer.org/cancer/types/cervical-cancer/about/key-statistics.html
- 3zur Hausen H. Papillomaviruses and cancer: from basic studies to clinical application. Nat Rev Cancer. 2002;2(5):342–50. doi: 10.1038/nrc 798 12044010 · doi ↗ · pubmed ↗
- 4Rodríguez AC, Schiffman M, Herrero R, Wacholder S, Hildesheim A, Castle PE, et al. Rapid clearance of human papillomavirus and implications for clinical focus on persistent infections. J Natl Cancer Inst. 2008;100(7):513–7. doi: 10.1093/jnci/djn 044 18364507 PMC 3705579 · doi ↗ · pubmed ↗
- 5Draft Recommendation: Cervical Cancer: Screening — United States Preventive Services Taskforce — uspreventiveservicestaskforce.org;. [cited 30-05-2025]. https://www.uspreventiveservicestaskforce.org/uspstf/draft-recommendation/cervical-cancer-screening-adults-adolescents
- 6Ramírez AT, Valls J, Baena A, Rojas FD, Ramírez K, Álvarez R, et al. Performance of cervical cytology and HPV testing for primary cervical cancer screening in Latin America: an analysis within the ESTAMPA study. Lancet Reg Health Am. 2023;26:100593. doi: 10.1016/j.lana.2023.100593 37766799 PMC 10520426 · doi ↗ · pubmed ↗
- 7Szarewski A, Mesher D, Cadman L, Austin J, Ashdown-Barr L, Ho L, et al. Comparison of seven tests for high-grade cervical intraepithelial neoplasia in women with abnormal smears: the Predictors 2 study. J Clin Microbiol. 2012;50(6):1867–73. doi: 10.1128/JCM.00181-12 22422852 PMC 3372127 · doi ↗ · pubmed ↗
- 8Gunes AC, Ozgul N, Turkyılmaz M, Kara F, Unlu F, Ayhan A, et al. Evaluation of colposcopy after the addition of human papillomavirus testing to the Turkish cervical cancer screening program. Cancer Med. 2023;12(24):21751–60. doi: 10.1002/cam 4.6740 37994572 PMC 10757080 · doi ↗ · pubmed ↗