Analysis of Autophagy-Related Signatures Identified Two Distinct Subtypes for Evaluating the Tumor Immune Microenvironment and Predicting Prognosis in Ovarian Cancer

Ovarian cancer (OC) is one of the most lethal gynecologic malignant tumors. The interaction between autophagy and the tumor immune microenvironment has clinical importance. Hence, it is necessary to explore reliable biomarkers associated with autophagy-related genes (ARGs) for risk stratification in OC. Here, we obtained ARGs from the MSigDB database and downloaded the expression profile of OC from TCGA database. The k-means unsupervised clustering method was used for clustering, and two subclasses of OC (cluster A and cluster B) were identified. SsGSEA method was used to quantify the levels of infiltration of 24 subtypes of immune cells. Metascape and GSEA were performed to reveal the differential gene enrichment in signaling pathways and cellular processes of the subtypes. We found that patients in cluster A were significantly associated with higher immune infiltration and immune-associated signaling pathways. Then, we established a risk model by LASSO Cox regression. ROC analysis and Kaplan-Meier analysis were applied for evaluating the efficiency of the risk signature, patients with low-risk got better outcomes than those with high-risk in overall survival. Finally, ULK2 and GABARAPL1 expression was further validated in clinical samples. In conclusion, Our study constructed an autophagy-related prognostic indicator, and identified two promising targets in OC.


INTRODUCTION
Ovarian cancer (OC) is one of the most lethal gynecologic malignant tumors (1). Due to the nonspecific symptoms in the early stage and the lack of effective screening techniques of the disease, a large number of patients are diagnosed at an advanced stage, of which the 5-year survival rate was less than 30%. Cytoreductive surgery and platinum-and paclitaxel-based chemotherapy are still the basic treatments for OC. Despite advances in combination chemotherapy, targeted therapy and intraperitoneal chemotherapy, 80% of OC patients initially respond to treatment, but most eventually relapse and ultimately develop into a chemotherapy-resistant disease; thus, no significant improvement in the prognosis of OC has been achieved (2)(3)(4).
The clinicopathological features of OC are predicted by the WHO classification and TNM staging system of tumor lymph node metastasis, which is also the key for selecting appropriate treatment. However, because of the heterogeneity of OC, there are obvious stratifications into histological or molecular subtypes, and the results may be significantly different even for patients with similar clinical features and treatment regimens. These observations showed that the clinicopathological features and current classification are not sufficient for prediction and risk stratification. Consequently, it is difficult to meet the needs of clinicians (5,6). Therefore, it is of great significance for improving the prognosis of OC to search for specific prognostic biomarkers and therapeutic targets with higher predictive value.
Tumor immunotherapy has become a promising treatment strategy, which aims to restore the immune response to fight against tumors. Immunotherapies such as immune checkpoint therapy (ICT), tumor vaccines, immune adoptive therapy and immunomodulators have been applied in many cancers. Many immunosuppressive receptors have been identified and studied in tumors; these studies have led to the development of therapies including, but not limited to, FDA-approved monoclonal antibodies that mediate clinically relevant immunostimulatory effects by suppressing immunosuppressive receptors, such as PD-L1, PD-1, CTLA-4, LAG3, TIGIT and BTLA (7)(8)(9). The application of immunotherapy has significantly changed the strategies and modes of treating OC and greatly improved the quality of life in some patients with OC. Pembrolizumab, nivolumab and avirumumab are anti-PD-1 or anti-PD-L1 monoclonal antibodies, and bevacizumab is a monoclonal antibody that binds to vascular endothelial growth factor (VEGF). These drugs have been successfully used to treat recurrent or drug-resistant OC (10)(11)(12). However, only some patients benefit from immune treatment, and some patients still show poor responses or resistance to immunotherapy. How to successfully identify which patients may benefit from immunotherapy and which patients may exhibit poor responses or resistance to immunotherapy is a clinically difficult problem. Therefore, screening subjects suitable for immunotherapy would help increase the success rate of treatment and benefit more patients.
Autophagy is an important immunomodulatory process in the tumor microenvironment that can maintain the homeostasis, activation and biological function of immune cells. Innate immune-mediated autophagy can be upregulated by activating innate immune receptors, including Toll-like receptors (TLRs) and nucleotide oligomeric domain (NOD)-like receptors (NLRs) (13). The adaptive immune response depends on the recognition of extracellular or intracellular peptide epitopes provided by major histocompatibility complex II (MHCII) and MHCI molecules and recognized by CD4 + T and CD8 + T cells, respectively. Autophagy provides ATP molecules for antitumor T cells to activate antigen-presenting cells (APCs) (14). Autophagy also plays a role in protecting cells and tissues from stress in normal physiological processes. However, inappropriate autophagy may lead to low antitumor immunity, affect infiltrating of immune cells, and inhibit the immune response to weaken immunotherapy (15). Recent studies have shown that autophagy is closely related to tumor immunotherapy and clinical prognosis of OC (16). Targeted autophagy may be a promising therapeutic strategy for improving the efficacy of immunotherapy and enhancing the immune response. However, it is still difficult to find effective and appropriate autophagy-related gene biomarkers to identify and evaluate tumor-specific cellular immune responses in tumor patients and to predict patient responses to immunotherapy.
In our study, we downloaded the expression profile of OC patients from the TCGA database and obtained autophagyrelated genes (ARGs) from the MsigDB database. Based on the ARGs, OC patients were successfully classified into two subtypes, and a risk model was established to assess the prognosis of OC. The system comprehensively evaluated the correlation among classification, risk model and TIM and explored the effect of autophagy on the regulation of the OC TIM. ARGs may be potential biomarkers and provide new ideas for immunotherapy. the TCGA data coordination center,307 ovarian cancer samples were included. This dataset shows the gene-level transcription estimates, as in log2(x+1) transformed RSEM normalized count. The genes are mapped onto the human genome coordinates using HUGO probeMap, refering to method description from University of North Carolina TCGA genome characterization center. For all the cohorts, only patients with available expression profiles, clinicopathological (including age, status, lympho vascular invasion, stage) data and survival data were included in the analyses. The eight OC samples and paired adjacent tissue samples were obtained from Third Xiangya Hospital of Central South University.

Identification of OC Subtypes
We obtained the autophagy gene set from the MSigDB (https:// www.gsea-msigdb.org/gsea/msigdb/search.jsp) KEGG_ REGULATION of AUTOPHAGY of molecular signatures; we obtained a total of 35 autophagy-related genes. Compared with the expression data of the genes in the TCGA, candidate genes whose expression level was too low (log2(x+1) <1) were excluded, and 19 ARGs were obtained. We performed k-means clustering based on the mRNA expression data of 19 autophagyrelated genes (17). Before performing k-means, a filtering procedure was conducted. We performed k-means ("kmeans" function in R) clustering and used the "NbClust", "cluster" and "factoextra" packages in R to determine the optimal number of OC subtypes. The values of k where the magnitude of the cophenetic correlation coefficient began to fall were chosen as the optimal number of clusters.

Estimation of Immune Infiltration
The inference of infiltrating cells in the tumor microenvironment (TME) was used to quantify the level of immune cell infiltration in the OC samples. Based on the immune cell marker genes provided by Bindea G,et al (18), we used the R language GSVA package and then used single-sample gene set enrichment analysis (ssGSEA) to quantify the levels of infiltration of 24 immune cells, including T lymphocytes, dendritic cells, natural killer cells, etc., into the sample based on the immune cell marker gene expression profile data in the TCGA-OV. According to the level of immune cell infiltration, we divided patients into a highinfiltration group and a low-infiltration group, observed the relationship between age, stage, grade, lymphatic metastasis and immune infiltration, and used a heat map display to observe various immune cells in the high-and low-risk groups. Finally, the Pearson correlation coefficient was used to calculate the correlation coefficient among immune cells, risk genes and immune cells.

Screening of Differentially Expressed Genes (DEGs) and Bioinformatics Analysis
To obtain the DEGs between cluster A and cluster B in the TCGA-OV cohort, the R package "limma" was used in the standard comparison mode. The DEG threshold was set to | logFC | > 1, and the significance criteria for identifying DEGs was set as an adjusted P value <0.05. Metascape (http:// metascape.org) is an online interactive website that helps biologists perform functional enrichment analyses on specific gene sets (19). We first introduced the 19 ARGs into Metascape and identified all the statistically enriched terms. The remaining significant terms were then hierarchically clustered into a tree based on k-statistical similarities among their gene memberships. Then, we analyzed all the transcripts based on the fold change (log2) obtained by difference analysis of the two different OC subtypes using Gene Set Enrichment Analysis (20) (GSEA, http:// software.Broadstitute.org/GSEA/) to evaluate the skewness of the two distributions of the selected genes in the list of ranked genes. Then we set the GSEA threshold for significantly enriched functional annotations to P value<0.05 and | Normalized Enrichment Score (NES) |≥1.

Identification and Validation of the Prognostic Gene Signature
The "glmnet" R package was utilized to carry out the LASSO COX regression. The "glmnet" R package was utilized to carry out the univariate Cox regression to screen 19 ARGs. Using the "glmnet" software package of R for LASSO Cox regression analysis, seven autophagy-related genes were screened to determine the best predictive model, and these genes were selected to further calculate the risk score of each patient (21,22): According to the median risk coefficient, the patients are divided into a high-risk group and a low-risk group. The Cox proportional hazard regression model includes risk score, age, grade, lymphatic invasion, and staging. The hazard ratio (HR) value distinguishes the prognostic predictors of risk genes and protective genes (HR>1 is a risk gene, HR<1 is a protective gene, p<0.05). Subsequently, Kaplan-Meier survival analysis was performed, and the sensitivity and specificity of the ROC curve were used to evaluate the prognostic performance of the signature. Circos is a visualization tool that can be used to identify and analyze the similarities and differences generated by genome comparisons, effectively display changes in genome structures, and generally display any other types of positional relationships between genome intervals (23,24). The mutations of 7 ARGs in the TCGA-OV cohort were downloaded from the cbioportal website (https://www.cbioportal.org/). We used Pearson correlation analysis to analyze the correlation between the risk score and common immune checkpoints and tried to use the risk score to accurately predict the effect of treatment.

Quantitative Real-Time PCR (qRT-PCR)
Total RNA was isolated from the OC samples using TRIzol reagent (Invitrogen). The PrimeScript RT (reverse transcription) Kit (TaKaRa Bio) was used to obtain cDNA. qRT-PCR was conducted using the AceQ qPCR SYBR Green Master Mix (Vazyme). b-actin was used as internal control.

Immunohistochemistry (IHC)
Tissues were derived from clinical specimens. The tissue sections were deparaffinized, hydrated, repaired and blocked with citric acid antigen. Subsequently, the tissue sections were probed with an ULK2 antibody (1:50, Omnimabs, OM294638) and GABARAPL1 antibody (1:200, Abcam, ab86497) at 4°C overnight. The sections were washed with PBS for 5 times, and then a secondary antibody was added and incubated at room temperature for 10 minutes. DAB and hematoxylin were added for visualization.

Statistical Analysis
All of the analyses were performed with R software (version 3.6.1, http://www.R-project.org). Univariate and multivariate Cox proportional hazard regression analyses were used to evaluate the relationship between the risk scores and OS. The area under the ROC curve (AUC) ("timeROC" package in R) was used to analyze the sensitivity and specificity of genotyping and gene signature risk scores in predicting survival rate. AUC can be used as an accuracy indicator of prognosis. All statistical P values were bilateral in all the analyses, and P < 0.05 was statistically significant. The primary prognosis endpoint was overall survival, and survival curves were estimated using the Kaplan-Meier method. The log-rank test was used to determine the significance of the difference. The "Surv-cut point" function that repeatedly tests all possible cut-off points to obtain the largest rank statistic was used to dichotomize the differential genes, and then, the largest log-rank statistic was selected to divide the patients into high and low subgroups to reduce the calculated batch effect. The Paired t-test was performed to analyze statistical significance of qRT-PCR data.

Overall Design of This Study
We have developed a flow chart to systematically describe our research ( Figure 1). The clinical data and corresponding gene expression profiles of OC patients were downloaded from the TCGA database. One patient without prognostic data was excluded. The ARG sets were downloaded from the MsigDB for follow-up analysis. The subtypes of OC were classified by kmeans unsupervised clustering and were divided into cluster A and cluster B. Then, the biological functions, metabolic pathways and signal transduction pathways with significant enrichment of differential genes were analyzed by the Metascape database, and the signaling pathways enriched in cluster A and cluster B were analyzed by GSEA. Seven ARGs were obtained through LASSO Cox regression analysis, and the risk prediction model of these ARGs was established and evaluated. The correlation between risk genes and immune cells was analyzed, and the expression of the ARGs in OC tissues and paracancerous tissues was verified by experiments.

Identification of Two Subtypes of OC
The 307 patients with OC clustered by k-means based on the mRNA expression of 19 ARGs. It was found that 2 was the optimal and stable number of clusters (Figures 2A, B). Most of the patients were enriched in cluster B, which correlated with a poor prognosis of OC ( Figure 2C). Cluster A tended to have better survival than cluster B. 535 differentially expressed genes (DEGs) between cluster A and cluster B were identified, as shown by a volcano map these genes included 108 upregulated genes (log2 fold change ≥ 1) and 427 downregulated genes (log2 fold change ≤ -1) ( Figure 2D). We analyzed the association among the 19 ARGs from the TCGA RNA-seq data and the clinicopathological features, including grade, stage, lymph node metastasis, survival status and overall survival time, of 307 OC patients ( Figure 2E).

Estimation of Cell Infiltration Into the TME
We quantified the level of immune cell infiltration to evaluate the immune landscape of cluster A and cluster B. First, we explored the relationship among cluster A, cluster B and gene expression of seven common immune checkpoints, including CD274, PDCD1LG2, CTLA4, LAG3, PDCD1, HAVCR2 and TIGIT. These genes were selected based on drug inhibitors that are currently in clinical trials or have been specifically approved in different tumor types. The results showed that the overall expression of the common immune checkpoints in cluster A was significantly higher than that in cluster B ( Figure 3A). Next, we analyzed the difference in the expression of common CD8 + T marker genes (CD8A, GZMB, CXCL9, CXCL10, PRF1, TBX21, and CD8B) in cluster A and cluster B and found that the expression of CD8 + T marker genes in cluster A was significantly higher than that in cluster B ( Figure 3B). In addition, ssGSEA was used to identify the abundance of tumor-infiltrating immune cells in cluster A and cluster B. We found that the levels of B cells, T cells, Treg cells, TFH cells, macrophages, aDCs, iDCs, CD56dim cells, Th1 cells, and Tgd cells in cluster A were significantly higher than those in cluster B, while the level of infiltrating NK cells in cluster A was significantly lower than that in cluster B ( Figure 3C). We quantified 24 kinds of immune cells, such as B cells, T cells, NK cells and macrophages, and drew a heatmap that included the relationship among age, stage, grade, lymphatic metastasis and immune infiltration. Through analyzing the heatmap, we observed that cluster A showed high immune infiltration, while cluster B, on the contrary, showed low immune infiltration ( Figure 3D). This result indicates that the immune response of cluster A is active, while that of cluster B is suppressed; thus, we speculate that OC patients in cluster A may have a better response to immunotherapy.

Identification of the Involved Signaling Pathways
DEGs were analyzed using Metascape to better understand their functional and pathway enrichment. First, all the statistical enrichment items were determined, the cumulative hypergeometric p-value and enrichment factor were calculated, and filtering was performed. Then, the remaining important terms were clustered into a tree in a hierarchical structure according to the Kappa statistical similarity between its gene members. The 0.3 kappa score was then used as a threshold to coerce the tree into term clusters. Next, we selected a subset of representative terms from this cluster and converted it into a network layout. The nodes of the same enrichment network are colored according to their p-values ( Figures 4A, B). From the bar graphs and network graphs, we found that the terms related to immune function were the most abundant; among these terms, lymphocyte activation, adaptive immune response and cytokine-mediated signaling pathway were the 3 terms with the highest significance levels. GSEA was used to execute gene ontology functioinf and pathway enrichment analysis. GSEA analysis showed that signaling pathways such as interferon gamma, interferon alpha, IL-6 JAK STAT3, Inflammatory, P53 pathway, and hypoxia were significantly enriched in the cluster A group ( Figure 4C and Table 1). Most of these pathways have been confirmed to be related to immunotherapy (25,26), suggesting that ARGs may affect tumor immune regulation and providing a certain experimental direction for further research on immunotherapy in the future.

Autophagy-Related Prognosis Classifier and Clinicopathological Characteristics of OC
We performed LASSO Cox regression analysis on 19 ARGs and further obtained 7 ARGs (ATG12, ATG4A, ATG4C, ATG5, GABARAPL1, IFNG, and ULK2) ( Figures 5A, B). A risk score to predict the prognostic value of these genes was calculated, and patients with OC were separated into the low-risk or high-risk groups. We found that with the risk score increased, the number of deaths increased ( Figure 5C). The ROC curve shows that the classifier has strong predictive ability, with an AUC of 0.63 in 3 years and 0.718 in 5 years ( Figure 5D). Univariate Cox regression showed that the risk core was a risk factor for DFS (HR>1, p<0.05), and it had a better prognostic effect than other clinical indicators ( Figure 5E). Kaplan-Meier analysis indicated that high-risk patients had significantly worse overall survival than low-risk patients ( Figure 5F). This result indicates that the risk model may serve as a promising indicator for evaluating the prognosis of OC patients and may be a powerful prognostic indicator.

Use of 7 ARGs in Survival Analysis of OC
We evaluated the prognostic values of the 7 ARGs in OC. As shown in Figure 6A,

Predicting the Effectiveness of ARGs and Correlation With Our Risk Score
The results of Pearson's correlation analysis indicated that our risk score was significantly correlated with the mRNA expression level of ARGs, namely, ATG12 (R =0.26, p =7.5e-07), GABARAPL1 (R =0.

Prognostic Signature Is Related to the TIM
The mutations in the 7 ARGs in the TCGA-OV cohort were downloaded from the cbioportal website (https://www.cbioportal. org/). The gene of highest mutation frequency is GABARAPL1, reaching 5%, of which the main type of mutation was amplification mutations. Following genes were ATG5, IFNG, and ATG4C, which were 3, 3, and 2.5%, respectively ( Figure 7A). We also found that the expression of IFNG in OC was more scattered and lower, while the expression of ATG12, ATG4A, ATG4C, ATG5, GABARAPL1, and ULK2 was relatively concentrated and higher ( Figure 7B). A correlation analysis of the 7 ARGs was performed to obtain Circos plots ( Figure 7C). We explored the correlation between 7 ARGs and immune cells ( Table 3).
Correlation analysis was also performed to analyze the 24 immune cells. The immune cell proportions were weakly to strongly correlated. T cells and cytotoxic cells showed the strongest positive correlation; macrophages also indicated a strongly positive correlation with iDCs. Eosinophils and Th17 cells showed a moderate negative correlation ( Figure 7D).
Additionally, we analyzed the correlation between 7 immune genes and immune cells and found that ATG5, ATG4C, ATG12, ULK2, and GABARAPL1 were positively correlated with NK cells,

Validation of the Gene Signature in Clinical Tissue Samples
To confirm the reliability of the identified gene signature, we examined the ULK2 and GABARAPL1 expression levels by qRT-PCR, WB, and IHC using 8 pairs of OC tumor tissues and paracancerous tissues. The results showed that the ULK2 and GABARAPL1 mRNA in the tumor tissues were significantly downregulated compared with those in the paracancerous tissues (p<0.05) ( Figure 8A). Moreover, the protein level of ULK2 and GABARAPL1 were detected using IHC in 8 OC clinical tissues that underwent cytoreductive surgery. The IHC analysis showed that ULK2 and GABARAPL1 were lowly expressed in the OC tissues ( Figures 8B, C).

DISCUSSION
The early symptoms of patients with OC can't be easily detectable. Due to difficulty in diagnosis at early stage, tumor recurrence, chemoresistance, and poor prognosis, improving clinical treatment is very challenging. Therefore, it is urgent to identify reliable biomarkers for early prognosis to facilitate the diagnosis and treatment of OC. Induction of autophagy has shown promise in the management of a wide range of illnesses, including neurodegenerative disorders, cardiovascular diseases, and rheumatic diseases (27)(28)(29). Accumulating evidence suggests that autophagy plays important roles in tumor chemotherapy and radiotherapy. In addition, autophagy affects TIM and tumor immunotherapy. Proper regulation of autophagy could enhance the immune response, leading to immunotherapy potentiation. Recently, a medical breakthrough was made in the field of cancer immunotherapy. However, how the immune system uses immune cells to eradicate tumors in a complex environment remains to be further investigated. Unfortunately, only a small proportion of cancer patients and only some cancers currently respond to immunotherapy (30). How to make more patients benefit from personalized immunotherapy and improve the efficacy of immunotherapy is the focus of tumor immunology. Immune cells play a protective role in the tumor microenvironment, which can lead to immune checkpoint inhibition and immune cell infiltration (31). Therapeutic strategies that target the immune microenvironment, such as immunosuppressive cells and immunosuppressive factors, can effectively prevent tumor cells from escaping immune surveillance. Therapies that block immune checkpoints show longer-lasting responses than traditional chemotherapy and have been approved by the FDA for the treatment of multiple cancer types, such as melanoma, nonsquamous cell carcinoma, OC, etc. (32)(33)(34). However, therapies that block immune checkpoints have low response rates in approximately 10%-30% cancers, which may be related to tumor mutational burden, PD-L1 expression level, IFN signaling and MHC-I loss. However, the current understanding of the low response rate of immunotherapy is still limited, and the regulatory mechanism is not clear.
Targeting the tumor microenvironment by immunotherapy may be a promising therapeutic strategy. Autophagy regulates innate and adaptive immune responses by regulating immune components, such as NK cells, T lymphocytes, B cells and DC cells (35). It was reported that autophagy could promote antigen presentation by MHCII and MHCI molecules on APCs. Autophagy provides an important antigen source for the loading of MHCII, thus activating antigen-specific CD8 + T cells (36). In addition, autophagy not only regulates cytokines to affect immunity but also increases the response of MHCI molecules to IFN-g to enhance cross-presentation (37). The defection of autophagy causes the accumulation of intracellular fat droplets, promotes the release of linoleic acid, and leads to the exhaustion of liver CD411 cells, thereby inducing an immunosuppressive tumor microenvironment and promoting tumor progression (38,39). Antibody-mediated responses to both T cell-dependent and T cell-independent antigens also require autophagy. Otherwise, these responses will cause endoplasmic reticulum stress and T cell death. Autophagy-related gene suppression in human pancreatic ductal adenocarcinoma 8988T cells can promote the expression of PD-L1, which is conducive to the establishment of an immunosuppressive tumor microenvironment (40).
There have been no comprehensive and systematic studies on the relationship among autophagy, the immune microenvironment and immunotherapy in OC. This study aims to identify suitable immunotherapy targets and effective prognostic biomarkers and to predict the response and prognosis of OC. Using k-means clustering, we successfully divided OC samples into two groups, cluster A and cluster B. Significant differences were found between cluster A and cluster B. Cluster A had a better prognosis for OS and had a higher immune infiltration level than cluster B. Next, we use Metascape to analyze the differential mRNA (FDR<0.05) of cluster A and cluster B for signal pathway enrichment analysis. Lymphocyte activation, adaptive immune response and cytokinemediated signaling pathway were the top three signaling pathways, and all of these pathways were closely related to the function of the immune system. We then used GESA for pathway enrichment and found that IFN-g, IFN-a, IL-6 JAK STAT3, inflammatory, P53 pathway, and hypoxia were significantly enriched in cluster A. Most of these pathways are closely related to tumorigenesis and the TME. For example, interferon plays a vital role in the immune response and has been studied in a variety of malignant tumors (41,42). IFNg signal transduction can promote Treg function in autoimmunity, and activation of the IFN-a signaling pathway leads to a more effective antiviral response and enhanced antitumor immunity (43). The JAK/STAT3 pathway is aberrantly activated in various cancers, including OC, and is involved in functional regulation of the tumor microenvironment (44). Inflammatory factors could induce immune suppression and mediate immune escape (45). The P53 tumor suppressor pathway plays a key role in tumor immunology and the homeostatic regulation of immune responses (46). Our classification is of great significance in guiding clinical work. Patients in cluster A have a better prognosis and significantly higher immune infiltration than those in cluster B, and the enriched signaling pathways are mainly related to the TIM. Thus, cluster A may obtain better clinical outcomes with immunotherapy and have a certain clinical reference value for OC.
Gene markers are widely used in modern clinical diagnosis. Traditional markers, such as CA125, CA199 and alpha-fetoprotein (AFP), are frequently used to diagnose OC, gastric cancer and liver cancer, predicting prognosis, and monitoring postoperative recurrence, respectively (47)(48)(49). Additional searches for more sensitive and specific markers are needed to improve clinical diagnosis and treatment. We established a 7 ARG prognostic risk signature by LASSO Cox regression and then divided patients into high-risk and low-risk groups. ROC curve and Kaplan-Meier survival analyses showed that the model was reliable and helpful to identify high-risk and low-risk patients. Univariate Cox regression analyses showed that the risk score was a more reliable indicator of DFS than other factors. The risk score was also closely related to risk genes and immune checkpoint molecules. Some risk genes have been confirmed to be closely related to the formation and development of OC. IFNG expression is related to the upregulation of XBP1, and lacking XBP1 selectively in T cells could promote antitumor immunity (50). P53 and RAS mutants regulate apoptosis and autophagy through the dysregulation of ATG12 and affect chemotherapy resistance (51). It was reported that ATG5 could promote the growth of OC cells in the peritoneal microenvironment, and inhibition of ATG5-induced autophagy sensitizes OC cells to olaparib and other PARP inhibitors (52,53). However, GABARAPL1 and ULK2 have not been reported in OC. Experimental verification found that ULK2 and GABARAPL1 are downregulated in OC, combined with OS analysis, the results show the higher expression of ULK2 indicated a better prognosis whereas the higher expression of GABARAPL1 showed a worse prognosis in OC. ULK2 may be a protective prognostic factor. Unexpectedly, GABARAPL1 may be a risk prognostic factor, the reason may be that the sample size was  too small for inferential statistics or that the characteristics of the autophagy gene itself caused this result, but the specific reasons need to be further explored. ULK2 and GABARAPL1 may be potential therapeutic targets for OC. The successful classification of OC in this study is conducive to screening patients suitable for immunotherapy, and 7 marker genes can be used as strong biological markers for the prognosis of OC. The research, however, has some limitations. Firstly, the number of ARG sets studied is relatively small, and there is no validation set. In addition, we verified the expression of some ARGs in tissues, but the biological functions of the ARGs need further experimental research.
In summary, the autophagy-immune-based gene risk signature might be helpful to guide the clinical treatment, evaluate the prognosis and predict the efficacy of immunotherapy. ARGs might be potential markers for the prognosis of OC immunotherapy.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The IRB of Third Xiangya Hospital, Central South University. The patients/participants provided their written informed consent to participate in this study.