Immunotherapeutic Value of MAP1LC3C and Its Candidate FDA-Approved Drugs Identified by Pan-Cancer Analysis, Virtual Screening and Sensitivity Analysis

Background: The autophagy pathway within the tumour microenvironment can be regulated to inhibit or promote tumour development. In the fight against tumour growth, immunotherapy induces an anti-tumour immune response, whereas autophagy modulates this immune response. A key protein in the autophagy pathway, microtubule-associated protein 1 light chain 3 (MAP1LC3), has recently become a hotspot for tumour research. As a relatively novel member, the function of MAP1LC3C in tumours still need to be investigated. Therefore, the goal of this study was to look into the possible link between MAP1LC3C and immunotherapy for 33 kinds of human malignancies by using pan-cancer analysis. Methods: High-throughput sequencing data from The Cancer Genome Atlas, Genotype-Tissue Expression Project and Cancer Cell Line Encyclopedia databases, combined with clinical data, were used to analyze the expression of MAP1LC3C in 33 types of cancer, as well as patient prognosis and neoplasm staging. Activity scores were calculated using ssGSEA to assess the MAP1LC3C activity in pan-cancer. Associations between MAP1LC3C and the tumour microenvironment, including immune cell infiltration and immunomodulators, were analyzed. Moreover, tumour tissue ImmuneScores and StromalScores were analyzed using the ESTIMATE algorithm. Additionally, associations between MAP1LC3C and tumour mutational burden/microsatellite instability, were investigated. Finally, based on the expression and structure of MAP1LC3C, the United States Food and Drug Administration (FDA)-approved drugs, were screened by virtual screening, molecular docking and NCI-60 drug sensitivity analysis. Results: Our study found that MAP1LC3C was differentially expressed in tumour and normal tissues in 23 of 33 human cancer types, among which MAP1LC3C had prognostic effects in 12 cancer types, and MAP1LC3C expression was significantly correlated with tumour stage in four cancer types. In addition, MAP1LC3C activity in 14 cancer types was consistent with changes in transcription levels. Moreover, MAP1LC3C strongly correlated with immune infiltration, immune modulators and immune markers. Finally, a number of FDA-approved drugs were identified via virtual screening and drug sensitivity analysis. Conclusion: Our study investigated the prognostic and immunotherapeutic value of MAP1LC3C in 33 types of cancer, and several FDA-approved drugs were identified to be highly related to MAP1LC3C and can be potential cancer therapeutic candidates.


INTRODUCTION
Autophagy is a process that occurs in all eukaryotes, involving in the capture by an autophagosome and transport to the lysosomes for decomposition and recycling (Ohsumi, 2012;Morishita and Mizushima, 2019). Autophagy has long been considered as a nonselective process that meets the needs of cell synthesis and metabolism; however, recent studies have proved the existence of selective autophagy pathways, specifically targeting damaged organelles, pathogens and unfolded proteins (Le Guerroué et al., 2017). To better cope with stresses in the tumour microenvironment, the autophagy pathways in various cell types can be regulated to inhibit or promote tumour development (Xia et al., 2021).
The immune system plays an important role in preventing tumour occurrence, development and metastasis and in tumour treatment. Immune surveillance is a function of the immune system that helps to recognise, kill and remove tumour cells; however, tumour cells can evade tumour immunity through immunosuppressive responses (Zou, 2005). Moreover, immunotherapy fights against tumours by arousing the antitumour immune response in the immune system (Zou et al., 2016). Recent studies have shown that autophagy is involved in the various biological processes of immune cells (Clarke and Simon, 2019) and can modulate tumour growth by regulating the immune response (Yang et al., 2018).
A common key factor in both selective and non-selective autophagy pathways is the ubiquitin-like protein binding system, including the ATG5-ATG12 binding system and ATG8-lipidation system (Tooze and Yoshimori, 2010). Human cells contain at least six ATG8 family members, which can be divided into two subfamilies: MAP1LC3A (LC3A), LC3B, LC3C and GABARAP, GABARAPL1, GABARAPL2 (Slobodkin and Elazar, 2013). Microtubule-associated protein 1 light chain 3 (MAP1LC3), a ubiquitinated protein of the ATG8-lipidation system, is involved in various pathophysiological processes of the body. The role and mechanism of MAP1LC3 subfamily in a tumour is complex, making it a hotspot for tumour research. However, research on MAP1LC3C, a relatively new member, is still lacking. There are few articles referring to MAP1LC3C and no articles that deeply investigate the function of this gene in tumours. Therefore, a pan-cancer study of MAP1LC3C is necessary.
This study is the first to analyse the prognostic and immunotherapeutic value of MAP1LC3C in various cancer types. The expression and activity of MAP1LC3C and its correlation with patient prognosis and tumour stage in 33 human cancer types, combined with clinical information, were analysed using online databases such as TCGA, GTEx and CCLE. The effects of MAP1LC3C on the immune microenvironment, tumour mutational burden (TMB), microsatellite instability (MSI) and two immunotherapy biomarkers were studied using CIBERSORT and ESTIMATE algorithms. Finally, based on the expression and structure of MAP1LC3C, virtual screening and drug sensitivity analysis of the FDA-approved drug library were conducted to identify therapeutic drugs. Therefore, this study aims to elucidate the role of MAP1LC3C in the immune system and its impact on cancer prognosis and immunotherapy, as well as to provide some references for therapeutic agents.

Data Collection
Transcriptome expression profiles, clinical data and mutation data of 33 cancer types in the TCGA database were downloaded through Xena (Goldman et al., 2020). Control gene data of six cancer types (ACC, LGG, LAML, OV, TGCT and UCS) were obtained from the GTEx database (Consortium, 2020). Sequencing data of cell lines were obtained from the CCLE database (Nusinow et al., 2020). As this study used open data, ethical approval was not required.

Identification of MAP1LC3C Differential Expression and Its Association With Clinical Characteristics
MAP1LC3C expression data were extracted from the TCGA and GTEx databases, and the gene expression levels between 33 cancer and normal samples were compared using R-package LIMMA (Ritchie et al., 2015). Additionally, the correlation between MAP1LC3C expression and tumour stage was explored. Combined with clinical information, the relationship between gene expression level and patient survival was calculated using univariate Cox regression. When the hazard ratio (HR) was greater than 1 (HR > 1), MAP1LC3C expression was considered as a risk factor. Kaplan-Meier (KM) analysis was performed to compare the overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI) and progression-free interval (PFI) of patients with cancer, which was stratified by median MAP1LC3C expression. p < 0.05 was considered statistically significant for the analyses.

Study of MAP1LC3C Activity
To further study the activity of MAP1LC3C in pan-cancer, genes that are closely related to MAP1LC3C as relevant genes of MAP1LC3C through the co-expression method were found, and the MAP1LC3C activity score of each sample was obtained using R-package GSVA and ssGSEA method (Hänzelmann et al., 2013). The differences in MAP1LC3C activity between the normal and tumour groups were studied, and the scores of MAP1LC3C activity in 33 cancer types were obtained.

Correlation Analysis of MAP1LC3C Expression and Immune-Related Factors
The CIBERSORT algorithm was used to calculate the immune cell invasion level of each tumour sample (Newman et al., 2019). A total of 100 permutations were run using the LM22 signature. p < 0.01 and |R| > 0.4 indicated significant correlation. Immunoscores and stromal scores were calculated for each sample using the ESTIMATE algorithm (Yoshihara et al., 2013). Additionally, the list of immune modulators (Supplementary Table S1), including immune inhibitors, stimulators and major histocompatibility complex (MHC) molecules, were downloaded from the TISIDB database (Ru et al., 2019).

VIRTUAL SCREENING AND MOLECULAR DOCKING
Structural information on 2858 FDA-approved and pharmacopeial drugs was downloaded from the TargetMol (Inc.A-Approved &, 2021). The spatial structure information of the MAP1LC3C protein (PDB 2NCN) was downloaded from the Protein Data Bank (PDB) database (Berman et al., 2000). The GHECOM algorithm was used to identify potential small molecule binding sites on the protein (Kawabata, 2010), and a docking pocket with volume 2,571 Å 3 was defined. UCSF DOCK 6.9 was used for virtual screening and molecular docking. Finally, PyMol was used to visualise the docked conformation, and Ligplus was used to analyse the interaction force (Wallace et al., 1995).

DRUG SENSITIVE ANALYSIS
The 23,808 identified RNA expression data and 23,255 analyzed drug data from the NCI-60 cell line were downloaded from the CellMiner database (Reinhold et al., 2012). To ensure the clinical practicality of the analysis results, FDA-approved drugs that had undergone clinical trial were selected to obtain a total of 792 drugs for the screening. MAP1LC3C expression was extracted and drug sensitivity (IC50) was calculated to obtain Pearson correlation coefficients between gene expression and different drugs, and the results were screened and visualized according to p < 0.05.

Differential Expression of MAP1LC3C in 33 Cancer Types
Detailed information on the 33 types of cancer included in this study is shown in Supplementary Table S2. The gene expression profiles of various tumour cell lines downloaded from the CCLE database, and the MAP1LC3C expression levels of 21 tissues according to their tissue sources were analyzed. MAP1LC3C showed inconsistent expression levels in different cell lines (p = 3.8E-10, Figure 1A), but in soft tissue tumour cells showed the highest expression levels. Considering that several tumours in the TCGA database do not include normal sample data, normal sample data from the GTEx database were integrated with the tumour sample data from the TCGA database to analyse the expression differences of MAP1LC3C in 33 tumour types ( Figure 1B). MAP1LC3C was significantly up-regulated in ACC, GBM, KIRC, KIRP, LAML, LGG, TGCT   Figure 3, TCGA data were used to analyse MAP1LC3C activity in 33 tumour types. The results showed that the transcription level was matched with MAP1LC3C activity. MAP1LC3C activity increased significantly in GBM, KIRC, BLCA, BRCA, CESC, COAD, HNSC, KICH, LIHC and LUAD tumour types, but decreased significantly in LUSC, PCPG, PRAD, READ, STAD, THCA and UCEC tumour types ( Figure 3A). 4 tumour types including UVM, MESO, DLBC and GBM showed relatively high activity ( Figure 3B).

PROGNOSTIC ROLE OF MAP1LC3C EXPRESSION
For the prognostic analysis, we selected clinical indicators, including overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI) and progression-free interval (PFI). OS was defined as the time from the date of diagnosis to death, regardless of the cause. In OS analysis, Univariate Cox regression identified high MAP1LC3C expression as a risk factor for COAD, KICH, KIRC, KIRP, LGG, LUSC and UCEC, but as a protective factor for UVM ( Figure 4A). KM analysis revealed that patients with high MAP1LC3C expression in LGG, LUSC, STAD and UCEC cancer types had lower OS rates than those with low MAP1LC3C expression. However, those with high MAP1LC3C expression in LAML and UVM had higher OS rates ( Figures  4B-G). In DSS analysis, unlike OS, patients who died from causes other than the specified disease were not counted. Univariate Cox regression indicated high MAP1LC3C expression as a risk factor for BRCA, COAD, KIRC, KIRP, LGG and UCEC, but as a protective factor for UVM ( Figure 5A). KM analysis demonstrated that patients with high MAP1LC3C expression in LGG and UCEC had lower DSS rates than those with low MAP1LC3C expression, while those with high MAP1LC3C expression in LUAD and UVM had higher DSS rates ( Figures 5B-E). In DFI analysis, patients who died from causes other than the specified disease were not counted. Univariate Cox regression identified high MAP1LC3C expression as a risk factor for KIRP, LGG, STAD and UCEC ( Figure 6A). KM analysis showed that patients with high MAP1LC3C expression in ESCA, STAD and UCEC had lower DFI rates than those with low MAP1LC3C expression ( Figures 6B-D). Unlike DFI, PFI was defined as progression or death from disease, again from any cause. In PFI analysis, Univariate Cox regression identified high MAP1LC3C expression as a risk factor for KIRC, KIRP, LGG, PCPG, STAD and UCEC, but as a protective factor for CHOL and UVM ( Figure 7A). KM analysis confirmed that patients with high MAP1LC3C expression in LGG, STAD and UCEC had lower PFI rates than those with low MAP1LC3C expression ( Figures 7B-D). Notably, KIRP, LGG and UCEC showed significant differences in all four analyses above, and high MAP1LC3C expression suggested a poor prognosis in all three types of cancers. Correlation Between MAP1LC3C Expression and Tumour Immunity Using CIBERSORT, the detailed immune cell composition of all patients in the TCGA database was calculated, and the correlation between 22 immune cells in 33 cancer types and MAP1LC3C expression was determined (Supplementary Table S3). The results revealed that a majority of immune cells were significantly correlated with MAP1LC3C expression. As shown in Figure 8, three types of immune cells in TGCT, two types in PAAD and one type in DLBC and BRCA were correlated with the expression of MAP1LC3C, respectively (p < 0.01 and |R| > 0.4). As shown in Figure 9A, 23 immune inhibitors were analyzed, and the expression of MAP1LC3C was significantly correlated with the immune inhibitors of various cancer types, including CHOL, LGG and PAAD. Correlation analysis of 45 immune stimulators revealed a significantly positive correlation between MAP1LC3C expression in CHOL and PAAD and many immune stimulators in Figure 9B. Additionally, correlation analysis of 21 MHCs revealed a significantly positive correlation between MAP1LC3C expression in CHOL, LGG, LIHC and LUSC and multiple MHC molecules in Figure 9C.
ESTIMATE algorithm was used to calculate the tumour tissue's immuno/stromal-scores and evaluate the relationship between the immuno/stromal-scores and MAP1LC3C expression. Figure 9D showed a significantly positive correlation of MAP1LC3C expression in PAAD, LUSC, CHOL, BLCA, LIHC, ESCA, BRCA, HNSC, PCPG and READ with immuno/stromal-scores (p < 0.01 and |R| > 0.5). The detailed results of ESTIMATE scores are summarized in Supplementary Table S4.

BIOLOGICAL FUNCTION OF MAP1LC3C
Given the differential expression of MAP1LC3C in many cancers, GSEA was used to assess the biological function of MAP1LC3C in 33 cancer types. In KIRC, MAP1LC3C shows enrichment in the following GO terms: GOBP_ENDOTHELIAL_CELL_MIGRATION, GOBP_REGULATION_OF_ENDOTHELIAL_CELL_MIGRATION, GOBP_RIBONUCLEOPROTEIN_COMPLEX_BIOGENESIS, GOBP_RIBONUCLEOPROTEIN_COMPLEX_BIOGENESIS and GOMF_RNA_BINDING_INVOLVED_IN_POSTTRANS CRIPTIONAL_GENE_SILENCING, and in the following  Figures 10D, 11F).

VIRTUAL SCREENING AND MOLECULAR DOCKING BASED ON MAP1LC3C STRUCTURE
Using virtual screening and molecular docking of the FDA-Approved and Pharmacopeia Drug Library, the top 10 drugs with the best docking score were obtained (  Figure 11A (Shears, 2001). Phytic acid is being investigated in clinical trials NCT01000233 for its value in the prevention of cardiovascular calcifications. Ceftaroline fosamil is a cephalosporin antibacterial agent for the treatment of the following infections caused by specified susceptible bacteria: Acute bacterial skin and skin structure infections and community-acquired bacterial pneumonia (Steed and Rybak, 2010). In summary, all of these drugs play an important role in the treatment of non-oncological diseases, but our studies suggest that they also have potential therapeutic value in oncology.

Correlation Between MAP1LC3C Expression and Drug Sensitivity
The examination of MAP1LC3C expression in NCI-60 cell lines showed a significant correlation between MAP1LC3C expression and drug sensitivity ( Table 2). In particular, as MAP1LC3C expression increased, 6-Thioguanine, By-Product of CUDC-305, 8-Chloro-adenosine, DIGOXIN, AT-13387 and Volasertib had a lower IC50 against cancer cells (Figures 12,  13). This implies that increased expression of MAP1LC3C

DISCUSSION
This study aimed to comprehensively investigate the differential expression of MAP1LC3C between the normal and tumour  LGG, LUAD, LUSC, PCPG and STAD. Additionally, MAP1LC3C expression was correlated highly with the tumour microenvironment, immune cells, immune modulators, TMB and MSI. The biological functions and signalling pathways related to MAP1LC3C expression were also identified. Finally, FDA-approved drugs were screened through virtual screening and drug sensitivity analysis. Immunotherapy represents a shift in treatment modalities for oncology, the goal is no longer to target the tumour itself, but to overcome the immune suppression caused by the tumour and its microenvironment, and then allow the immune system to target and kill the cancer cells (Billan et al., 2020). About 10 years ago, the first immune checkpoint inhibitor, ipilimumab, a monoclonal antibody targeting CTLA-4, was approved for the FDA and it marked the beginning of the immunotherapy era (Squibb, 2020a). Two additional immune checkpoint inhibitors targeting PD-1 (pembrolizumab and nivolumab) were subsequently approved, and all three were initially approved for unresectable or metastatic melanoma (Squibb, 2020b;Merck, 2020). More immune checkpoint inhibitors have subsequently been approved for a variety of cancers. As a result, immunotherapy is increasingly being used in the clinical management of tumours and our study also focuses on the analysis of the value of MAP1LC3C in immunotherapy to provide new insights into future immunotherapy regimens.
Identifying abnormally expressed genes and tumour-specific targets or features for personalized treatment in different cancers could increase the possibility of cure or remission for patients (Andre et al., 2014). Therefore, the use of TCGA and GTEx databases for pan-cancer analysis helps to determine the differential expression and role of MAP1LC3C in various types of cancer (Cao and Zhang, 2016;Cava et al., 2018). Furthermore, a thorough pan-cancer analysis can be performed in cell lines using the CCLE database to assess gene expression, which may have implications for future cell experiments. Consistent with previously reported results, differentially expressed MAP1LC3C has certain prognostic values in some cancers, especially COAD, LGG and LUAD. MAP1LC3C is a marker that indicates poor prognosis in patients with colon cancer, low-grade glioma and lung adenocarcinoma (Xu et al., 2020;Guo et al., 2021;Wang et al., 2021). Generally, the protein expression level better reflects the tissue activity (Mo et al., 2021). Due to the lack of relevant data on protein expression levels in public databases, it is difficult to study MAP1LC3C protein expression level. However, MAP1LC3C activity score in various cancer types was obtained using the ssGSEA method. By comparing the transcription level and activity score, the transcription levels of some cancers (BLCA, BRCA, CESC, COAD, GBM, HNSC, KICH, KIRC, LUAD, LUSC, PCPG, READ, STAD and THCA) were matched with MAP1LC3C activity, indicating that the transcription levels represented activation or repression of MAP1LC3C. If the transcription levels could not be matched with MAP1LC3C activity in certain cancers, it indicates that it may be due to post-transcriptional protein level modifications or protein metabolism, which affect MAP1LC3C expression. In our study, there were no cancers with mismatches.
Tumour tissue contains not only tumour cells but also immune cells. Immune cells that infiltrate the tumours can profoundly influence tumour development and anti-cancer therapy (Dalton et al., 1993). Therefore, the quantification of immune cells has extraordinary significance. Recently, immunotherapy has shown an increased efficacy in the treatment of tumours (Finn, 2008). This study reports that the expression level of MAP1LC3C is related to cancer immunity. A strong correlation between MAP1LC3C and macrophages M2, B cells naive and plasma cells was observed in TGCT. In PAAD, a strong correlation was noted between MAP1LC3C and macrophages M0 and B cells naive. Correlation analysis of 23 immune inhibitors showed that MAP1LC3C expression was significantly correlated with various immune inhibitors agents of different cancer types, especially CHOL, LGG and PAAD. Correlation analysis of 45 immune activators showed a significantly positive correlation of MAP1LC3C expression in CHOL and PAAD with many immune activators. The correlation between MAP1LC3C expression and 21 MHCs was also analysed. Human leukocyte antigen (HLA) is the expression product of human MHC, which is the most complex polymorphic system in the human body (Norman et al., 2017). It is worth noting that MHC is closely related to human immune response, immune regulation and few pathological state generation (Xu et al., 2006;Wang et al., 2013). The results indicated that most cancers were positively correlated with HLA. ESTIMATE is a tool for analysing tumour purity and stromal and immune cells' presence in the tumour (Yoshihara et al., 2013). ESTIMATE algorithm generates four final scores: stromal score (indicating the presence of stromal cells in the tumour tissue), immune score (indicating the invasion of immune cells in the tumour tissue), ESTIMATE score and tumour purity score. Our results revealed that MAP1LC3C expression in PAAD, LUSC, CHOL and BLCA had a significantly positive correlation with matrix fraction and immune fraction. In ESCA and READ, MAP1LC3C expression was significantly positively correlated with immune score, and in LIHC, BRCA, HNSC and PCPG, MAP1LC3C expression was significantly positively correlated with stromal scores.
Gene mutation is the main cause of cancer (Martincorena and Campbell, 2015). Specific gene mutations can be used to predict patient prognosis and treatment outcome (Sanz-Garcia et al., 2017). There is a likelihood that more neoantigens are formed with more somatic mutations in a tumour, which can help the adaptive immune system in recognising and detecting cancer. TMB provides a useful estimate of the tumour-neoantigen load , and the level of TMB affects the production of immunogenic peptides, thereby influencing the response of patients to immune checkpoint inhibitors (Havel et al., 2019). Additionally, MSI is a powerful mutant phenotype caused by DNA mismatch repair defects (Yamamoto and Imai, 2019), and it is an important indicator for predicting tumour occurrence and development (Li et al., 2020). MSI is also used as an FDAapproved biomarker for immunotherapy (Lemery et al., 2017). Our study reports that MAP1LC3C is highly negatively correlated with these two immunotherapy biomarkers (TMB and MSI) in a few cancers, such as DLBC, HNSC, LIHC, LUAD, LUSC and STAD. These results suggest that MAP1LC3C may influence the immunotherapy response of these six cancer types.
Furthermore, virtual screening is a commonly used computational technique for drug designing. Virtual screening can be divided into two categories: structure-based virtual screening and ligand-based virtual screening (Lavecchia and Di Giovanni, 2013;Forli, 2015). Virtual screening based on the three-dimensional structure of the receptor (target protein) FIGURE 13 | Boxplot of drug sensitivity in high and low expression groups of MAP1LC3C (*p < 0.05).
Frontiers in Pharmacology | www.frontiersin.org March 2022 | Volume 13 | Article 863856 15 was adopted to screen for drugs that had a good interaction with the MAP1LC3C protein from the FDA-Approved and Pharmacopeia Drug Library. Furthermore, a drug sensitivity analysis using data from the NCI-60 cell line was performed. On the one hand, these results demonstrate the feasibility of MAP1LC3C as a drug target and, on the other hand, provide a reference for the development of therapeutic drug regimens. Notably, the drug with better results in both the virtual screen and the drug sensitivity analysis is Volasertib, a Plk1 inhibitor, that has reached phase III clinical trials for adult acute myeloid leukaemia patients ineligible for intensive remission induction therapy (Döhner et al., 2016). Given the better results of Volasertib in both drug sensitivity analysis and molecular docking against MAP1LC3C, it is worth trying to explore its therapeutic value in other cancers. In general, all the screened drugs can be considered as potential therapeutic agents for various cancer types. However, further in vivo studies need to be conducted.

CONCLUSION
To the best of our knowledge, this study is the first report to investigate the prognostic and therapeutic value of MAP1LC3C in 33 types of cancer. Our results suggested that MAP1LC3C can be a valuable prognostic biomarker for certain types of cancer, and showed high correlation with important immunological indexes in certain cancers. This will facilitate us to understand the role of MAP1LC3C in the immune system and lay a solid theoretical foundation for future immunotherapy. The FDA-approved drugs identified using virtual screening and drug sensitivity analysis could be potential cancer therapeutic agents, thereby paving the way for future cancer treatment research. The bioinformatics method was used in this study to provide relatively preliminary results and in future in-depth studies are needed to clarify the association between MAP1LC3C and cancer treatment.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
Conceptualization, YB. and XZ methodology, XZ software, XZ validation, YB XS. and QZ formal analysis, XZ investigation, KL, SL, SZ, TL, LL, NB and SH; resources, XZ data curation, XZ original draft preparation, XZ review and editing, XZ, YB and XS visualization, XZ supervision, XZ, YB and XS project administration, XZ, YB and XS funding acquisition, YB and XS All authors have read and agreed to the published version of the manuscript.