Three human CC cell lines, HeLa, SiHa, and C-33 A, were obtained from the American Type Culture Collection (ATCC, United States) and cultured in Dulbecco’s Modified Eagle Medium (DMEM) (Hyclone, United States) supplemented with 10% fetal bovine serum (FBS) (BI, United States) at 37°C with 5% CO₂.

Cell viability was assessed using the Cell Counting Kit-8 (CCK-8, Sparkjade) assay. GA was dissolved in methanol at a concentration of 500 mM and then diluted to different concentrations (0 or control, 100 and 200 µM) using appropriate amounts of methanol. At the logarithmic growth stage, the cells were trypsinized, neutralized, and then seeded at a density of 5 × 10³ cells in 96-well plates containing complete DMEM. After 12–24 h, the culture medium was replaced with diluted GA. The absorbance of the plates was measured at 24, 48, and 72 h posttreatment at 450 nm. The data were analyzed using GraphPad Prism 9.0.0 software.

HeLa cells were cultured with the methanol control and a medium containing 200 μM GA (control and GA groups, respectively) for 3 days. When the cells reached a fusion of 60%–80%, they were collected and washed twice with PBS. Total RNA was extracted and purified from six samples (n = 3/group) using TRIzol Reagent (Invitrogen). RNA quantification and sequencing were performed using the BGISEQ-500 platform (BGI, China). The raw sequencing data were deposited in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the accession number PRJNA1169526.

A total of four mRNA expression datasets obtained using human CC cell lines [GSE7803 (ZHAI et al., 2007), GSE9750 (Scotto et al., 2008), GSE63514 (den Boon et al., 2015), and GSE527 (WONG et al., 2003)] were selected and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). All datasets included analyses of the differences between normal cervical and CC tissue samples. GSE7803 is based on GPL96, [HG-U133A] Affymetrix Human Genome U133A array platform, which contains 10 normal squamous cervical epithelium and 21 invasive squamous CC samples. GPL96 was also used as the platform for GSE9750, which contains 24 normal cervical and 33 CC samples. The platform for GSE63514 is GPL570, [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, which contains 24 normal cervical epithelium and 28 cervical squamous epithelial cancer specimens. GSE527 is the platform for GPL355, Human 10 K cDNA Array, which contains eight normal cervical and 25 CC samples. The RNA-seq results of HeLa cells treated with GA were provided by BGI.

GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) is an online analytical tool within the GEO dataset designed to analyze the DEGs between normal cervical and CC tissue samples. Genes were selected according to the following criteria: |logFC| ≥ 1, P-value < 0.05, and adj. P-value < 0.05.

The gene expression profiling interactive analysis (GEPIA2) (http://gepia2.cancer-pku.cn/#degenes) online database was used to analyze the DEGs in TCGA as follows: CESC, |Log2FC|Cutoff: 1, q-value Cutoff: 0.05, Differential Methods: ANOVA. The criteria for screening DEGs based on the RNA-seq results were as follows: DESeq2: |log2FC| ≥ 0.5, Q-value ≤ 0.05.

First, the significantly upregulated or downregulated genes among the four datasets of the GEO database were combined and duplicated, and the overlapping genes were removed. Then, a combined analysis of the results from the GEO and TCGA databases was used to screen common upregulated and downregulated genes between the CC and normal cervical samples. The duplicated genes were removed. Finally, by combining the results of the GEO and TCGA databases and RNA-seq, we identified genes that were upregulated or downregulated in CC compared with normal cervical samples, while the expression patterns of these genes were the reversed in GA-treated CC cells. The final differentially expressed genes are referred to as the “DEGs” in the following text for convenience.

The identified DEGs were subjected to Gene Ontology (GO, biological process, cellular component, and molecular function) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the Metascape database (https://metascape.org/gp/index.html). The enrichment criteria were as follows: Min overlap ≥3, P-value cutoff ≤0.05, min enrichment ≥1.5, and default values for the rest.

The Search Tool for the Retrieval of Interacting Genes (STRING, https://cn.string-db.org/) was employed to analyze and visualize the interaction relationship between proteins. The final set of DEGs was uploaded to STRING using the following settings: Meaning of network edges: confidence, Minimum required interaction score: high confidence (0.700), Network display options: disable structure previews inside network bubbles, hide disconnected nodes in the network.

Gene sets with interaction relationships were imported into the Cytoscape software to create a visible network. First, the Gentiscape2.2 plug-in was employed to calculate the network and topological characteristics of each node. The genes corresponding to the nodes with degree value ≥Mean + SD and Betweenness value ≥Mean + SD nodes were referred to as hub and bottleneck genes, respectively. With the MCODE plug-in, the parameters were set as degree cutoff ≥3, K-core ≥4, and default values for the rest. The core genes that constitute the stable network structure were screened out. Finally, the potential therapeutic candidate genes that were in common among the hub, bottleneck, and core genes were identified using the “Calculate and draw custom Venn diagrams” online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

The total RNA was reverse transcribed to cDNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific), which was used as a template for qRT-PCR amplification. The PCR mixtures were prepared to a final volume of 20 µL using the PowerUpTM SYBRTM Green Master Mix (Applied Biosystems) with four wells for each sample. The PCR assay was performed using the ABI Q5 real-time PCR platform (ABI Life Technologies), and the cycling parameters were as follows: 1 cycle of 50°C for 2 min, 1 cycle of 95°C for 2 min, 40 cycles of 95°C for 15 s and 60°C for 1 min. Fold changes were calculated and normalized via the ΔΔCt method using the glyceraldehyde 3-phosphate dehydrogenase gene as the internal normalization control. The primers used were as follows: (5ʹ-3ʹ): CDC20 (F: CTG​GAT​CAA​AGA​GGG​CAA​CTA, R: GGC​AGA​GTG​ACT​GGT​CAT​ATT); DLGAP5 (F: GTT​GTG​CAG​CCT​GTA​ATG​CC, R: TAG​CAG​CTC​TTG​TGA​CTG​GC); KIF20A (F: TGC​TGT​CCG​ATG​ACG​ATG​TC, R: AGG​TTC​TTG​CGT​ACC​ACA​GAC); and GAPDH (F: CAA​TGA​CCC​CTT​CAT​TGA​CC, R: GAC​AAG​CTT​CCC​GTT​CTC​AG).

The candidate genes were submitted to the MirDIP (http://ophid.utoronto.ca/mirDIP/index_confirm.jsp), miRDB (https://mirdb.org/), and ENCORI (https://rnasysu.com/encori/) online databases to predict their upstream interaction with miRNAs. Then, the reliability of the miRNAs was verified using TarBase V 9.0 software (https://dianalab.e-ce.uth.gr/tarbasev9). The miRNAs were submitted to the ENCORI online tool for analyzing the relationships between miRNA and lncRNA (parameters setting: Clip-data ≥2 datasets, Degradome Data ≥1), and miRNA and circRNA (parameter setting: Clip-data ≥2 datasets, Degradome Data ≥3, and default value for the rest).

The relationships between the upstream lncRNAs and circRNAs with miRNAs, miRNAs with candidate genes, and the candidate genes with pathways that met the screening criteria were evaluated. The file was imported into the Cytoscape software to establish a visual lncRNA/circRNA-miRNA-mRNA-pathway regulatory network.

The data were analyzed using GraphPad Prism 9.0.0 software (La Jolla, CA, United States) (https://www.graphpad.com). The control and GA groups were compared using the Student’s t-test. A significant difference (P < 0.05) was indicated between the groups based on the means +SEM.

The CCK-8 assay was used to assess the proliferation of HeLa, SiHa, and C-33 A cells treated with different GA concentrations (0, 100, or 200 µM) was assessed. The results indicated that GA significantly inhibited the proliferation of CC cells (Figure 1).

Four raw expression microarray datasets of CC (GSE7803, GSE9750, GSE63514, and GSE527) were analyzed using GEO online. After filtering out 64 genes with opposite expression trends, we screened 5,234 DEGs between normal cervical and CC tissue samples, including 2,806 upregulated and 2,428 downregulated genes (Figure 2A).

In the TCGA database, 6,072 genes were screened out (three duplicates and one gene with opposite expression removed), including 2045 upregulated and 4,027 downregulated genes. Combining the results of both databases, we screened the genes with the common expression trend, which included 629 upregulated and 630 downregulated genes (Figure 2B).

As GA significantly inhibited the proliferation of CC cells, we treated HeLa cells with GA and then performed an RNA-seq analysis to identify the target genes. A total of 2,343 DEGs were screened out, including 1,501 upregulated and 842 downregulated genes (Figure 2C).

Combining the results from the GEO and TCGA databases with RNA-seq, we identified the upregulated and downregulated genes in CC compared with normal cervical samples. After treatment with GA, the expression of these genes was reversed, with downregulated genes becoming upregulated and vice versa. Finally, we screened out 127 DEGs, including 40 upregulated and 87 downregulated genes, in the GEO and TCGA databases, which were downregulated and upregulated after treatment with GA, respectively, showing an opposite trend (Table 1).

The DEGs were uploaded to the Metascape database for GO (biological processes, cellular components, and molecular functions) and KEGG pathway enrichment analysis.

A total of 221 biological processes, 82 cellular components, and 63 molecular functions were enriched. The top 20 clusters of biological processes are shown in Figure 3A, such as negative regulation of organelle organization, regulation of nuclear division, cell morphogenesis, elastic fiber assembly, and positive regulation of organelle organization. The top 18 clusters of cellular components, as shown in Figure 3B, included the extracellular matrix, basal plasma membrane, neuromuscular junction, nuclear envelope, and cytoplasmic side of the membrane. As shown in Figure 3C, the top 19 clusters of molecular function included extracellular matrix structural constituents, anaphase-promoting complex binding, growth factor binding, integrin binding, and organic anion transmembrane transporter activity.

Furthermore, 36 KEGG pathways were enriched, and the top 10 clusters with their representative enriched terms included protein digestion and absorption, the Pl3K-Akt signaling pathway, “Glycine, serine and threonine metabolism,” human T-cell leukemia virus 1 infection, and EGFR tyrosine kinase inhibitor resistance (Figure 4).

The protein–protein interaction relationships corresponding to the DEGs were analyzed and visualized using the STRING online tool. According to the setting parameter criteria, 87 proteins without interaction were hidden, and a protein interaction network with 61 edges of the remaining 40 node-proteins was constructed (Figure 5).

The protein–protein interaction network was imported into the Cytoscape software with the MCODE plug-in. The eight core genes constituting the stable structure of the network were screened out (Figure 6A; Table 2). Using the Gentiscape2.2 plug-in, eight hub and seven bottleneck genes were screened out (Table 2). The core, hub, and bottleneck genes, which were in common and are critical to the network, were considered the candidate genes. These included cell division cycle 20 (CDC20), disks large homolog associated protein 5 (DLGAP5), and kinesin family member 20A (KIF20A) (Figure 6B; Table 2). We used real-time PCR to detect the RNA expression levels of the candidate genes in the HeLa control and GA groups. The results showed that the expression of three candidate genes (CDC20, DLGAP5, KIF20A) matched the RNA-seq analysis (Figure 6C).

Furthermore, the three candidate genes were submitted to the MirDIP, miRDB, and ENCORI online databases to predict their upstream miRNAs. Reliability was verified using TarBase V9.0 software, and 13 miRNAs were screened (Table 3). Subsequently, these miRNAs were submitted to the ENCORI database and four upstream lncRNAs and 42 circRNAs were obtained after removing duplicate genes.

All relationships, including the candidate genes with their associated KEGG pathways, the upstream lncRNAs and circRNAs with the miRNAs, and the miRNAs with the candidate genes, were imported into the Cytoscape software for constructing and visualizing the lncRNA/circRNA-miRNA-mRNA-pathway regulatory network. This network comprised 63 nodes (3 candidate genes, 1 KEGG pathway, 13 miRNAs, 4 lncRNA, 42 circRNAs) and 95 edges (Figure 7).