The original datasets were downloaded from geo database to compare the mRNA expression and miRNA expression between cervical cancer and normal tissues. The expression profile data of GSE7410, GSE63514, GSE86100 and GSE9750 were downloaded based on different platforms (Table 1). The mRNA expression data of cervical cancer were downloaded from TCGA database for preprocessing. Clinical samples related to cervical cancer were selected. The data set includes 307 cervical cancer samples, 307 normal samples and corresponding clinical data. Then, Sangerbox tool 2.0 (http://sangerbox.com/Index) corrects and normalizes the original expression data background of each GEO dataset, and eliminates batch effects and other irrelevant variables.

The differentially expressed mRNAs and miRNAs between cervical cancer and normal tissues were identified by R software limma software package. The GEO dataset filtering criteria are: log2|Fold Change|>1, adjusted p-value< 0.05 and false discovery rate (FDR) < 0.05. The screening criteria of the TCGA dataset were: log2|Fold Change|>2, adjusted p-value< 0.05 and FDR<0.05. The ggplot2 and pheatmap packages were used to draw differential mRNAs and miRNAs volcano maps and heatmaps respectively. Then, a Venn diagram was used to find overlapping differentially expressed mRNAs.

To further analyze differentially expressed mRNAs, GO and KEGG enrichment analyses were performed using the DAVID online tool 6.8 (https://david.ncifcrf.gov/home.jsp) and Metascape (http://metascape.org/gp/index.html). GO enrichment analysis mainly annotates the biological process (BP), cytological component (CC) and molecular function (MF) of genes. KEGG enrichment analysis mainly predicts the signal pathways that may be involved. p < 0.05 was considered statistically significant.

PPI network interactions were identified and constructed using the STRING online database 11.5 (http://string-db.org). An interaction confidence≥0.4 was considered significant. Then, the data were imported into Cytoscape software 3.8.0 for reprocessing. Subsequently, the molecular complex detection (MCODE) plug-in was used to filter hub modules in the PPI network. Then, GO and KEGG analyses of mRNAs in hub modules were performed using Metascape (http://metascape.org/gp/index.html).

The clinical characteristic data of cervical cancer patients from TCGA were downloaded. First, we subtyped the downloaded cervical cancer samples, and then analyzed the survival of each subtype and each pathological parameter. Then, univariate Cox proportional hazards regression analysis was used to identify potential genes highly related to overall survival and to screen mRNAs related to the prognosis of cervical cancer (p < 0.05). Then, the data were further screened by multivariate Cox proportional hazards regression. A risk assessment model was constructed with key prognostic mRNAs as dependent variables to evaluate the clinical value of key prognostic mRNAs in predicting patient survival. The risk assessment model formula is: risk score = expression of mRNA1 × β_mRNA1 + expression of mRNA2 × β_mRNA2 +……+ expression of mRNAn × β_mRNAn. β represents the multivariate Cox regression coefficient, which is obtained from the multivariable Cox proportional hazards regression model of each gene. According to the median risk score, patients with cervical cancer were divided into high and low-risk groups. Then the survival rates of the high and low-risk groups were analyzed. At the same time, a receiver operating characteristic (ROC) curve was constructed.

First, we preprocessed the GSE9750 data, removed the samples with incomplete clinical data, removed the outliers, and performed WGCNA in the Sangerbox tool 2.0 (http://sangerbox.com/Index). The soft threshold was set as β = 1, and a scale-free network was constructed. Then, we transformed the adjacency matrix into a topological overlap matrix (TOM). Next, we performed hierarchical clustering to identify modules, and each module contained at least 20 genes (min module size = 20). Subsequently, characteristic genes were calculated, hierarchical clustering modules were used, and similar modules were merged (abline = 0.25). Finally, we calculated the correlation between the module and clinical data to determine meaningful clinical modules.

We performed GSEA using normalized RNA SEQ data obtained by TCGA. GO and KEGG analyses of candidate genes were performed with Sangerbox tool 2.0 (http://sangerbox.com/Index) to study the possible biological functions of candidate genes. Adjusted p value < 0.05 and FDR <0.05 were considered to be statistically significant.

HeLa cells were purchased from Beijing Beina Chuang lian Biotechnology Research Institute (Beijing, China). The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Invitrogen, United States). At the same time, 10% fetal bovine serum (Invitrogen, United States), 100 U/ml penicillin (Sigma, United States) and 100 g/ml streptomycin (Sigma, United States) were added. The cells were cultured in 5% CO₂ at 37 C in a constant-temperature incubator. The eighth passage HeLa cells were used in this study. According to the instructions for HeLa cells, Lipofectamine 3,000 reagent (Invitrogen, United States) was used to transfect 50 nm siRNA (si-TYMS; RiboBio, China), pCDH-TYMS (RiboBio, China), 50 nm miR-197-3p mimic (mir100000 227-one to five, RiboBio, China), 100 nm miR-197-3p inhibitor (mir200000 227-one to five, RiboBio, China) and their NC control for 24 h.

A total of 30 cervical cancer tissues and 30 adjacent normal tissues were collected (Supplementary Material S1). This study was approved by the hospital ethics committee and was carried out in accordance with the Declaration of Helsinki. In addition, each patient provided written informed consent. The cervical cancer tissue was cut off during surgery and immediately stored at 80°C until use.

Total RNA was extracted from tissues and cells using RNeasy Mini Kit (Qiagen, Germany). Then, the concentration of extracted RNA was determined with Nanodrop 2000 (Invitrogen, United States). The same amount of RNA was reverse transcribed into cDNA with Goldstar™ RT6 cDNA synthesis kit ver2 (TsingKe, China). Master qPCR mix (SYBR Green 1) (Beyotime, China) was used for qPCR detection. The expression of GAPDH was used as the endogenous control. Primer 5.0 software was used to design primers 5.0 (TYMS forward primer 5'-3': ACT​TGT​GCA​GAT​TAT​TCA​GGA​C, TYMS reverse primer 5'-3': ATT​CTT​CTG​TCG​TCA​GGG​T; GAPDH forward primer 5'-3':CAT​TTC​CTG​GTA​TGA​CAA​CGA, GAPDH reverse primer 5'-3': GGG​TCT​TAC​TCC​TTG​GAG​G). The primer sequences were sent to Sangon (China) for synthesis. For miRNA detection, U6 was used as the endogenous control, and the primers were purchased from RiboBiology (miR-197-3p: mqps0000765-1-100, RiboBio, China). QPCR detection was performed using bulge loop miRNA QRT PCR Starter Kit (RiboBio, China). MiRNA and mRNA expression levels were calculated by 2^−ΔΔCt.qRT-PCR data were analyzed by GraphPad Prism, and each reaction was repeated 3 times. A t test was used in both groups (p < 0.05), and one-way ANOVA was used in more than two groups (p < 0.05).

Logarithmic growth HeLa cells were inoculated into 96-well plates at 4,000 cells/well per well. 24 h after transfection, CCK-8 reagent (10 μl per well, Beyotime, China) was added to all cells, and the cells were incubated for 3 h. Then, a multi-function microplate (Thermo, United States) was used to measure the absorbance at 450 nm. Each experiment was conducted three times.

HeLa cells at the logarithmic growth stage were inoculated into 96-well plates at 1×10⁵ cells per well. 24 h after transfection, the fused cell monolayer was scratched with 200 μl sterile pipette tip and then added to serum-free DMEM. The scratch was recorded by microscopy at 0 and 24 h, and the scratch closure rate was evaluated by ImageJ software.

HeLa cells transfected for 24 h were collected for single-cell suspension, and 1×10³ cells were added to each upper chamber (Corning, United States) containing matrix gel. DMEM with 20% fetal bovine serum was added to lower chamber. After 24 h, the cells were fixed with 4% paraformaldehyde (Beyotime, China) and stained with 1% crystal violet (Beyotime, China). Cells were observed and counted under a light microscope.

24 h after transfection, the HeLa cells were digested with trypsin to prepare single cell suspension. Then, the cells were gently washed with precooled PBS 3 times, and the cell concentration was maintained at 3 × 10⁵ cell/ml. The cells were then treated with annexin V-FITC & PI apoptosis detection kit (Sangon, China). After treatment, the number of apoptotic cells analyzed by Beckman flow cytometry (Beckman, United States).

HeLa cells were inoculated into 96-well plates and cotransfected with 100 ng of the dual-luciferase reporter constructs pmiR-RB-Report-TYMS 3'UTR (RiboBio, China) and miR-197-3p mimic (RiboBio, China). After incubation for 48 h, the supernatants of the cells were collected and the fluorescence value was determined using dual-luciferase reporter kit (Promega, United States).

48 h after the transfection of HeLa cells, the culture medium was removed. After PBS rinsing, the cells were lysed with protein lysate (Beyotime, China), and the total proteins were collected. Then the protein concentration was detected by BCA (Beyotime, China) protein quantitative kit. After adding proper amount of SDS loading buffer (Beyotime, China), denaturation was carried out in boiling water at 100°C for 5 min. Then 12% SDS-PAGE electrophoresis (Beyotime, China) was performed. Protein bands were transferred to PVDF membranes (Beyotime, China) by western transmembrane system. Subsequently, the PVDF membrane was sealed in 5% skim milk powder (Beyotime, China) and incubated for 4 h. After washing with TBST (Beyotime, China), the PVDF membrane was incubated with primary antibody at 4°C overnight. After washing, the membrane was incubated with the secondary antibody at room temperature for 60 min. Finally, the PVDF membrane was stained with western TMB substrate (Beyotime, China; P0211) and scanned. The antibodies used were β-actin (Beyotime, China; AF5003), TYMS (Abcam, United States, ab108995) and horseradish peroxidase labeled goat anti-rabbit IgG (Beyotime, China; A0208).

All data were processed by SPSS 25 and GraphPad Prism 8. The measurement data are expressed as mean ± standard deviation (SD). A t test and one-way ANOVA were performed between two groups and multiple groups, respectively. p < 0.05 was considered statistically significant.

Three cervical cancer expression profile datasets GSE7410, GSE63514 and GSE86100 were downloaded from the GEO database. Table 1 shows the details of the three GEO datasets. GSE7410 and GSE63514 screened for differentially expressed mRNAs between cervical cancer and normal tissues, and GSE86100 screened for differentially expressed miRNAs. Based on screening criteria, GSE7410 and GSE63514 obtained 1,138 mRNAs (709 down-regulated and 429 up-regulated) and 1,009 mRNAs (641 down-regulated and 368 up-regulated), respectively (Figure 1A, Figures 1B–D, Figures 1A,B, Figure 1B, Supplementary Material S2 and Supplementary Material S3). GSE86100 obtained 166 miRNAs (74 down-regulated and 92 up-regulated) (Figures 1A–C, Figures 1B,C and Supplementary Material S4). In addition, 822 mRNAs (567 down-regulated and 255 up-regulated) were obtained (Figures 1A–D, Figures 1A,B and Supplementary Material S5) from TCGA-CESC. Subsequently, Venn diagram was used to compare mRNAs screened by GSE7410, GSE63514 and TCGA-CESC, and 34 overlapping differentially expressed mRNAs were obtained, including 22 up-regulated genes and 12 down-regulated genes (Figure 1C, Supplementary Table S1).

To explore the function and pathway of overlapping genes, we performed GO and KEGG pathway analyses. GO analysis showed that, in terms of biological progression, overlapping genes were mainly enriched in the regulation of mitotic nuclear division, T cell chemotaxis, regulation of nuclear division, leukocyte chemotaxis, positive regulation of the release of sequestered calcium into the cytosol, lymphocyte chemotaxis and mitotic nuclear division (Figure 2A). In terms of molecular function, overlapping genes were mainly enriched in CXCR chemokine receptor binding, chemokine receptor binding, chemokine activity, receptor regulator activity, G protein-coupled receptor binding, receptor ligand activity and cytokine activity (Figures 2A,B). Pathway enrichment analysis showed that overlapping genes interact with the Toll-like receptor signaling pathway, the p53 signaling pathway, viral protein interaction with cytokine and cytokine receptors, chemokine signaling pathways, the cell cycle and cellular senescence were related (Figure 2B).

We constructed PPI networks for 34 overlapping genes using the SRTING online database 11.5. The PPI network contains 22 nodes and 48 edges (Figure 2C). Then, two hub subnetworks were filtered through the MCODE plug-in from the PPI network (Figures 2A–C and Figures 2B,C). Subnetwork one contains the genes CXCL1, CXCL9, CXCL10, CXCL13, STAT1 and TLR2. Subnetwork two contains the genes ANLN, CCNB1, CHEK1, RRM2, TYMS and UBE2C. Based on the TCGA-CESC dataset, we mapped the expression profiles of 12 hub genes in cervical and normal tissues (Supplementary Figure S1). GO and KEGG analysis results showed that genes of subnetwork one were closely related to immunity, stress chemotaxis, angiogenesis, Toll-like receptor signaling pathway and chemokine signaling pathway (Supplementary Figure S2A). Genes in subnetwork two were mainly enriched in mitosis, the cell cycle and the p53 signaling pathway (Supplementary Figure S2B).

First, overall survival (OS) data of cervical cancer patients were obtained from TCGA database. The survival curves of various subtypes and pathological parameters of cervical cancer were analyzed (Supplementary Figure S3). Subsequently, three genes were identified by univariate cox proportional risk regression model (p < 0.05) (Supplementary Figure S4, Table 2). Two further prognostic genes: RRM2 and TYMS, were selected by a multivariate Cox proportional risk regression model (Table 2). These two genes were also independent prognostic genes (p < 0.05). Risk prognosis model was established based on multivariate Cox screening results. The risk score was calculated as follows: mRNA risk score = (0.256)×expression (RRM2)+(-0.333)×expression (TYMS). Patients in TCGA-CESC were divided into high and low-risk groups based on the median risk score. The risk score of the TCGA-CESC dataset is shown in Figure 3. Figure 3A shows the details of the risk score. Kaplan-Meier survival analysis showed that the prognosis of high-risk group was worse than that of low-risk group (Figure 3C). The AUCs at 1 year, 3 and 5 years were 0.75, 0.72 and 0.69 respectively (Figure 3B). These results show that prognostic gene markers (RRM2 and TYMS) have good survival prediction ability, indicating that the model can effectively predict the prognosis of cervical cancer patients.

The coexpression network was constructed by (WGCNA) to further identify genes strongly associated with cervical cancer. First, the optimal soft threshold β = 1 was obtained through analysis, and the scale-free network was constructed (Figure 4A). Finally, two gene modules were obtained for further analysis (Figure 4B). From the heatmap of the topological overlap map, it can be seen that there was a high correlation between the genes in the module (Figure 4C). As shown in Figure 4D, the turquoise module was identified as the most specific module with a correlation coefficient of 0.84 (p = 2E-18). Therefore, this module was selected as clinically important for further analysis. By comparing the hub genes screened by multivariate Cox analysis with the genes in the turquoise module, we obtained a common key gene TYMS (Supplementary Table S2). This suggests that TYMS may play a very important role in the progression of cervical cancer.

To investigate the expression of TYMS in cervical cancer, qRT-PCR was used to detect the expression levels of TYMS in cervical cancer tissue samples and adjacent tissues. The results showed that TYMS was significantly down-regulated in cervical cancer (Figure 4E), which was consistent with GEO database analysis results. Subsequently, we downloaded the expression data of TYMS and the clinical characteristics of cervical cancer patients from TCGA-CESC, and evaluated the relationship between the expression level of TYMS and various clinicopathological parameters. The analysis showed that down-regulation of TYMS expression was significantly correlated with histological grading (p < 0.001), clinical stage (p < 0.001) and patient age (p < 0.001) (Figure 4F). These results suggest that cervical cancer patients with low expression of TYMS are more likely to have more advanced tumor status, grade and stage than those with high expression of TYMS.

To further analyze and explore the potential biological functions of TYMS, we performed GSEA on TYMS. As shown in Figure 5A, angiogenesis, myogenesis, inflammatory response and protein levels were significantly positively correlated with TYMS expression at GO concentrations. DNA repair, G2M checkpoint and apoptosis were negatively correlated with TYMS expression (Figures 5A,B). KEGG pathway analysis showed that the pathways with a significant positive correlation with TYMS expression included fatty acid metabolism, the PPAR signaling pathway, fatty acid metabolism, arachidonic acid metabolism and histidine metabolism. The p53 signaling pathway, cell cycle and pyrimidine metabolism were significantly negatively correlated with TYMS expression (Figures 5A,B and Figure 5B). These results suggest that essential cell cycle control, amino acid metabolic pathways and regeneration processes are closely related to TYMS expression in patients with cervical cancer. In addition, we also analyzed the relationships between TYMS expression and immune neoantigens, tumor mutational burden, microsatellite instability and methylase expression. The results showed that the expression of TYMS was negatively correlated with tumor mutational burden (Figure 5C) and microsatellite instability (Figure 5D), and positively correlated with the expression of four methyltransferases (Figure 5E, DNMT1: red, DNMT2: blue, Dnmt3A: green and DNMT3B: purple). These results also show that low-level TYMS may easily leads to cell mutation, resulting in tumor formation or deterioration.

To further verify the analysis results of the biological function of TYMS, gain-of-function analysis was performed on HeLa cells transfected with the TYMS vector. QRT-PCR results showed that TYMS was effectively expressed in HeLa cells after transfection (Figure 6A). Then, the cell proliferation was detected by CCK-8 assay. The results showed that compared with the control group, HeLa cells transfected with si-TYMS showed significant growth promotion. In contrast, TYMS overexpression showed the opposite effect (Figure 6B). The results of wound healing and transwell assay showed that silencing TYMS promoted cell migration and invasion, and overexpression of TYMS inhibited cell migration and invasion (Figure 6C, Figure 6D, Figure 6E and Figure 6F). However, flow cytometry showed that overexpression of TYMS significantly promoted the amount of apoptosis in HeLa cells, while the apoptosis rate was relatively low when TYMS was knocked out (Figure 6G and Figure 6H). These results suggest that TYMS may be involved in the development of cervical cancer.

To determine the potential mechanism of TYMS in cervical cancer, we analyzed the potential miRNAs targeted by TYMS using miRWalk and starbase online tools (Supplementary Material S6 and Supplementary Material S7). The analysis results were compared with the screened differentially expressed miRNAs, and a hub miRNA miR-197-3p was obtained (Figure 7A). The 3'UTR binding site sequences of miR-197-3p and TYMS are shown in Figure 7B. To verify the correlation between miR-197-3p and TYMS, luciferase reporter gene detection was performed. The results showed that miR-197-3p overexpression caused a significant decrease in luciferase activity in the wt-TYMS group. While, overexpression of miR-197-3p failed to inhibit luciferase activity in the mut-TYMS group (Figure 7C). In addition, qRT-PCR and western blot results showed that overexpression of miR-197-3p could reduce the expression of TYMS in HeLa cells. Similarly, the expression of TYMS was increased in HeLa cells after transfection with miR-197-3p inhibitor (Figure 7D, Figure 7E and Figure 7F). These results suggest a direct targeting relationship between TYMS and miR-197-3p.

To further confirm the regulatory relationship between miR-197-3p and TYMS, we detected the mRNA level of miR-197-3p in cervical cancer tissues. Detection results showed that miR-197-3p was up-regulated in cervical cancer tissues (Figure 7G). Then, correlation analysis was conducted on the expression levels of miR-197-3p and TYMS, and the results showed a significant negative correlation between miR-197-3p and TYMS (Figure 7H), which was consistent with the results of bioinformatics analysis (Figure 7I). The results of proliferation assay showed that overexpression of miR-197-3p significantly enhanced the proliferation activity of HeLa cells. In contrast, the miR-197-3p inhibitor inhibited HeLa cell proliferation (Figure 7J). In HeLa cells cotransfected with miR-197-3p inhibitor and si-TYMS, cell proliferation activity was significantly lower than that in the TYMS knockout group (Figure 7J). These results suggested that down-regulation of miR-197-3p could reverse the effect of the TYMS inhibitor on HeLa cell viability and suggest that miR-197-3p directly targets TYMS and negatively regulates its expression.