To evaluate the differentially expressed HPV-related gene signature between tumor and normal tissues, the online web server GEPIA2 was applied to analyze the gene expression data from TCGA cohorts and GTEx projects. Differential signature score analysis revealed that the HPV-related gene signature was significantly reduced in the tumor tissues of eight cancer cohorts as compared to the expression levels in corresponding non-tumor tissues, including Cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), Breast invasive carcinoma (BRCA), Acute Myeloid Leukemia (LAML), Ovarian serous cystadenocarcinoma (OV), Skin Cutaneous Melanoma (SKCM), Testicular Germ Cell Tumors (TGCT), Uterine Corpus Endometrial Carcinoma (UCEC), Uterine Carcinosarcoma (UCS). In the cohort of CESC, 306 tumor tissues and 13 normal tissues were compared (Figure 1A). A pan-cancer analysis with 10,967 tumor samples in 32 TCGA cohorts revealed that a low frequency of ALDH1A2 (1.4%), OSR2 (5%), GRIA4 (3%), GATA4 (4%), and IRX4 (5%) was genetically altered. The high frequency of genetic alterations in the five genes was observed in several cancers, including esophageal adenocarcinoma (30.77%), bladder urothelial carcinoma (29.68%), uterine carcinosarcoma (28.07%), lung adenocarcinoma (27.56%), lung squamous cell carcinoma (27.52%) (Figure 1B). 13.8% of 297 cervical squamous cell carcinoma cases revealed that genetic alterations in the five-gene signature with a low frequency of ALDH1A2 (1.3%), OSR2 (1%), GRIA4 (4%), GATA4 (1%), and IRX4 (7%) (Figure 1C). Together, the gene expression level of HPV-related gene signature was significantly lower in cervical cancer as compared to the level of normal tissues. In contrast, a rather low frequency of genetic alterations was observed.

To investigate the promoter methylation data of HPV-related gene signature in cervical cancer, we applied an interactive web-portal UALCAN to perform analyses for TCGA-CESC. Box whisker plots revealed that promoter DNA methylation levels (average beta values) of ALDH1A2, GATA4, GRIA4, and IRX4 were significantly elevated in primary tumors as compared to normal tissues (Figure 2A). It is generally accepted that gene promoter hypermethylation may downregulate transcription of the affected gene. We observed a significant inverse correlation between gene promoter methylation and relative transcript levels for ALDH1A2, OSR2, and IRX4 (Figure 2B). Differential expression analysis revealed that the levels of ALDH1A2 and OSR2 were significantly reduced in the tumor tissues as compared to the normal tissues (Figure 2C). Thus, our findings suggest that ALDH1A2, OSR2 gene transcription is regulated by gene promoter methylation, which is suppressed in the primary tumors. In addition, promoter DNA methylation levels (average beta values) of ALDH1A2, OSR2, GATA4, GRIA4, and IRX4 were compared among the patients with different clinical stages. Box whisker plots showed that promoter DNA methylation levels of ALDH1A2 in tumor clinical stage IV were significantly elevated than any other subgroup with stage I/II/III (Figure 3A). Concerning average beta values of GATA4, GRIA4, and IRX4, the patients with tumor stage IV had the highest level as compared to the subgroup patients with stage I/II/III tumor (Figures 3C–E). By contrast, OSR2 promoter DNA methylation level in patients with stage I tumors was significantly higher as compared to the patients with stage III/IV tumors (Figure 3B). In summary, promoter methylation of HPV-related gene signature is significantly associated with the CC clinical stage.

Next, we asked whether HPV-related methylation-based gene signature is associated with the clinical prognosis of patients with cervical cancer. A heat map based on ALDH1A2, OSR2, GATA4, GRIA4, and IRX4 transcript levels were generated by the SurvExpress web tool, and two clusters with either low or high prognostic risk were defined by prognostic index (PI) and Cox fitting (20, 21) (Figure 4A). The expression level of the five genes was compared between the high-risk group (n=86) and the low-risk group (n=175), which revealed the levels of the five genes were significantly reduced in the high-risk group (Figure 4B). As expected, the age of patients and the number of death events in the high-risk group are significantly higher as compared to the low-risk group (Figures 4C, D). Chi-square analysis showed the low-risk group is significantly associated with younger age (P=0.002) and lower T status (p=0.006) (Supplementary Table S1). Furthermore, the Kaplan-Meier plot and Log-Rank analysis exhibited a significant difference in overall survival between the high and low-risk groups (Figure 4E; Log-rank = 16.446, P < 0.001). In summary, these findings demonstrated that HPV-related methylation-based gene signature predicts the clinical outcome for CC patients.

To confirm the prognostic value of the candidate gene ALDH1A2 in the protein level, we determined ALDH1A2 protein levels by immunohistochemical (IHC) staining in primary tumor tissues from patients with cervical cancer, and we found the heterogeneous cytoplasmic expression of ALDH1A2 in tumor tissues (Figure 5A). The patients were divided into two subgroups depending on different immunoreactivity scores (IRS) with ALDH1A2^high (n = 13) and ALDH1A2^low (n = 37) for further analysis. The patients with clinical stage I/II had a statistically higher IRS of ALDH1A2 as compared to the patients with stage III/IV (Figure 5B), which was in line with our in silico analysis findings. Survival analysis revealed that patients with ALDH1A2^low expression tumors had a significantly worse disease-specific survival (DSS) as compared to the patients with ALDH1A2^high expression pattern (Figure 5C). Taken together, our results showed low expression of ALDH1A2 serves as an unfavorable risk factor and a prognosticator to identify patients with cervical cancer at high risk for treatment failure.

Mutation and CNA analyses were conducted on the TGCA PanCanAtlas datasets using the cBioPortal database version 3.6.12 (Center for Molecular Oncology at MSK, New York, NY, USA). DNA methylation data and gene expression data from TCGA-CESC were accessed through cBioPortal tool. Correlation analysis was performed between methylation beta values and log-base-2-transformed gene expression data with Pearson’s correlation with a Bonferroni correction to the P-values. Expression in various cancers and their normal tissue counterparts were analyzed using the Gene Expression Profiling Analysis (GEPIA) (Beijing, China). GEPIA2 provides RNA sequencing data from TCGA of tumor samples with paired adjacent TCGA and Genotype-Tissue Expression (GTEx) normal tissue samples. TCGA and GTEx RNA-Seq expression datasets in GEPIA are based on the UCSC (University of California, Santa Cruz) Xena project (34), which are recomputed based on a uniform bioinformatic pipeline to eliminate batch effects. To compare expression data, data are normalized by quantile-normalization or other two additional normalization strategies. UALCAN database (http://ualcan.path.uab.edu) was used to compare the promoter methylation of HPV-related gene signature between normal and primary CC tissues, as well as methylation beta values of each gene among different tumor stages.

We downloaded the five genes expression profile of TCGA-CESC from cBio Cancer Genomics Portal (http://cbioportal.org) in December 2020. Low-risk and high-risk groups of patients with CC were generated depending on the prognostic index (PI) and Cox fitting. A heat map representation of the five genes expression values was made by supervised clustering with ClustVis (21), a freely available web server at http://biit.cs.ut.ee/clustvis/.

Patients with primary CC diagnosed and treated between 2015 and 2020 were comprised in the retrospective study cohort. Tumor tissue sections from formalin-fixed paraffin-embedded (FFPE) were obtained at the Department of Pathology, The First Affiliated Hospital of University of Science and Technology of China after approval by the ethics committee. Biopsies of non-surgically treated patients, as well as samples of patients who underwent tumor surgery, were included. All subjects gave written informed consent for data collection as it is a standard procedure in our department. Patients with suspicious clinical findings who underwent diagnostic endoscopy and/or patients before tumor surgery with a histologically confirmed diagnosis of CC were asked to consent. The protocol was approved by the Ethics Committee of the First Affiliated Hospital of USTC (Ethic No 2021/BLK02) in accordance with the declaration of Helsinki. Experimental treatment procedures were not part of this study. The patients were treated according to the guidelines for cervical cancer. The final analysis was based on 50 patients with CC (Supplementary Table S2). Clinical and therapeutic follow-up of the cohort was assessed retrospectively. The cohort did not include recurrences or specimens from surgery.

Immunohistochemical staining was carried out in paraffin-embedded tissue from patients with primary CC. Tissue sections were deparaffinized and rehydrated by the following our previous study (35) and incubated with primary antibody ALDH1A2 (1:500, HPA010022, Sigma-Aldrich) at 4°C overnight. Tissue sections were incubated with a biotinylated anti-rabbit secondary antibody (Boster) for 2 h. Streptavidin-HRP (Boster) and Peroxidase Substrate (DAB) solution (MXB Biotechnologies, Fuzhou, China) were added to the tissue sections for signal development. Tissue slides were scanned using the DMS-10 Scan System (D-metrix, Suzhou). Three experienced pathologists analyzed scanned slides by using the D-metrix Viewer software (version 1.0.0). Semiquantitative analysis was performed according to the number of stained tumor cells (score A; 1, no positive cells; 2, less than 33% positive cells; 3, between 34% and 66% positive cells; 4, more than 66% positive cells), and according to the staining intensity (score B; 1, no staining; 2, weak staining; 3, moderate staining; 4, high staining). Both values were multiplied to calculate the final immunoreactivity score (IRS) for defined subgroups (final expression scores ranged 1-16).

All statistical analysis was done using GraphPad Prism version 9.1 (GraphPad Software, USA) and IBM SPSS Statistics version 22. Differences of clinicopathological features between the low-risk and high-risk groups were compared using Chi-square Test or Fisher’s exact test. Significance of the difference was calculated by student’s t-test. Disease-specific survival (DSS) was estimated as the time from the date of primary tumor diagnosis to the date of CC-related death within the follow-up period (events) or to the date of CC unrelated death or without progression (censored). Survival difference was assessed by Kaplan-Meier plot and Log-Rank analysis. A P-value <0.05 was considered statistically significant.