The dataset associated with CC [The Cancer Genome Atlas (TCGA)-CC] was obtained by searching for the keywords ‘cervical cancer’ in TCGA (https://portal.gdc.cancer.gov/) database. The dataset (accession no. GSE56363) was obtained from the Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/geo/) database. The inclusion criteria were as follows: i) Patients with CC who received DDP; and ii) survival information or response status to DDP were recorded.

The expression profile of CC and clinical data were obtained by searching for ‘cervical cancer’ in the TCGA database from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/), and consisted of 178 CC tissues and three adjacent non-tumor tissues. Based on the clinical information of the patients in the TCGA-CC dataset, there were 43 patients with both the response status to DDP and overall survival (OS) status recorded. These patients were selected as the discovery set (TCGA-CC₁ set; Table I) to identify the genes associated with DDP resistance. The 40 samples, which only recorded the OS of patients receiving DDP were used as the validation set (TCGA-CC₂ set; Table I) to support the association of genes with DDP resistance. To exclude the prognostic association of the genes, the 95 patients that did not receive treatment were selected as the control set (TCGA-CC₃ set; Table I) for survival analysis (24). GSE56363 consisted of 12 CC samples with complete response to DDP and 9 CC samples with non-complete response to DDP.

RNA-sequencing data were downloaded from TCGA via the GDC Data Portal (https://portal.gdc.cancer.gov/), which had been detected using the Illumina HiSeq 2000 platform. The fragments per kilobase of transcript per million mapped read values were log2-scaled plus 1 for gene expression level measurements.

To identify key genes associated with CC cell survival, the CRISPR-Cas9 screening data of CC cell lines were downloaded from the DepMap portal (https://depmap.org/portal/) by selecting ‘Version: DepMap Public 21Q2’ and ‘CRISPR_gene_effect’ sections. The database recorded the gene essentiality scores [CRISPR-Cas9 gene knockout scores (CERES)] of genes in CC cell lines, which indicated the influence of knockout genes on the proliferation in CC cell lines (25,26). The lower the CERES score, the greater the effect after the gene knockout.

To validate the association of genes with DDP resistance, the gene expression profiles of CC cell lines and their half-maximal inhibitory concentration (IC₅₀) values for DDP drugs were acquired from the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org; release-8.2) database (27) by selecting the ‘Cell Line Gene Expression Data’ and ‘Drug Sensitivity Data’ sections.

Relevant literature was used to identify 30 m⁶A-associated genes (28–31), including 11 methyltransferases, two demethylases and 17 reader proteins (Table II).

HeLa, a human CC cell line, was purchased from Macgene Biotechnology (https://www.macgene.com/). HeLa cells were routinely cultured in Dulbecco's modified Eagle's medium (DMEM; Wuhan Servicebio Technology Co., Ltd.), which was supplemented with 10% fetal bovine serum (FBS; Zhejiang Tianhang Biotechnology Co., Ltd.). Cells were grown at 37°C and 5% CO₂ under humidified conditions and passaged upon reaching 80–90% confluency.

Cell viability was investigated using the Cell Counting Kit-8 (CCK-8; cat. no. C0038; Beyotime Institute of Biotechnology) assay. Cells were seeded at a density of 1×10⁴ cells/ml in a 96-well plate at a volume of 100 ml/well. Various concentrations (0–100,000 nM) of DDP (cat. no. P4394; Sigma-Aldrich; Merck KGaA) were introduced into the culture medium, with a three-fold gradient to systematically probe the cytotoxic effects. After a 96-h incubation at 37°C, cell viability was quantified using the CCK-8 assay and measuring the absorbance, which was used to calculate the cell survival rate. The subsequent data were fitted to a dose-response curve to determine the IC₅₀ of cell proliferation. The equation used to calculate inhibition (%) was: Inhibition (%)=[(A_c-A_s)/(A_c-A_b)] ×100. ‘A_s’ and ‘A_b’ represent the absorbance of the experimental wells and the wells with the highest concentration, respectively. ‘A_c’ represents the absorbance of the control wells.

HeLa cells were initially treated with 1 µM DDP which was increased to 2 µM after ~2 months and treatment was continued at this concentration for another 4 months until stabilization, resulting in DDP-resistant cells (HeLa/DDP). Subsequently, HeLa/DDP cells were seeded at a density of 5×10⁵ cells/well into 6-well plates and maintained in culture medium containing 2 µM cisplatin at 37°C. Next, HeLa/DDP cells were cultured in the presence of increasing concentrations of DPP (cat. no. P4394; Sigma-Aldrich; Merck KGaA) to establish the IC₅₀. The drug sensitivity of the cells were quantified by determining the IC₅₀ using a cell viability assay. The resistance index (RI) was calculated as the ratio of the IC₅₀ of the resistant cells to the IC₅₀ of the parental cells, which served as a measurement of the relative resistance. An RI >3 indicated that the resistant cell line was less sensitive to the drug compared with the parental cell line.

WB was used to detect HYOU1 protein levels in three independent experiments. HeLa and HeLa/DDP cells were harvested and lysed in Whole Protein Extraction kit (cat. no. WLA019, Wanleibio Co., Ltd.) for 5 min. The supernatant was centrifuged at 4°C and 10,005 × g for 10 min and the protein concentration was determined using a bicinchoninic acid kit. Following this, 40 µg of protein from the supernatant was loaded per lane on a 10% gel and SDS-PAGE was carried out before the proteins were transferred to a PVDF membrane. Subsequently, the membrane was blocked with blocking buffer (cat. no. WLA066; Fast Blocking Western; Wanleibio Co., Ltd.) for 1 h at room temperature and then incubated with either the HYOU1 (cat. no. R383157; 1:500; Chengdu Zen-Bioscience Co., Ltd.) or the β-actin (cat. no. WL01372; 1:1,000; Wanleibio Co., Ltd.) primary antibody overnight at 4°C. The membranes were then rinsed with TBST (0.15% Tween20; Wanleibio Co., Ltd.) and incubated with a secondary antibody (cat. no. WLA023; 1:5,000; Goat Anti-Rabbit IgG/HRP; Wanleibio Co., Ltd.) for 45 min at 37°C. Subsequently, the membrane was washed with TBST six times and visualized using Ultrasensitive ECL Chemiluminescence Kit (cat. no. WLA006; Wanleibio Co., Ltd.) (32). The total protein concentration obtained was 2 µg/µl. The intensity of each band was quantified using Gel-Pro-Analyzer software (version 4.0; Media Cybernetics, Inc.).

All small interfering RNA (siRNA), with a final concentration of 50 nM, were transiently transfected into HeLa/DDP cells using Lipofectamine^®™ 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) for 20 min to form transfection complexes at 37°C. Following a 6-h incubation, the transfection medium was replaced with fresh growth medium. DDP was added the next day and the culture was continued for 48 h at 37°C. Transfection efficiency was semi-quantified using WB. The siRNA sequences (Wanleibio Co., Ltd.) used were as follows: HYOU1 sense: 5′-AAGCUGCUGCGUGAGGCUAAUC-3′; anti-sense: 5′-GAUUAAGCCUCACGAGCAGCUU-3′; HYOU1 siRNA-2 sense: 5′-AGCUGGGGAAGAACAUCAAU-3′; anti-sense: 5′-AUUGUUCUUCCCAUCAUCG-3′; and siRNA negative control (NC) sense: 5′-AUAAACAUCGACUCAAU-3′; anti-sense: 5′-AUUGAGCUCGAUUGUUAU-3′.

An unpaired student's t-test was used to identify differentially expressed genes (DEGs) between tumor and normal samples. OS was defined as the time from the date of initial surgical resection to the date of mortality or last contact (censored), which was truncated to 60 months. As the number of responders and non-responders may not be equal, the ‘surv_cutpoint’ algorithm was used to determine the optimal cut-off to distinguish between the high and low expression levels of genes. Survival curves were drawn using the Kaplan-Meier method and statistically compared using the log-rank test. A univariate Cox regression model was used to analyze the association between clinical factors and OS. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated using Cox regression models.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using gene set enrichment analysis (GSEA) from the web-based gene set analysis toolkit (WebGestalt; http://www.webgestalt.org) (33), with a cut-off value of <0.05 for the false discovery rate (FDR). The m⁶A sites of genes were predicted using a sequence-based RNA adenosine methylation site predictor (SRAMP) program (http://www.cuilab.cn/sramp/) (34) by inputting the sequences of the genes. The RMBase version 2.0 platform (http://rna.sysu.edu.cn/rmbase/) (35), a comprehensive resource for RNA modification data verified using methylated RNA immunoprecipitation sequencing, m⁶A-sequence and/or m⁶A-crosslinking immunoprecipitation arrays, was used to validate whether the predicted m⁶A sites underwent m⁶A modification. Subsequently, the interaction probabilities between predicted m⁶A site sequence motifs and the protein sequence of a m⁶A-associated gene were retrieved using the RNA-Protein Interaction Prediction (RPISeq) database (http://pridb.gdcb.iastate.edu/RPISeq/) (36). This database calculated the interaction probabilities using random forest (RF) and support vector machine (SVM) methods.

The correlation between gene expression levels and IC₅₀ values for DDP in the GDSC database was estimated using Pearson correlation analysis and the ggplot2 package in R (https://cran.r-project.org/web/packages/ggplot2/index.html) was used to visualize the results. Comparisons between two groups were analyzed using the unpaired student's t-test. Comparisons among multiple groups were analyzed using one-way analysis of variance (ANOVA) and Tukey's test. P-values were adjusted using the Benjamini-Hochberg procedure for multiple testing (37) to control for the FDR. FDR <0.05 for multiple testing or P<0.05 was considered to indicate a statistically significant difference.

Based on CRISPR-Cas9 screening data from CC cell lines, 699 genes were identified with potential impact on cell proliferation in CC cell lines, in which the CERES scores were <-1 in >75% of CC cell lines. Compared with normal samples, 3,309 DEGs were identified in 43 samples with CC derived from the TCGA-CC₁ set (unpaired student's t-test; FDR <0.05 and log₂(FC)>0; Fig. 1A). Furthermore, 401 DDP resistance genes were identified in the non-response group compared with those in the response group (unpaired student's t-test; P<0.05 and log₂(FC)>0; Fig. 1B). Three genes, including HYOU1, nonsense-mediated mRNA decay factor (SMG5) and ankyrin repeat and LEM domain containing 2 (ANKLE2), were selected as they were significantly upregulated in samples with CC and in the non-response group when compared with normal samples and the response group, respectively (Fig. 1C).

The area under the curve of HYOU1, SMG5 and ANKLE2 for predicting the response and non-response status was 0.802, 0.815 and 0.775, respectively (Fig. 1D-F). Finally, for each gene, the mean expression level was used to stratify patients into high- and low-expression groups and a survival analysis was performed. The results showed that there was no significant difference in the OS between the two groups for the three genes [HYOU1 (high vs. low expression, 19 vs. 24; log-rank P=0.5412; HR=1.50; 95% CIs, 0.40–5.62), SMG5 (high vs. low expression, 23 vs. 20; log-rank P=0.5557; HR=1.51; 95% CIs, 0.38–6.07) and ANKLE2 (high vs. low expression, 20 vs. 23; log-rank P=0.6183; HR=1.40; 95% CIs, 0.37–5.22); Fig. S1].

It was hypothesized that the mean value may not be suitable for distinguishing patients with different responses to DDP. Therefore, the ‘surv_cutpoint’ algorithm was used to re-determine the optimal threshold for HYOU1 expression levels, which was 4.9094. Survival analysis using the TCGA-CC₁ set showed that patients with high HYOU1 expression levels (>4.9094) had a significantly reduced OS compared with patients with low HYOU1 expression levels (<4.9094) following DDP treatment (high vs. low expression, 10 vs. 33; log-rank P=0.0456; HR=3.59; 95% CIs, 0.94–13.67; Fig. 1G). Similarly, the ‘surv_cutpoint’ algorithm was used to re-determine the optimal thresholds for SMG5 and ANKLE2, which were 4.6751 and 3.0310, respectively. However, high or low SMG5 expression levels (threshold, 4.6751) and ANKLE2 expression levels (threshold, 3.0310) did not indicate a significantly different OS in the TCGA-CC₁ set [SMG5 (high vs. low expression, 30 vs. 13; log-rank P=0.0761; HR=291,560,949.96; 95% CIs, 0-infinity (inf); Fig. 1H) and ANKLE2 (high vs. low expression, 35 vs. 8; log-rank P=0.2170; HR=229,801,985.83; 95% CIs, 0-inf; Fig. 1I)]. Therefore, HYOU1 was selected for follow-up analyses as a gene associated with the survival of CC cells and DDP resistance.

In the TCGA-CC₂ set, the ‘surv_cutpoint’ algorithm was used to determine the optimal threshold for HYOU1, which was 4.9094. The 8 patients with high HYOU1 expression levels (>4.9094) demonstrated a significantly reduced OS compared with the 32 patients with low HYOU1 expression levels following DDP treatment (log-rank P=0.0012; HR=7.09; 95% CIs, 1.81–27.70; Fig. 2A). Using the TCGA-CC₃ set, high and low HYOU1 expression levels did not indicate a significantly different OS in patients that did not receive DDP treatment (high vs. low expression, 19 vs. 76; log-rank P=0.6254; HR=1.49; 95% CIs, 0.30–7.38; Fig. 2B). Additionally, according to the GDSC database, the expression levels of HYOU1 were significantly positively correlated with the IC₅₀ values of DDP in CC cell lines (Pearson's correlation analysis; P=0.0384; r=0.58; Fig. 2C).

To validate the effect of HYOU1 on the DDP resistance of CC, HeLa/DDP cells were constructed. The parental HeLa cells exhibited an IC₅₀ of 1.65 µM. By contrast, the resistant cells had an IC₅₀ of 15.51 µM, corresponding to an RI of 9. The IC₅₀ values were determined using dose-response curves generated from cell viability assays (Fig. 2D). Using western blotting, the protein bands revealed an increased HYOU1 expression level in HeLa/DDP cells across three experiments compared with that in HeLa cells (Fig. 2E) and the semi-quantification values in Table SI further elucidates this. The results showed that the protein expression of HYOU1 was significantly increased in HeLa/DDP cells compared with that in parental HeLa cells (unpaired student's t-test; P=0.0002; Fig. 2F). To confirm the efficacy of HYOU1 knockdown, knockdown efficiency was assessed. Using WB analysis, a significant reduction in protein expression levels of HYOU1 was observed in the knockdown groups (one-way ANOVA; P<0.001; Fig. S2), indicating the success of HYOU1 knockdown. Based on this effective knockdown, it was further revealed that HYOU1 knockdown significantly reduced the viability of DDP treated cells compared with the control (one-way ANOVA; P<0.001; Fig. 2G). These results suggest that high HYOU1 expression levels are associated with resistance to DDP.

To further investigate the function of HYOU1, 2,952 genes that significantly correlated with the expression of HYOU1 were identified (Pearson correlation analysis; FDR <0.05; |r|>0.3). These genes were notably enriched in 12 functional pathways (GSEA; FDR <0.05; Fig. 2H). Among these functional pathways, ‘protein processing in endoplasmic reticulum’ and ‘N-glycan biosynthesis’ were significantly enriched in genes that positively correlated with HYOU1 and were involved in DDP resistance (9,13) (Fig. 2I). These results suggest that upregulated expression of HYOU1 is associated with the accumulation of unfolded proteins, and may enhance the stress response in the ER and induce DDP resistance.

Previous preliminary studies report that m⁶A modifications are present in almost all types of RNA molecules in the cell, and regulate the transcriptome to influence RNA splicing, translation, export, localization and stability (18–20). To investigate whether the expression of HYOU1 was regulated by m⁶A modification, the online tool SRAMP was used to predict m⁶A modification sites on HYOU1. This revealed six HYOU1 sequence motifs with high confidence (Fig. 3A; Table III).

The correlation between m⁶A-associated genes and the expression of HYOU1 using the TCGA-CC₁ set was analyzed and eight m⁶A-associated genes were found that significantly correlated with the expression of HYOU1 (Pearson's correlation analysis; FDR <0.05; |r|>0.4; Fig. 3B; Table II). Among these genes, the expression of EIF3A was significantly upregulated in the non-response group compared with that of the response group (unpaired student's t-test; P=0.0399; FC=1.07; Fig. 3C). Furthermore, the ‘surv_cutpoint’ algorithm was used to determine the optimal thresholds for EIF3A, which was 5.2442. Survival analysis indicated that patients with high EIF3A expression levels (>5.2442) had a significantly reduced OS compared with patients with low EIF3A expression levels (<5.2442) following DDP treatment using TCGA-CC data integrated with TCGA-CC₁ and TCGA-CC₂ sets (high vs. low expression, 35 vs. 48; log-rank P=0.0310; HR=2.81; 95% CIs, 1.05–7.48; Fig. 3D). In an independent dataset of patients with CC (GSE56363), the expression of EIF3A was significantly increased in the non-response group compared with the response group (unpaired student's t-test; P=0.0228; FC=1.04; Fig. 3E).

Sequence docking prediction analyses with the RPISeq database confirmed, with high probabilities and confidence, that the EIF3A reader may bind with the six m⁶A site motifs of HYOU1 (interaction probabilities >0.5; Table IV; Fig. 3F), including the ‘3294’, ‘8651’, ‘10147’, ‘10786’, ‘11220’ and ‘11607’ sites. Furthermore, searching for the HYOU1 gene on the RMBase version 2.0 platform revealed that the m⁶A site (‘3294’) of HYOU1, which exhibited a high probability of binding with EIF3A, was modified by m⁶A modification (Table III).