RNA-seq reads data, methylation 450 k data, and clinical data of patients with cervical carcinoma from TCGA were acquired from Broad Institute Firehose (https://xenabrowser.net/datapages/). ARGs were obtained from the Human Autophagy database (HADb http://www.autophagy.lu/project.html). Samples without clinical survival status follow-up data were excluded from the study. Overall, a total of 296 CESC samples were enrolled in the present study. The 31 core ARGs involved in the further analysis are described in Table 1. These ARGs, identified as important regulators of autophagy, were essential for autophagosome formation and were involved in the further result analysis. Rstudio and R programming(4.03) were the principal tools for analyzing data throughout the study.

The optimal cutoff values of ARGs based on the expression Fragments Per kilobase per Million mapped reads (FPKM) value or methylation Beta-value, the survival time, and the survival status were identified by the “surv_cutpoint” function in “survminer” package as described before (Wang et al., 2021). The cutoff value of the OS significant genes is listed in Supplementary Figure S1. The patients with relative values above or below the optimal cutoff were considered as high or low groups for each gene expression level, methylation level, or EMT score. Kaplan–Meier (KM) survival analysis with log-rank test was then used for investigating the OS difference between the abovementioned high and low groups. The correlation between the OS characteristics of CESC patients and each ARG expression was explored by a univariate Cox regression analysis. A multivariate Cox regression analysis was performed to investigate whether each ARG is an independent OS influence factor.

The patients with a relative value above or below the optimal cutoff were considered as high or low group for each gene expression level, methylation level, or EMT score. KM survival analysis with log-rank test was then used for investigating the overall survival (OS) difference between the abovementioned high and low groups.

By “GSVA” R package, the EMT score for evaluating the EMT status of each CESC patient was quantified by single-sample gene-set enrichment analysis (ssGSEA) based on previously reported EMT signature genes (Choi et al., 2010; Hanzelmann et al., 2013). The correlation between ATG5 and EMT score or between ATG5 and EMT signature related genes is calculated by Pearson’s correlation coefficient and plotted in heatmap by the “circlize” R package. The EMT score density between normal and tumor tissues and the survival difference between high and low groups were also analyzed.

The genes co-expressed with ATG5 in CESC were screened according to the Pearson correlation coefficient (|cor| > 0.3, p < 0.05) and plotted in the volcano plot by the “ggplot2” package. Heatmaps of the top 50 negative or the top 50 positive ATG5 expression correlated genes were plotted by the “pheatmap” R package. Metascape (Zhou et al., 2019) was employed to gain insights into the biological functions of those co-expressed genes of ATG5.

The methylation site information of ATG5 in TCGA CESC methylation data based on Illumina Infinium Human Methylation450 BeadChip platform was mapped by an annotation file retrieved from the hg19 GPL16304 legacy annotation file (Price et al., 2013). The methylation level of each site in ATG5 was grouped into low and high by the survminer package as mentioned above. The OS significant associated sites were screened through KM survival analysis. Besides, MEXPRESS was used for visualizing expression, DNA methylation, and clinical parameters based on TCGA data (http://mexpress.be). The correlation between the methylation level of each site of ATG5 and the expression of ATG5 is calculated by the “corrplot” package (Wei and Simko, 2021).

Waterfall plots of ARGs mutation were produced by the oncoplot function in the “maftools” R package (Mayakonda et al., 2018) to visualize the mutation status such as classification and frequency of mutation types, frequency of variant types, and frequency of SNV classes for CESC samples. The mutational exclusion and co-occurrence mutations among those ARGs were analyzed by Fisher’s exact test (Mayakonda et al., 2018). The variant allele frequency (VAF) distribution of ARGs was also plotted.

ALL CESC patients were grouped based on different clinical staging information (pathological—T/N/M, histologic—Grade/Stage) or the age at diagnosis (<45 years, 45–69 years, and >69 years). Subtypes of CESC were grouped into squamous cell neoplasms or adenocarcinoma based on the clinical data. The OS difference between the high- or low-ATG5 expression groups was evaluated in KM survival analysis.

Hela, Ca Ski (human cervical cancer cell lines), and PANC-1 cells (pancreatic cancer cell line) were purchased from the American Type Culture Collection (ATCC, Manassas, VA, United States). Hela cells were incubated in DMEM supplemented with 10% FBS (Gibco, Grand Island, NY). Ca Ski and PANC-1 cells were incubated in RPMI 1640 containing 10% FBS (Gibco, Grand Island, NY). PANC-1 cells were used as a positive control for ATG5 expression (Yang et al., 2015). Cells are cultured in an incubator at 37°C with 5% CO₂.

According to the manufacturer’s instructions (GeneChem Corporation, Shanghai, China), Hela and Ca Ski cells were transfected with the recombinant lentiviral vector LV-ATG5-RNAi labeled with green fluorescent to establish ATG5-knockdown cells and named as KD. In brief, the cells were digested at the logarithmic growth phase, then resuspended in the complete medium at a concentration of 3–5 ×10⁴ cells/ml, and seeded in six-well plates. At 20% confluence, cells were transfected with the recombinant lentivirus for 16 h. Subsequently, a conditioned medium containing transfection reagent was refreshed. At 72 h post-infection, the positive expression rate of GFP was used to evaluate the transfection efficiency under a fluorescence microscope (Olympus Corporation, Tokyo, Japan). siRNAs specific for ATG5 were obtained from GenePharma (Shanghai, China), and the detailed sequences were as follows: ATG5-RNAi (98917-1) sense: CCTTTCATTCAGAAGCT GTTT, ATG5-RNAi (98918-1) sense: CCTGAACAGAA TCATCCTTAA, ATG5-RNAi (98919-1) sense: GAT​TCA​TGG​AAT​TGA​GCC​AAT, NC-RNAi sense: TTC​TCC​GAA​CGT​GTC​ACG​T.

Cells were collected and lysed by Trizol reagent (Invitrogen, Carlsbad, CA). Next, cDNA was obtained by using M-MLV reverse transcriptase (Promega, Madison, WA). Using a LightCycler 480 II RT-PCR System (Roche, Basel, Switzerland), qRT-PCR was performed as follows: 95°C for 30 s, and then 40 cycles of 95°C for 5 s and 60°C for 30 s. The sequences of primers amplified in this study were as follows: ATG5: forward 5′-AAGAGTGTTTATTCG TCGGT-3′, reverse 5′-ATC​ACA​GCT​TAG​TGT​TCC​CT-3′. ACTB: forward 5′-GCG​TGA​CAT​TAA​GGA​GAA​GC-3′, reverse 5′-CCACGTCACA CTTCATGATGG-3′. Beta-actin (ACTB) gene was used as an internal reference.

In short, cells were harvested and washed with PBS twice, and then lysed in RIPA buffer (Beyotime, Shanghai, China) on ice for 15 min. Then, BCA Protein Assay Kit (Beyotime, Shanghai, China) was used to detect the protein concentration in the supernatant from crushed cells by ultrasonication. The final concentration was adjusted to 2 μg/μl by fresh RIPA solution, followed by protein denaturation. 10% SDS-PAGE was applied to separate the protein samples and then transferred onto PVDF membranes (Millipore, Burlington, MA). After that, membranes were incubated with TBST containing 5% skim milk to block non-specific antigens. Next, the membranes were incubated with primary antibodies and HRP-conjugated second antibodies in turn. Finally, the protein bands were visualized by the electrochemiluminescence (ECL) detection system. The primary antibodies used include mouse anti-ATG5 antibody (1:500, sc-133158, Santacruz, Santa Cruz, CA), mouse anti-ERK antibody (1:2,000, #9107,CST, Danvers, MA, United States), rabbit anti-NFκB p65 antibody (1:3,000, #8242, CST, Danvers, MA, United States), rabbit anti-mTOR antibody (1:1,000, #2983, CST, Danvers, MA, United States), rabbit anti-p-ERK antibody (1:1,000, #4376, CST, Danvers, MA, United States), rabbit anti-p-NFκB p65 antibody (1:500, ab76302, Abcam, Cambridge, MA, United States), rabbit anti-p-mTOR antibody (1:500, ab109268, Abcam, Cambridge, MA, United States), rabbit anti-P38 antibody (1:3,000, #8690, CST, Danvers, MA, United States), rabbit anti-AKT antibody (1:500, #4691, CST, Danvers, MA, United States), mouse anti-E-Cadherin antibody (1:500, #14472, CST, Danvers, MA, United States), rabbit anti-p-P38 antibody (1:5,000, #4631, CST, Danvers, MA, United States), rabbit anti-p-AKT antibody (1:1,000, #4060, CST, Danvers, MA, United States), rabbit anti-N-Cadherin antibody (1:500, ab18203, Abcam, Cambridge, MA, United States), mouse anti-Twist antibody (1:100, ab50887, Abcam, Cambridge, MA, United States), rabbit anti-c-Myc antibody (1:1,000, ab32072, Abcam, Cambridge, MA, United States), rabbit anti-Snail antibody (1:500, #3879, CST, Danvers, MA, United States), rabbit anti-β-Catenin antibody (1:2,000, #8480, CST, Danvers, MA, United States), rabbit anti-Slug antibody (1:500, #9585, CST, Danvers, MA, United States), rabbit anti-p-β-Catenin antibody (1:500, #2009, CST, Danvers, MA, United States), rabbit anti-MMP-2 antibody (1:500, #40994, CST, Danvers, MA, United States), mouse anti-Fibronectin antibody (1:500, MAB 1918, R&D Systems, Minneapolis, MN, United States), rabbit anti-Vimentin antibody (1:500, #5741, CST, Danvers, MA, United States), and rabbit anti-MMP-9 antibody (1:500, #13667, CST, Danvers, MA, United States). Mouse anti-β-actin (1:5,000, sc-69879, Santacruz, Santa Cruz, CA) and mouse anti-GAPDH (1:10,000, ab37168, Abcam, Cambridge, MA, United States) were used as an internal control.

Cells were seeded in the upper chamber coated with (invasion) or without (migration) Matrigel (Corning Costar, Cambridge, MA, United States) at a density of 1 × 10⁵ cells/well. In invasion assay, serum-free medium was added into the lower chamber. In the migration assay, a medium containing 30% FBS was added to the lower chamber. The noninvasive/non-migrated cells were removed by a cotton swab after 48-h incubation. Then, the cells were fixed in 4% paraformaldehyde for 0.5 h and stained with 0.5% crystal violet solution, followed by taking photographs through the microscope.

All molecular biology assays were repeated in triplicate. Experimental data are presented as the mean ± SD. R (version 4.03) and Rstudio software were applied to perform statistical analysis. The differences between the groups were analyzed via Student’s t-test and Mann–Whitney U test. Univariate and multivariate Cox regression analysis were performed to identify the factors influencing prognosis. Survival analysis was done using the KM method and the log-rank test. p < 0.05 was designated as statistically significant.

A total of 309 sequencing data and 312 methylation data in TCGA CESC cohort were obtained from the TCGA data. The median age of those 306 participants was 46.57 years (min 20.95 years and max 89.04 years). The median follow-up time for CESC patients with survival information was 1.87 years (range 0–17.55 years), and among them, 225 were alive and 74 were dead during the follow-up time.

High expression levels of ATG4C, ATG5, BECN1, and WIPI1 were associated with poor prognosis in CESC, while CESC patients with overexpression of ATG3, ATG4A, ATG4B, ATG4D, ATG7, ATG9B, ATG13, GABARAP, and GABARAPL2 have longer OS and higher 5-years survival rate (Figure 1). Furthermore, the univariate Cox regression results indicated that ATG5, ATG4D, GABARAP, ATG4A, and ATG3 were significantly associated with OS (Figure 1N). ATG5 was identified as a risk prognostic predictor in CESC patients (HR 1.7; 95% CI, 1.0–2.8, p = 0.047), while another four ARGs were identified as favorable prognostic factors (HR < 1, p < 0.05). Subsequently, the multivariate Cox results showed that ATG5, ATG4D, and ATG4A were independent prognostic factors affecting OS of patients with CESC (Figure 1O). Among them, high expression of ATG5 or low expressions of ATG4D and ATG4A were associated with unfavorable prognosis. These results indicated that ATG5 among these ARGs was the most harmful factor affecting the prognosis of CESC patients. According to the KM analysis, the 5-years survival rate for high- and low-expression groups of ATG5 is 0.486 (0.375–0.631) and 0.782 (0.708–0.863) respectively. Furthermore, KM survival analysis was also performed to show that ATG5 expression could effectively predict the prognosis of patients with CESC with different clinical stages, pathological grades, and ages at diagnosis (Figures 2A–F). For patients with cervical squamous cell carcinoma, patients with ATG5 overexpression have significantly shorter OS than the ATG5 down-expression group (Figure 2G). For cervical adenocarcinoma, due to the limited sample size of patients included in the analysis, the p-value of the km survival curve is not significant, but patients with ATG5 down-expression showed a trend towards a longer survival with no death during the follow-up times (Figure 2H).

Based on the TCGA database, we analyzed the co-expressed genes of ATG5 in 304 CESC patients through the Pearson correlation coefficient. As shown in Figure 3A, red dots represented 5,666 genes positively associated with ATG5, while blue dots represented 2,114 genes negatively associated with ATG5. The heat map was used to show the top 50 genes positively or negatively correlating with ATG5 (Figures 3B,C). Next, GO enrichment analysis was performed to investigate the biological pathway of the top ATG5-related co-expressed genes. ATG5 co-expressed genes are correlated with ribosome, cytoplasmic, C1 complex, TRBP-containing complex, glucosamine-containing compound catabolic process, and so on. There are concerns that ATG5 also participated in the regulation of the immune effector process and regulation of T-cell-mediated immunity (Figure 4A). We further found that the ATG5 co-expressed genes involved in the regulation of the immune effector process were A2M, C1QA, C1QB, C1QC, RPS19, MAP3K7, TYROBP, STX7, and IL20RB. MAP3K7, STX7, and IL20RB were involved in the regulation of T-cell-mediated immunity. Moreover, ATG5 expression was inversely associated with MAP3K7 and STX7 and positively related to A2M, C1QA, C1QB, C1QC, RPS19, TYROBP, and IL20RB (Figures 4B–J). Furthermore, it was important to investigate the relation between ATG5 and EMT-related gene signature. The result is represented in Figure 5. We found that ATG5 was significantly correlated with EMT scores quantified by ssGSEA and EMT signature-related genes (Choi et al., 2010) such as DPH3, CD24, KRT8, CTNNAL1, EPCAM, TSPAN13, MYB, TPD52, ECHDC1, CLSPN, USP33, TTK, TOM1L1, GCA, COMMD8, CGAS, OSTM1, and EMT score. Tumor tissues had higher EMT scores than normal tissues, and a high EMT score predicted shorter survival time for CESC patients.

Genetic alteration results showed that 16.15% of patients own mutations of these ARGs. Among these genes, a total of 16 genes have a mutation rate ≥1%, of which ATG2B and RB1CC1 are the most frequently mutated genes (4%). Missense and splice are the two most common types of mutations (Figure 6A). ATG5 was altered in 3 (1.03%) of the 291 CESC patients. Figure 6B shows the mutation type and domain of ATG5. Moreover, concomitant occurrence of ATG5 mutation and ULK2 mutation was found in CESC (Figure 6C). Compared with most mutated ARGs, ATG5 owns a relatively low VAF, except ATG10. MEXPRESS was used for visualizing expression, DNA methylation, and clinical parameters based on TCGA data, and showed the dramatic relationship between ATG5 methylation status and OS in CESC patients (Figure 7A). According to the annotation file of Methylation450 BeadChip platform, there are 19 methylation sites of ATG5, namely, cg24211478, cg14951955, cg00333154, cg07541085, cg15336269, cg07747241, cg01786360, cg16950012, cg22876713, cg10566763, cg26178529, cg01644481, cg02711268, cg14660784, cg18206,736, cg20059537, cg21148531, cg04389950, and cg10850197. As shown in Figure 7B, there was a negative correlation between ATG5 expression and the methylation sites cg02711268, cg18206,736, cg07541085, cg24211478, cg26178529, and cg10850197. The significant sites of those methylation sites in ATG5 on OS in CESC patients were shown in KM including beneficial OS relevance methylation sites such as cg00333154, cg15336269, cg02711268, and cg21148531 and the harmful OS-related sites such as cg24211478, cg10566763, cg26178529, and cg10850197 (Figures 7C–J).

ATG5 mRNA showed higher expression in cervical cancer cell lines Hela and Ca Ski compared with the cervical epithelial cell line CRL2614 (Figure 8A). The recombinant lentiviral vector containing ATG5-RNAi labeled with green fluorescence was used to downregulate ATG5 expression of Hela and Ca Ski cells. Cells presented with green fluorescence indicated the success of transfection (Figure 8B). RT-qPCR showed that the levels of ATG5 mRNA were significantly downregulated in Hela-KD (knockdown) and Ca Ski-KD cells compared with NC (negative control) cells (Figure 8C, p < 0.01). In addition, the sequence with the highest knockdown efficiency was selected for use in the next study. The knockdown efficiency was further confirmed by Western blot (Figure 8D).

By Transwell migration assay, ATG5 knockdown demonstrated a markedly negative effect on the migratory capacity of Hela and Ca Ski cells compared with NC (Figure 9). Subsequently, the role of ATG5 in invasion was assessed by a Transwell invasion assay. Likewise, the data of the Transwell invasion assay showed that downregulation of ATG5 dramatically suppressed the invasion in cervical cancer cells (Figure 9). In brief, these results exhibited the repression of ATG5 in the migration and invasion of cervical cancer cells in vitro.

Further investigation found that knockdown of ATG5 in Hela cells downregulated expressions of interstitial phenotypic markers (N-cadherin and Vimentin), EMT signaling pathway-related proteins (P-ERK, P-NFκBp65, and P-mTOR), extracellular matrix metalloproteinases (MMP-9), CC-specific interstitial phenotypic marker (Fibronectin), and EMT-related transcription factors (Twist), but upregulated epithelial marker E-cadherin (Figure 10). However, other markers including NFκBp65, mTOR, P38, AKT, P-P38, P-AKT, MMP-2, β-Catenin, Slug, P-β-Catenin, c-Myc, and Snail were not affected by ATG5 knockdown. Generally, downregulation of ATG5 could reverse the EMT process by specific pathways in cervical cancer cells and resulted in attenuation of migration and invasion.