Autophagy regulation patterns characterization identifies immune phenotypes and immunotherapy response in head and neck squamous cell carcinoma (HNSCC)

doi:10.21203/rs.3.rs-91606/v1

Download PDF

Primary research

Autophagy regulation patterns characterization identifies immune phenotypes and immunotherapy response in head and neck squamous cell carcinoma (HNSCC)

https://doi.org/10.21203/rs.3.rs-91606/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Autophagy degraded and recycled cytoplasmic components to maintain cellular homeostasis under stress conditions, which was recognized as double-edged sword in oncogenesis and novel target in cancer treatment. However, comprehensive analysis of the relationship between autophagy regulation and immunity has not been reported yet.

Methods

Unsupervised consensus clustering algorithm was used to identify autophagy regulation patterns. LASSO cox regression algorithm was used to build a scoring system (ATGscore) to represent the individual autophagy regulation pattern. Then integrated analysis of autophagy regulation patterns and ATGscore was performed.

Results

We have successful depicted five autophagy regulation patterns and established a scoring system (ATGscore) to represent it, which was shown to be significantly correlated with TIME infiltration, immune phenotypes, molecular subtypes, and genetic variation, etc. in 1165 head and neck squamous cell carcinoma (HNSCC) patients. Moreover, ATGscore was an independent prognostic factor and potent predictor for clinical response to immune-checkpoint inhibitors (ICIs) targeting immunotherapy.

Conclusion

Understanding the molecular characteristics of autophagy regulation patterns in HNSCC could help us to depict the underlying mechanism of tumour immunity and lay a solid foundation on combination of autophagy targeting therapies and immunotherapies for clinical application in HNSCC.

Cancer Biology

head and neck squamous cell carcinoma

tumour microenvironment

immune phenotype

tumour mutation burden

immunotherapy

Head and neck squamous cell carcinoma is one of the most common cancers with nearly 500,000 new cases and 350,000 deaths, which ranks 11th cancer-related deaths worldwide each year [1]. HNSCC is heterogeneous disease which not only originates from the different sites of the upper aero-digestive tract including oral cavity, oropharynx, larynx or hypopharynx, but also caused by distinct aetiologies. Smoking and alcohol abuse were recognized as the classical risk factors for HNSCC. Currently, infection with high-risk human papillomaviruses (HPVs) was considered as another vital cause for these tumours, especially in the oropharynx and occurring particularly in the Western world [2]. The patients with HPV-negative (HPV-) and HPV-positive (HPV+) diseases displayed different socioeconomic profiles, clinical presentation, molecular profiles, as well as the prognosis of the tumours [3–6]. Although combination therapy improved the prognosis of HNSCC patients, nevertheless, the five-year survival rate of patients with HNSCC is still fewer than 50% [7, 8]. Thus, lack of rapidly improvement in patient survival and personalized treatment approaches prompted us to further study the molecular landscape of HNSCC.

Autophagy is a constitutively conserved physiologic catabolic process to harvest energy and nutrients from cellular components by degrading and recycling damaged organelles including macromolecules and organelles [9–11]. Under normal conditions, autophagy actives at basal rates to maintain cellular homeostasis, but could respond to survival demands to induce nonselective bulk degradation or selectively target cytoplasmic constituents via cargo-specific autophagy receptors [12, 13]. Thus, autophagy was involved in numerous biological and pathological processes, including neurodegenerative diseases, aging and cancer [14]. However, autophagy plays a dual role in oncogenesis through a context-dependent manner [15]. Under different tumour type, pathological stage, genetic background, and even therapeutic regimen, distinct autophagy patterns will be sustained to either promote or inhibit tumorigenesis [16]. Lack of autophagy could trigger the accumulation of genotoxic cellular waste and induce the genetic and chromosomal alteration, subsequently facilitating precancerous cells transformation and mature cancer cell formation [17]. However, excessive autophagy continuously recycle the remodelling components and refuel energy supply, which allow cancer cells to escape the damage from the immune system or targeting drugs and promoting tumour progression [18, 19]. All of these indicated the autophagy is double-edged sword in tumorigeneses, which has been already verified in preclinical cancer models [20]. In head and neck squamous cell cancers, some reported that increased levels of cytoplasmic p62, which indicated autophagy inhibition, was correlated with reduced overall and disease specific survival [21]. While other groups demonstrated that poor clinical outcome patients exhibited higher levels of LC3-II, which suggested the re-activation of autophagy [22]. However, Zhou et al showed that radiation-induced autophagy could increase the survival of CNE-2 cells, which was further counteracted by autophagy inhibition with chloroquine, resulting in enhanced cell death [23]. We noticed that in HNSCC there also existed an autophagy paradox where played a controversial role in tumorigenesis, but the underlying mechanism has not been reported.

Currently, we are experiencing a strikingly shift from chemotherapy and radiation anti-tumour broadly towards immune-checkpoints inhibitors (ICIs) targeting immunotherapy more precisely that regulate immune response against tumours. Recent studies have reported that autophagy significantly controls the immune response. Autophagy could not only activate receptors such as toll-like receptors (TLRs) and nucleotide oligomerization domain (NOD)-like receptors (NLRs) to induce natural killer T (NKT) activation, cytokine production, and phagocytosis in innate immunity, but also provide an abundance of antigen for MHC class II molecules including HLA molecular toward dendritic cells for cross-priming to CD8 + T cells in adaptive immunity [24, 25]. Moreover, autophagy could facilitate, promote or inhibit the proliferation and differentiation of a variety of immune cells or secretion a wide range of cytokines to modulate the tumour immune microenvironment (TIME) homeostasis. Conversely, certain cytokines and immune cells also exhibit a significant impact on autophagy function [26, 27]. Furthermore, immune-checkpoints including Indoleamine 2,3 dioxygenase(IDO), CTLA-4 and PD-1 were shown to regulate tumour immune tolerance through autophagy pathways. Accumulating evidences have showed that autophagy could interact with TIME immune cells and cytokines to strengthen or attenuate the immunotherapy response, which has provided us targets for combination of immunotherapy [28]. However, definite correlation between autophagy and tumour immune microenvironment has not been study yet, a comprehensive and systematic analysis of them were urgently needed.

In present study, we integrated the transcriptional and genetic profiles of several cohorts to systematically analyse autophagy regulation patterns and established an autophagy phenotype related signature (ATGscore) to represent it in individual. As a result, we found distinct autophagy regulation patterns were highly correlated with TIME infiltration, immune phenotypes, molecular subtypes, as well as clinicopathological characteristics, and ATGscore was a robust independent prognostic and predictive factor for clinical outcome of ICIs immunotherapy.

Data collection and processing

Publicly available transcriptional datasets for HNSCC were systematically searched. Then eleven datasets including E-MTAB-1328, GSE39366, GSE41613, GSE42743, GSE65858, E-TABM-302, GSE6791, GSE30784, GSE40774, GSE84846 and TCGA-HNSCC were enrolled (Additional file 1: Table S1). After filtering out the samples without complete prognosis information, we finally got six datasets (E-TABM-302, GSE6791, GSE30784, GSE40774, GSE84846 and TCGA-HNSCC) with 1165 samples for further study. The microarray datasets were downloaded from Gene-Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (www.ebi.ac.uk/arrayexpress/). The raw signal data in microarray were processed via RMA algorithm background correction, log2 transformation, quantile normalization and annotation by the package “Affy” in R [29]. Relative expression of each gene symbol was annotated as highest when several probes mapped to single gene symbol. Then “ComBat” algorithm of package “sva” in R, which reduce the likelihood of batch effects of non-biological technical biases from each dataset [30], was utilized to merge the five microarray datasets (E-TABM-302, GSE6791, GSE30784, GSE40774, and GSE84846) as a meta-HNSCC cohort. The level 3 fragments Per Kilobase per Million (FPKM) data of TCGA-HNSCC dataset were downloaded from the TCGA Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/). Transcripts per kilobase million (TPM) values were transferred from the FPKM values to represent the relative expression of each gene symbol, which are more similar to gene expression from microarrays and more comparable between samples [31]. Detailed information of clinicopathological characteristics for TCGA-HNSCC dataset could be found in Additional file 2: Table S2. The somatic mutation data processed with MuTect2 algorithm for TCGA-HNSCC cohort was downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov/) using the package “TCGAbiolinks” in R [32]. The total number of mutations counted in the whole exon territory was set as 38 Mb according to a previous study [33]. Moreover, we also enrolled IMvigor210 (mUC) cohort, which included patients with metastatic urothelial cancer receiving PD-L1 inhibitor atezolizumab, to validate the results we found in HNSCC. The raw gene expression and clinical data were retrieved using the package “IMvigor” in R (http://research-pub.gene.com/IMvigor210CoreBiologies) [34]. The detailed clinical information of IMvigor210 (mUC) cohort could be found in Additional file 3: Table S3. Data were analysed with the R (version 3.5.3) and Bioconductor packages.

Unsupervised consensus clustering for autophagy regulation patterns

Autophagy related genes (ATGs) were obtained from Human Autophagy Database (HADb, http://www.autophagy.lu/). Then ATGs were submitted to unsupervised consensus clustering algorithm (K-means) based on Euclidean distance and Ward's linkage for analysis to identify distinct autophagy regulation patterns (ATGclusters) [35]. The number of clusters was determined with a consensus clustering method in TCGA-HNSCC cohort and further validated in meta-HNSCC cohort. The package “ConsensuClusterPlus” in R was applied to performed this procedure and repeated 1000 times to guarantee the stability of classification.

Differentially expressed genes (DEGs) between autophagy regulation patterns

Based on the distinct expression of ATGs, patients with HNSCC were successfully classified into five autophagy regulation patterns through consensus clustering algorithm. DEGs were screened out by comparing the patterns with different function annotation by using the package “edgeR” in R [36], which implements an empirical Bayesian approach to evaluate the change of gene expressions in different groups. The significance criteria for determining DEGs were set as a false discovery rate (FDR) < 0.05 and |log2FC| > 1.0.

Estimation of infiltrating immune cells

Single cell gene set enrichment analysis (ssGSEA), which evaluated the variation in pathway and biological process activity in the single sample of gene expression dataset, was utilized to estimate relative amount of infiltration immune cells in HNSCC TIME by using package “GSVA” in R [37]. The basic unit for ssGSEA was gene set, which is consisted of genes that share common biological function, chromosomal location, or genetic regulation [38]. TIME infiltration immune cell types such as innate immune cells (dendritic cell, eosinophil, mast cell, macrophage, natural killer cell, neutrophil, etc.) and adaptive immune cells (B cell, T cell, T helper cell, CD8+ T cell, regulatory T (Treg) cell and cytotoxic cell, etc.) were obtained and from the study of Bindea and Charoentong [39, 40]. The normalized enrichment score (NES) from the ssGSEA was regarded as relative amount of TIME infiltrating immune cell in HNSCC.

Functional annotation and pathway enrichment analyses

The difference on biological process between distinct autophagy regulation patterns was analysed through ssGSEA algorithm. The gene sets of “ h.all.v7.1.symbols” and “c2.cp.kegg.v7.1.symbols” were downloaded from MSigDB database of Broad Institute (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) to evaluate the regular biological pathway enrichment. In spite of HALLMARK and KEGG signalling pathway, a series of gene sets were curated to represent biological processes related with immune activation, stromal activation, DNA damage repair, and immune checkpoints, which was obtained from Mariathasan et al, including (a) Angiogenesis; (b) Antigen processing machinery; (c) Base excision repair; (d) CD8 T-effector; (e) Cell cycle; (f) Cell cycle regulators (g) DNA damage repair 1 and 2; (h) DNA replication; (i) Epithelial-mesenchymal transition (EMT) 1, 2 and 3; (j) Fanconi anemia; (k) Homologous recombination; (l) Immune-checkpoint; (m) Mismatch repair; (n) Nucleotide excision repair; (o) Pan-fibroblast TGFb response signature (Pan-F-TBRS); (16) WNT targets [34, 41, 42]

Construction of autophagy phenotype related signature

In order to define autophagy phenotype in individual patient, a set of scoring system was needed. The DEGs between autophagy regulation patterns were submitted to univariate cox analysis with cut-off p-value < 0.01 to generate candidate prognostic DEGs. For dimension reduction, we then performed LASSO cox regression algorithm on prognostic DEGs and construct an optimal autophagy phenotype related signature (ATGscore) through package “glmnet” in R [43]. The optimal values of the penalty parameter λ were determined through 10 cross-validations. The ATGscore of each sample was defined by the relative expression of candidate prognostic DEGs within the model and their cox coefficient. The ATGscore = ∑_(i=1)^n (coefi × Expri), where Expri is the relative expression of the DEG in the signature for patient i and coefi is the LASSO cox coefficient of DEG i.

Statistical analyses

The statistical significance of variables between two groups or more that two groups were by Wilcoxon tests or Kruskal-Wallis tests respectively. The differences between survival curves for each group were determined by Kaplan-Meier analysis with the log-rank test using the package “survminer” in R. The distance between different parameters was computed by Pearson and distance correlation analyses. Independent prognostic factors were identified through univariate and multivariate cox proportional hazard analysis and visualized with the package “forestplot” in R. ATGscore was integrated with other independent factors to establish the nomogram and calibration curves through packages “rms”, “nomogramEx” and “regplot” in R. According to Iasonos' suggestion, DCA was used to evaluate the clinical utility of our established nomogram [44]. The mutation landscape of patients with TCGA-HNSCC cohort was visualized with waterfall plot through packages “maftools” [45] and “complexheatmap” [46] in R. Contingency tables such as the ICIs targeting immunotherapy response were analysed by two-sided Fisher’s exact tests. All statistical analyses were performed with R software 3.5.3. Statistical significance was set at p < 0.05.

Characterization of autophagy regulation patterns in HNSCC

A workflow and design of this study were shown in Fig. 1a. TCGA-HNSCC cohort was used as the training set to identify the autophagy regulation patterns through package “ConsensuClusterPlus” in R. We found that ATGs could successful cluster five autophagy regulation patterns with high stability in TCGA-HNSCC cohort, including 175 cases in pattern A, 56 cases in pattern B, 143 cases in pattern C, 191 cases in pattern D, and 33 cases in pattern E, which was also termed as ATGclusterA-E respectively (Fig. 1b). Unsupervised hierarchical clustering demonstrated that the ATGs were significantly differential expressed among ATGclusterA-E in TCGA-HNSCC cohort (Fig. 1c). Moreover, Kaplan-Meier survival curves showed ATGclusterA-E displayed a complete different survival beneficial that patients with ATGclusterB and E have a better prognosis than other clusters (Log-rank test, p = 0.0054; Fig. 1d).

Tumour immune microenvironment (TIME) landscape in distinct autophagy regulation patterns

As autophagy played a dual role in tumour immune microenvironment (TIME), here we aimed to make clear the association between TIME infiltration and autophagy in HNSCC. The landscape of TIME was calculated via ssGSEA algorithm and showed via cluster heat map (Additional file 4: Table S4). The relative amount of TIME immune cells was strikingly different in distinct ATGclusters as follows (Fig. 2a): ATGclusterA exhibited high infiltration with almost all immune cells; While ATGclusterB was remarkably rich in effector immune cells including activated B, activated CD4, CD8 T cells, activated dendritic cells (DCs) and cytotoxic cells, but less infiltrated with immunosuppressive cells such as regulatory T cells (Treg cells), macrophages and mast cells; ATGclusterC/D displayed low infiltration with all immune cells; ATGclusterE was characterized by high infiltration in activated CD8 and cytotoxic cells but low infiltration in regulatory Treg cells, macrophages and mast cells.

Characteristics of biological process in distinct autophagy regulation patterns

The biological process of distinct autophagy regulation patterns was explored. As shown in Additional file 5: Fig. S1a-b and Additional file 6-7: Table S5-6, five distinct ATGclusters were enriched in different HALLMARK and KEGG signalling pathways respectively: ATGclusterA was annotated with pathways associated with immune infiltration and stromal activation such as interferon a/g response, inflammatory response, epithelial mesenchymal transition (EMT) and TGFb signalling pathway, etc.; ATGclusterB was highly featured with immune activation and DNA damage repair (DDR) including antigen processing and presentation, T/B cell receptor signalling pathway, cytokine-cytokine receptor interaction and natural killer cell mediated cytotoxicity (Additional file 5: Fig. S1C); ATGclusterC was enriched in carcinogenic pathways including MTORC1 and TP53 signalling pathways; ATGclusterD/E were both prominently related to DDR related pathways activation.

Distinct autophagy regulation patterns exhibit different immune phenotypes

We found that ATGclusterA was highly infiltrated but patients with this pattern showed a worse prognosis. Recent study have determined three immune phenotypes of tumours: desert, excluded, and inflamed. Immune inflamed phenotype was characterized as high infiltration with immune cells, while immune desert phenotype displayed the opposite situation. Moreover, immune excluded phenotype were featured with infiltration of abundant immune cells, but most immune cells were located in the stroma surrounding the core tumour niche rather than penetrate their parenchyma, which was considered cytotoxic T-cell suppressive [47]. Then we enrolled a specific gene sets from Mariathasan et al to investigate the enrichment of key signalling pathways associated with immune phenotypes. From the TIME landscape and function annotation analysis, ATGclusterA was recognized as immune-excluded phenotype as stromal related signalling pathways including angiogenesis, epithelial-mesenchymal transition (EMT), WNT target, and pan-fibroblast TGFb response signalling pathways (pan-F-TBRS) were strikingly activated, which could hamper the beneficial effect of highly immune cells infiltration (Fig. 2b and Additional file 8: Table S7). Furthermore, ATGclusterB/E was remarkably induced in immune activation signalling pathway including antigen processing machinery, CD8 T effector, and immune checkpoint but deactivated in stromal related signalling pathways, which was the characteristics of immune-inflamed phenotype (Fig. 2b and Additional file 8: Table S7). The results from TIME landscape and pathway enrichment indicated that ATGclusterC/D was more likely to be immune-desert phenotype (Fig. 2b and Additional file 8: Table S7). Then major markers representing for immune phenotypes related signalling pathways were subsequent curated as follows: immune activation: CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF; immune-checkpoints: CD274 (PD-L1), CD80, CD86, CTLA4, IDO1, LAG3, PDCD1 (PD-1), PDCD1LG2 (PD-L2), TIGIT, HAVCR2 (TIM-3) and TNFRSF9; MHC molecules: HLA-A/B/C, HLA-DMA/DMB, HLA-DOA/DOB, HLA-DPA1, HLA-DPB1/DPB2, HLA-DQB1, HLA-DRA, HLA-DRB6, HLA-E/F/G; TGFβ/EMT signalling pathway: ACTA2, COL4A1, PDGFRA, SMAD9, TGFB1/2/3, TGFBR1/2/3, TWIST1/2, VIM, and ZEB1/2. Surprisingly, the expression of above markers displayed a similar distribution as the biological process and pathways enrichment among distinct ATGclusters, which again validated that distinct autophagy regulation patterns exhibited differential immune phenotypes (Fig. 2c and Additional file 5: Fig. S1d-f).

Tumour somatic mutation in distinct autophagy regulation patterns

The relationship between somatic mutation and autophagy was also measured. We found that ATGclusterD exhibited highest tumour mutation burden (TMB), while ATGclusterB was associated with lowest TMB (Fig. 2d). Moreover, top 30 highly variant mutant genes were utilized to plot the somatic mutation landscape among distinct ATGclusters in patients with HNSCC. In Fig. 2e, ATGclusterD displayed highest mutation rate of top 30 mutant genes, especially for TP53, which was identified as key gene in tumorigenesis of HNSCC [48]. But ATGclusterB only showed small amount of TP53 mutation, which was consistent with the TMB calculation. All above results demonstrated that autophagy regulation patterns both correlated with TIME infiltration and tumour mutation landscape, which underlies the indispensable role of autophagy in HNSCC development.

Validation of autophagy regulation patterns in meta-HNSCC cohort

Unsupervised consensus clustering also identified five ATGclusters in meta-HNSCC cohort (Fig. 3a). The cluster heat map showed that the transcriptional profile of ATGs were differential distributed among five ATGclusters (Fig. 3b). Moreover, Kaplan-Meier survival curves demonstrated that the prognosis of five distinct ATGclusters were strikingly different, which patients with ATGclusterB/E lived longer and patients with ATGclusterA/C/D was associated with poorer survival (Log-rank test, p = 0.0022; Fig. 3c).

Furthermore, a similar TIME landscape among five distinct ATGclusters of TCGA-HNSCC was determined in meta-HNSCC cohort as follows (Fig. 3d and Additional file 9: Table S8): ATGclusterA/B/E was more infiltrated with immune cells, while ATGclusterC/D was less infiltrated. Compared with effector immune cells (CD4, CD8 T cells, and cytotoxic cells), ATGclusterA exhibited an abundance of immunosuppressive cells (Treg cells, macrophages and mast cells). The opposite situation was found in ATGclusterE that the amount of effector immune cells was robustly higher against immunosuppressive Treg cells, macrophages and mast cells. A similar phenomenon was seen in ATGclusterB which the ratio of cytotoxic T lymphocyte (CTL) including activated CD8 T cells and cytotoxic cells to Treg cells was highest among five ATGclusters. ATGclusterC/D displayed lowest infiltration of effector immune cells among five ATGclusters.

Then function annotation of ATGclusters in meta-HNSCC cohort showed a similar enrichment of KEGG and HALLMARK signalling pathways with TCGA-HNSCC cohort. ATGclusterA was both activated in immune infiltration and stromal activation related pathways, while ATGclusterB/E was only annotated with immune activation signalling pathways. ATGclusterC/D was both remarkably enhanced in DDR and EMT/ TGFb related signalling pathway (Additional file 10: Fig. S2a-b and Additional file 11-12: Table S9-10). Moreover, ssGSEA of specific gene sets also showed a similar trend. ATGclusterA was both enriched in immune and stromal activation signalling pathways, which could be recognized as immune-excluded phenotype. ATGclusterB/E was more prominently enhanced in immune activation gene sets such as CD8 T effector, and immune checkpoint, which was more likely to immune-inflamed phenotype. And ATGclusterC/D was slightly annotated with DDR related signalling pathways but less enriched in immune infiltration signalling pathways, which was the characteristic of immune-desert phenotype (Fig. 3e and Additional file 13: Table S11). The expression distribution of immune phenotypes related markers displayed a similar trend to function annotation among five ATGclusters (Additional file 10: Fig. S2c-f). All of these again validated that autophagy played an indispensable role in the immune regulation within TIME and might distinguish immune phenotype in HNSCC.

Establishment of autophagy phenotype related signature (ATGscore)

Above findings revealed that there exists distinct autophagy regulation patterns in HNSCC and demonstrated an essential role of autophagy in shaping TIME landscapes, but all analyses were conducted in cohorts with patient population. Due to the heterogeneity of individuals, it was urgent to construct a set of scoring system to quantify the autophagy related pattern in individual patient with HNSCC. By comparing the transcriptomic profiles from main autophagy regulation patterns, we have obtained 5734 phenotype-related meta-DEGs (Additional file 14: Table S12). Then meta-DEGs were submitted to univariate cox regression analyses, 383 of 5734 meta-DEGs which were significantly correlated with prognosis were finally identified as autophagy phenotype candidate genes (Additional file 15: Table S13). Furthermore, LASSO cox regression analysis was utilized for dimension reduction on these genes to construct an autophagy phenotype related signature (ATGscore) which was representative for autophagy related pattern in individuals. At last, 11 genes were selected to establish ATGscore and the formula was as follows: ATGscore = ACTL10*(-0.0729) + C19orf57*(-0.0292) + CHAD*(-0.0801) + FCN2*(-0.3865) + FGB*(0.1961) + GPR174*(-0.0437) + HSF5*(-0.0126) + SERPINA5*(0.1589) + SRPX*(0.0322) + ZNF541*(-0.033) + ZNF831*(-0.3632).

Kaplan-Meier survival curves showed that patients with low ATGscore lived longer than patients with high ATGscore when they were grouped at best cut-off in TCGA-HNSCC cohort (Log-rank test, p < 0.00001; Fig. 4a-b). Then we found that ATGscore was differential distributed among five ATGclusters that ATGclusterA with poor prognosis showed the highest median score while ATGclusterB with good prognosis exhibited the lowest median score (Kruskal-Wallis test, p < 2.2e-16; Fig. 4c). Recently, a comprehensive molecular landscape has been constructed for HNSCC by The Cancer Genome Atlas Network, which classified HNSCC into atypical, basal, classical, and mesenchymal molecular subtypes [49]. Then we surprised to find that low ATGscore group was concentrated on the atypical subtype, which was correlated with low malignancy and good outcome, while other molecular subtypes with high malignancy and bad outcome, significantly accumulated in high ATGscore group (Fig. 4d-e). The correlation matrix showed that ATGscore and autophagy candidate genes was significantly correlated with immune cells which ATGscore was robustly negative associated with almost immune cells, indicating its role in TIME regulation (Fig. 4f and Additional file 16: Table S14). Then ssGSEA of KEGG and HALLMARK gene sets demonstrated that low ATGscore group was enriched in immune infiltration related pathways, but high ATGscore group was induced in EMT/TGF-b signalling pathways (Additional file 17: Fig. S3a-b and Additional file 6-7: Table S5-6). We also detected the relationship between ATGscore and intrinsic gene expression gene sets closely linked to the stromal activation program and immune activation procedure. Function annotation demonstrated that immune activation gene sets such as antigen processing machinery, CD8 T effector, and immune checkpoint were strikingly induced in low ATGscore group, while stromal relevant signalling pathways were enhanced in high ATGscore group (Fig. 4g). ATGscore was found to be negative correlated with immune activation signature and positive correlated with stromal activation signature through Pearson correlation analysis (Fig. 4h and Additional file 18: Table S15). Then the expression of immune phenotype marker genes showed similar trend with results of ssGSEA (Additional file 19: Fig. S4a-d).

Then the role of ATGscore was validated in EMTAB1328, GSE41613, and GSE42743 meta-cohorts, which contained eleven autophagy phenotype candidate genes. Cluster heat map demonstrated that TIME infiltration distribution among five ATGclusters was similar to meta-HNSCC cohort (Fig. 5a and Additional file 20: Table S16). Moreover, we could clearly see that the amount of effector immune cells including CD8 T cells and cytotoxic cells were robustly accumulated in low ATGscore group when compared with high ATGscore group in ATGclusterA/B/E, which was characterized as high immune infiltration (Fig. 5a). Kaplan-Meier survival curves demonstrated that patients with low ATGscore had a better prognosis than patients with high ATGscore (Log-rank test, p = 0.0196; Fig. 5b). Then function annotation of HALLMARK and KEGG gene sets showed that EMT/TGF-b signalling pathways was highly enriched in high ATGscore group, while low ATGscore group was featured with immune infiltration and activation related signalling pathways enrichment (Additional file 21: Fig. S5a-b and Additional file 22-23: Table S17-18). Furthermore, low ATGscore group exhibited upregulation of immune activation, immune checkpoints, and human leukocyte antigen (HLA) relevant markers and downregulation of TGF-b/EMT relevant markers, while high ATGscore displayed the opposite distribution (Additional file 24: Fig. S6a-d). We also found that ATGscore was negative correlated with the expression of immune checkpoints (Fig. 5c). Moreover, ATGscore was significantly positive correlated with angiogenesis, EMT, Pan-F-TBRS, and WNT targets signatures, and negative correlated with antigen processing machinery, CD8 T-effector, and immune-checkpoint signatures (Fig. 5d-e and Additional file 25-26: Table S19-20). All above results strongly suggested that ATGscore is a good representative of autophagy regulation patterns and is competent at distinguishing immune phenotypes in HNSCC.

ATGscore can be utilized as an independent prognostic factor in HNSCC

We next want to make clear the association between ATGscore and clinical characteristics in HNSCC. The results showed that high ATGscore group were more likely to be patients with more advanced pathological TNM stage, HPV negative, and smoker, as well as treatment outcome of stable disease (SD)/ progression disease (PD), with tumour, and dead status, but did not show a correlation with alcohol history. The opposite patterns could be seen in low ATGscore group (Fig. 6a-f, Additional file 27: Fig. S7a-e, and Additional file 2: Table S2). The distribution changes between ATGscore and clinical features were visualized with alluvial diagram (Fig. 6g).

As gene mutation and TMB were remarkably different in five ATGclusters, we next investigated whether they were still different between low and high ATGscore groups in TCGA-HNSCC cohort. The waterfall plots revealed that TP53 and CDKN2A were differential mutated between high and low ATGscore groups (Fig. 6h). Kaplan-Meier survival curves demonstrated that the prognosis of patients with TP53 mutation were robustly worse than patients with TP53 wild type (WT) (Fig. 6i). We were amazed to find that patients with TP53 mutations were more likely to be higher ATGscore (Fig. 6j). Although mutation of CDKN2A was not correlated with patients’ survival, patients with high ATGscore were also more likely to be CDKN2A mutation (Fig. 6k-l). Unfortunately, we found no difference in TMB between the two groups with high and low ATGscore (Additional file 27: Fig. S7f). This might be caused by the mutation data in TCGA-HNSCC cohort as the high TMB patients were correlated with worse prognosis, which is against our common sense (Additional file 27: Fig. S7g).

As highly associated with the above clinical traits in HNSCC, we sought to determine whether ATGscore was responsible for prognosis prediction independent of these clinicopathological characteristics. Patients were divided into different subgroups according to their features. Stratification survival analyses demonstrated that ATGscore could efficiently predict the prognosis of HNSCC patients in almost all the subgroups from the aforementioned clinical features (Additional file 28-29: Fig. S8-9). Furthermore, univariate and multivariate cox regression analyses were enrolled to figure out independent prognostic role of ATGscore in HNSCC. ATGscore, together with other clinical features, including age, gender, histological grade, HPV status, smoking status, histologic grade, and alcohol history, TMB, pathological T stage, pathological N stage, pathological TNM stage, were enrolled as covariates to conduct the analysis. The results demonstrated that the ATGscore, pathological T stage, pathological N stage, and gender were independent factors that could be utilized to predict the prognosis of HNSCC patients (Additional file 30: Fig. S10a-b and Additional file 31: Table S21). Then we constructed a nomogram by integrating above independent prognostic factors to serve as a clinically relevant quantitative method for clinicians to predict mortality in patients with HNSCC (Additional file 30: Fig. S10c). By using the nomogram, each patient would get a total point by adding the points for each prognostic parameter. The patients obtained the higher total points correspond to a worse clinical outcome. The calibration plots and DCA demonstrated that our nomogram had displayed a similar performance to that of an ideal model and had high potential clinical utility (Additional file 30: Fig. S10d-g). All of these demonstrated that ATGscore could not only estimate the autophagy related pattern of individual patients, characterized immune phenotypes, as well as further acted as an independent prognostic factor in HNSCC.

ATGscore was potent to predict clinical response to ICIs immunotherapy

Recently, ICIs immunotherapy has emerged as a major breakthrough in the treatment of solid tumours. Above findings demonstrated that ATGscore could not only distinguish the prognosis of patients, but also determine TIME infiltration and immune phenotype, which indirectly confirmed its potential role in predicting clinical response to ICIs treatment. Then IMvigor210 (mUC) cohort, which was consist of the patients with metastatic urothelial cancer receiving PD-L1 inhibitor with atezolizumab, were introduced to study the role of ATGscore in evaluating immunotherapeutic benefits. Kaplan-Meier survival curves showed that patients with low ATGscore exhibited a significantly better clinical outcome compared with patients with high ATGscore (Log-rank test, p < 0.001, Fig. 7a).

Previous study reported that ICIs treatment responders were more likely to be patients with genomically unstable (GU) subtype in Lund classification system and TCGA II subtype TCGA classification system by using the IMvigor210 (mUC) cohort. Here in cluster heat map, we found that low ATGscore group exhibited an abundance of activated B, activated CD4, CD8 T cells, and cytotoxic cells, but less infiltrated with Treg, macrophages and mast cells in patients with GU and TCGA II subtypes (Additional file 32: Fig. S11A). Moreover, ssGSEA of specific gene sets showed antigen processing machinery (APM), CD8 T effector, and immune checkpoints signature which represented the immune activation were strikingly enriched in low ATGscore group (Fig. 7b). The immune activation related markers showed a similar trend to function annotation (Fig. 8e-f and Additional file 32: Fig. S11b). However, we found EMT signature was slightly activated in low ATGscore group, but Pan-F-TBRS and WNT-target signatures were strikingly enriched in high ATGscore group (Fig. 7b and Additional file 32: Fig. S11c). And most EMT/TGF-b signalling pathways related markers were not differential expressed between two groups, which indicated that ATGscore might not distinguish stromal pattern in IMvigor210 (mUC) cohort. The results of the correlation matrix are almost identical to above findings (Fig. 7c). As IMvigor210 (mUC) cohort contained a complete information of immune phenotypes, we were amazed to find that the patients with immune-inflamed phenotype exhibited lowest ATGscore, while patients with immune-desert phenotype displayed high ATGscore (Kruskal-Wallis test, p = 2.2e-6, Fig. 7d). The numbers of patients with immune-inflamed phenotype were nearly two times in low HPXscore group when compared with that in high HPXscore group (Fisher’s exact tests, p < 0.00001, Fig. 7e). All of these definitely confirm our findings and hypothesis that autophagy was able to shape the TIME infiltration and distinguish the immune phenotype.

Although ATGscore showed no correlation with TMB in TCGA-HNSCC cohort, here we were surprised to find that there existed a negative correlation between TMB and ATGscore in IMvigor210 (mUC) cohort (Wilcoxon test, p = 0.036, Fig. 8a and Fisher’s exact tests, p = 0.033886, Fig. 8b). Kaplan-Meier survival curves revealed that survival benefit of patients with high TMB was superior to patients with low TMB (Log-rank test, p < 0.001, Fig. 8c). Moreover, we combine the information of ATGscore and TMB to find patients with low ATGscore as well as high TMB exhibited a tremendous survival advantage to all other subgroups (Log-rank test, p < 0.001, Fig. 8d). Then we comprehensively explore the role of ATGscore in the whole IMvigor210 (mUC) cohort as it contained so many information associated with immunotherapy response. The alluvial diagram was utilized to show distribution changes of these features according to ATGscore of individual patient in IMvigor210 (mUC) cohort (Fig. 7f). We found that GU subtype and TCGA II subtype, which displayed high somatic mutation and more likely responded to ICIs treatment, demonstrated lowest ATGscore when compared with other molecular subtypes in the Lund and TCGA classification system respectively (Kruskal-Wallis test, p = 8e-8, Additional file 32: Fig. S11d; Kruskal-Wallis test, p = 2.6e-5, Additional file 32: Fig. S11f). And numbers of patients with GU subtype and TCGA II subtype in low ATGscore group were over two times of that in high ATGscore group (Fisher’s exact test, p = 0.010413, Additional file 32: Fig. S11e; Fisher’s exact test, p = 0.116173, Additional file 32: Fig. S11g). Furthermore, the correlation between ATGscore and immune checkpoint PD-L1 located on tumour cells (TC) or immune cells (IC) were also measured. Surprisingly, we found that most of patients with IC2, which was correlated with better clinical outcome of immunotherapy, concentrated on low ATGscore group, while IC0 strikingly accumulated in high ATGscore group (Fisher’s exact test, p = 0.000136, Fig. 8g). Moreover, patients with IC2 exhibited lowest ATGscore and patients with IC0 showed highest ATGscore (Kruskal-Wallis test, p = 5.2e-9, Fig. 8h). Then we found no difference in ATGscore among patients with TC0-TC2 (Kruskal-Wallis test, p = 0.34, Fig. S11h; Fisher’s exact test, p = 0.728113, Fig. S11i). In addition, we were delighted to find that patients with low ATGscore were more likely to be ICIs treatment responders (CR/PR), while patients with high ATGscore tended to non-responders (SD/PD), (Wilcoxon test, p = 0.0046, Additional file 32: Fig. S11j; Kruskal-Wallis test, p = 0.00037, Fig. 8j). Moreover, the numbers of patients with ICIs immunotherapy responders (CR/PR) was over two times in low ATGscore group to that in high ATGscore group (Fisher’s exact test, p = 0.0055, Additional file 32: Fig. S11k; Fisher’s exact test, p = 0.006009, Fig. 8i). All our findings implied that there exist distinct autophagy regulation patterns in tumours, which could shape the TIME infiltration and immune phenotypes, as well as predict the clinical outcome of ICIs immunotherapy.

In recent years, the key topics on cancer research have transferred from focusing on the tumour itself to studying the interaction between the tumour and its surrounding environment, which is commonly referred to tumour microenvironment (TME). Tumour microenvironment, where the tumour cells arise and live, includes not only the tumour cells themselves, but also their surrounding stroma, microvasculature, and a variety of cells including immune cells, stromal cells, fibroblasts, and various other cells, as well as the biological molecules such as cytokines and chemokines they secreted [50]. They cooperated with each other to form chronic inflammatory, immunosuppressive, and tumour promotion environment so that tumour cells can escape from immune surveillance and survive against attack of effector immune cells [51].

Under normal circumstances, autophagy degraded and recycled cytoplasmic components to maintain protein synthesis and other necessary metabolic functions, which is considered to be an endogenous defence mechanism [52–54]. But autophagy could also promote or inhibit tumour progression through many ways depending different backgrounds, which is thought to be a “double-edged sword” in tumours. Autophagy has been reported to regulate the interaction of tumours cells with various substances within surrounding milieu, especially components from the immune system, containing B and T lymphocytes, dendritic cells (DCs), macrophages, and natural killer (NK) cells, as well as the cytokines and immunoglobulins they released [27]. Conversely, immune cells and their secreted cytokines and antibodies could trigger dysfunction of autophagy, which induced or supressed tumorigenesis. However, the autophagy-mediated regulation of the immune system might strengthen or attenuate the effects of immunotherapy. Autophagy was reported to enhance the effect of immunotherapy by exposing antigens to antigen-presenting cells (APCs) and CTL to initiate and execute the tumour recognition and elimination. Otherwise, a growing number of studies have demonstrated that autophagy could also attenuate the effect of immunotherapy by inducing the activation of immunosuppressive atmosphere to hamper the effector immune cells to kill tumour cells [55]. These all indicated that appropriate induction or inhibition of autophagy may propose a prospective therapeutic strategy when combined with chemotherapy, radiotherapy, and immunotherapy. But comprehensive and systematic analysis of correlation between autophagy and tumour microenvironment has not been fully identified yet, which impedes the clinical development of autophagy based activators or inhibitors.

In this study, we gathered ATGs to determine autophagy regulation patterns through unsupervised consensus clustering. Five distinct autophagy regulation patterns (ATGclusters) were identified with differential ATGs expression, survival benefit, tumour immune microenvironment (TIME) infiltration, and functional annotation. In TCGA-HNSCC cohort, ATGclusterA/B exhibited high TIME infiltration, while ATGclusterC/D/E was relatively less infiltrated with immune cells. Kaplan-Meier survival curves showed that the prognosis of patients with ATGclusterB/E were robustly better than patients with ATGclusterA/C/D, which was not consistent with the findings in TIME infiltration. We noticed that ATGclusterA was both infiltrated with effector immune and immunosuppressive cells, while ATGclusterB/E was highly infiltrated with activated CD8 T cells and cytotoxic cells and less infiltrated with Treg cells, macrophages, and mast cells. So we inferred that highly infiltrated immunosuppressive cells could counteract and impair effector immune cells to distinguish and eradicate abnormal tumour cells. As the opposite function of CTL and Treg cell in tumour immunity, a combined assessment of CTL and Treg infiltration have comprehensively studied. The ratio of CTL/Treg has finally been recognized as the independent prognostic factor in many tumour types [56, 57]. Here in our study, we found that CTL/Treg ratio was significantly higher in ATGclusterB/E than that in ATGclusterA/C/D, which well explains the mismatch of immune infiltration and survival analysis. Moreover, function annotation demonstrated ATGclusterA was not only enriched in signalling pathways associated with inflammation, but also induced in angiogenesis, EMT, and pan-fibroblast TGFβ response signalling pathways, which indicated the stromal status was relatively activated in this pattern. The stromal status, which could reversibly shift between “loose” and “dense” status, is the key checkpoint for appropriate localization and migration of T cells into tumour parenchyma. As activated in stromal status, effector immune cells cannot effectively penetrate the stroma surrounding core tumour islets, leaving a large number of CTL trapped in the extracellular matrix, thus failing to execute their antitumor role of immune surveillance and eradication, allowing the tumour to continue to progress in ATGclusterA, which was the characteristic of immune-excluded phenotype [58]. In addition to lack of activated and priming T-cell, ssGSEA of KEGG, HALLMARK and specific signatures all revealed that immune tolerance and ignorance was fully induced in ATGclusterC/D, which was more likely to be to immune-desert phenotype [59]. Moreover, we also found despite highly infiltrated with effector immune cells, the antigen processing machinery, CD8 T effector, and immune checkpoint signature, which were representative for immune activation, was also strikingly enriched in ATGclusterB/E, which was featured as immune-inflamed phenotype. Moreover, the expression of immune activation, stromal activation, MHC molecules, and immune checkpoints was distributed as the same trend to function annotation. Recently, many teams have defined the non-inflamed or immune-suppressed tumours with immune-excluded and immune-desert phenotypes as “cold” tumour and tumours with immune-inflamed phenotype as “hot” tumour, which might be responsible for clinical response of immune-checkpoints inhibitors (ICIs) immunotherapy [60, 61]. Next, we systematically collected transcriptome datasets of HNSCC across GEO and ArrayExpress database and merged them as the meta-HNSCC cohort. We were amazed to find that all above results could be precisely validated in the meta-HNSCC cohort. All of these suggested that distinct autophagy regulation patterns associated with signalling pathway enrichment, TIME infiltration, and immune phenotypes do exist in HNSCC, which provides the possibility of the combination of autophagy activator or inhibitors with ICIs.

Accumulated evidence demonstrated that tumour patients with immune-inflamed phenotype will obtained a durable responses and better overall survival when receiving with ICIs treatment [59]. But we also noticed that not all patients benefit from ICIs-targeting therapy, with an estimated response rate only modestly above the historical 10% response rate to traditional chemotherapies [41]. Immune tolerance to these tumours is still a major impediment in cancer immunotherapy. To improve the efficacy of immunotherapy and prolong survival from immunotherapy, we need to clarify the underlying mechanisms of immune tolerance, which disenabled the role of effector immune cells. Many important factors has been identified to influence the immune tolerance, such as hypo-infiltration of effector immune cells into tumour parenchymal, imbalance between effector immune and immunosuppressive cells, abnormalities in MHC molecules function and expression, and lack of tumour antigen or epitopes exposure leading to failure of antigen processing to T cells [62]. As highly correlated with TIME infiltration and immune phenotypes, the association between autophagy regulation pattern and tumour mutation load was further investigated. Previous study demonstrated that the recognition of neo-antigens exposure, which mainly triggered by somatic nonsynonymous mutations, was essential for initiating the antigen processing and activation of the adaptive tumour immunity cascade. Moreover, tumour mutation burden (TMB), which could be easily assessed to replace the overall neo-antigen detection, has been determined as a potential biomarker for predicting the clinical outcome of ICIs treatment [63, 64]. In TCGA-HNSCC cohort, we found that TMB was differential distributed among five distinct autophagy regulation patterns, where ATGclusterB showed lowest TMB and ATGclusterD showed highest TMB, which is not consistence with the results from TIME infiltration and function annotation. But next we were surprised to find patients with high TMB were associated with worse prognosis than patients with low TMB in TCGA-HNSCC cohort. This indicated that TMB is harmful in the TCGA-HNSCC cohort which was against our common sense. However, we then found that DNA damage repairing signalling pathways were significantly enriched in ATGclusterB/E, which might be another cause of good prognosis in this pattern.

Then autophagy phenotype related genes, which was exacted from the DEGs among distinct autophagy regulation patterns, were submitted to LASSO cox regression analysis to establish a set scoring system to evaluate and quantify autophagy regulation pattern in individuals, which was referred to the autophagy phenotype related signature (ATGscore). ATGscore was found to be differential distributed among five distinct autophagy regulation patterns. Moreover, in TCGA-HNSCC, microarray-HNSCC, and IMvigor210 (mUC) cohorts, we were pleasantly to find that low ATGscore group was significantly enriched in immune activation relevant signalling pathways and deactivated in stromal relevant signalling pathways, which was characteristic of “hot” tumour. Thus the opposite phenomenon was seen in high ATGscore group, which referred to “cold” tumour. In addition, we validated these in IMvigor210 (mUC) cohort to find that immune-inflamed phenotype exhibited lowest ATGscore, while immune-desert and immune-excluded phenotypes displayed higher ATGscore, which indicated the successful of model construction. All above results revealed that ATGscore was not only a reliable tool to assess autophagy regulation pattern, but also was potent to effectively evaluate TIME infiltration and immune phenotype in individual patient. In addition, we found patients with good clinical outcomes, as well as low malignancy clinicopathological traits and molecular subtypes more likely to accumulate in low ATGscore group, while the opposite patterns were observed in high ATGscore group. In IMvigor210 (mUC) cohort, we were amazed to find that patients with GU and TCGA II molecular subtypes, as well as IC2 phenotypes, which was reported that more likely to have an good response to ICIs-targeting immunotherapy, were robustly concentrated in the low ATGscore group, and rarely observed in with high ATGscore group [34, 41]. Moreover, we found that TMB was not different between high and low ATGscore groups. But we next found that high ATGscore group were more likely to be mutation in TP53, which again indicated the role of ATGscore in prediction the ICIs immunotherapy response as mutation status of TP53 could have predictive value in immunotherapy in patients with HNSCC [65]. Although TMB did not work well in TCGA-HNSCC cohort, here in IMvigor210 (mUC) cohort we found that ATGscore was negative correlated with TMB and patients with low ATGscore group were more likely to be patients with high TMB, which was consistent with findings in GU molecular subtype with high mutation load. Moreover, Kaplan-Meier survival curves showed that combination of ATGscore and TMB could significantly improve the predictive value when compared with that of TMB or ATGscore alone where patients with low ATGscore and high TMB exhibited best prognosis. Finally we found patients with low ATGscore were more likely to benefit from ICIs treatment and ICIs targeting immunotherapy responders exhibited a lower ATGscore. According to all above findings, autophagy regulation patterns and ATGscore was found to be significantly correlated with three main factors: pre-existing activated CTL or immunoreactivity, activation of the EMT/TGFβ signalling pathway or stromal status, and tumour neo-antigen or TMB levels, to influence the clinical outcome of ICIs immunotherapy.

In conclusion, we comprehensively and systematically assessed distinct autophagy regulation patterns and established a set scoring system ATGscore which could represent them and was associated with the TIME infiltration, immune phenotypes, molecular subtypes, genetic variation, clinical outcome of ICIs targeting immunotherapy, etc. More importantly, this study has yielded novel insights of combination of autophagy-based inducer or inhibitor with various therapeutic strategies such as immunotherapy for clinical application in HNSCC.

DCA: decision curve analysis; DEGs:differential expression genes; ECM:extracellular matrix; EMT:epithelial-mesenchymal transition; FDR:false discovery rate; GEO:Gene Expression Omnibus; GSEA:Gene Set Enrichment Analysis; GSVA:Gene Set Variation Analysis; HR, hazard ratio; HNSCC:head and neck squamous cell carcinoma; ICIs:immune-checkpoint inhibitors; LASSO:Least absolute shrinkage and selection operator; mUC:metastatic urothelial cancer; ssGSEA:single-sample gene set enrichment analysis; TCGA:The Cancer Genome Atlas; TIME:tumour immune microenvironment; TMB:tumour mutation burden; TME:tumour microenvironment; TNM:Tumour Node Metastasis.

Acknowledgements

The authors are grateful for the invaluable support and useful discussions with other members of the Department of Stomatology. We thank the TCGA, GEO, and ArrayExpress databases as well as IMvigor210 package for the availability of the data.

Authors’ contributions

BM made substantial contributions to conception and design of the research. BM, HL and MZ carried out data collection and analysis. BM, RC and RY wrote the paper. BM, HL, MZ, RC and RY edited the manuscript and provided critical comments. All authors read and approved the final manuscript.

Data accessibility

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. The data that support the findings of this study were derived from the following resource: TCGA GDC data portal (https://portal.gdc.cancer.gov/), GEO) (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress database (www.ebi.ac.uk/arrayexpress/) as well as package “IMvigor” in R (http://research-pub.gene.com/IMvigor210CoreBiologies)

Funding

This work was supported by grants from Beijing Natural Science Foundation (Grant No.7194279), Beijing Shijitan Hospital Natural Science Foundation (Grant No. 2019-C09), National Natural Science Foundation of China (Grant No. 81802550), China Postdoctoral Science Foundation (Grant No. 2019M660041) and Beijing Postdoctoral Research Foundation (Grant No. ZZ2019‐04).

Compliance with ethical standards

Conflict of interest

The authors declare that they have no conflict of interest.

Informed consent

All the analysed data used in our research were collected from a public database, such as TCGA, GEO, and ArrayExpress as well as IMvigor210 package and so on; therefore, informed consent was not required for this analysis.

Ethical approval

Since this was a retrospective study and all the data were collected from a public database TCGA, GEO, and ArrayExpress as well as IMvigor210 package and so on, therefore, ethical approval was not required.

Siegel RL, Miller KD, Jemal A: Cancer Statistics, 2017. CA Cancer J Clin 2017, 67(1):7-30.
Castellsague X, Alemany L, Quer M, Halec G, Quiros B, Tous S, Clavero O, Alos L, Biegner T, Szafarowski T et al: HPV Involvement in Head and Neck Cancers: Comprehensive Assessment of Biomarkers in 3680 Patients. Journal of the National Cancer Institute 2016, 108(6):djv403.
Braakhuis BJ, Snijders PJ, Keune WJ, Meijer CJ, Ruijter-Schippers HJ, Leemans CR, Brakenhoff RH: Genetic patterns in head and neck cancers that contain or lack transcriptionally active human papillomavirus. Journal of the National Cancer Institute 2004, 96(13):998-1006.
Rietbergen MM, Brakenhoff RH, Bloemena E, Witte BI, Snijders PJ, Heideman DA, Boon D, Koljenovic S, Baatenburg-de Jong RJ, Leemans CR: Human papillomavirus detection and comorbidity: critical issues in selection of patients with oropharyngeal cancer for treatment De-escalation trials. Annals of oncology : official journal of the European Society for Medical Oncology 2013, 24(11):2740-2745.
Slebos RJ, Yi Y, Ely K, Carter J, Evjen A, Zhang X, Shyr Y, Murphy BM, Cmelak AJ, Burkey BB et al: Gene expression differences associated with human papillomavirus status in head and neck squamous cell carcinoma. Clin Cancer Res 2006, 12(3 Pt 1):701-709.
Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, Scharrer U, Halama D, Bertolini J, Bauer U et al: The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. International journal of cancer 2015, 137(12):2846-2857.
Gupta S, Kong W, Peng Y, Miao Q, Mackillop WJ: Temporal trends in the incidence and survival of cancers of the upper aerodigestive tract in Ontario and the United States. International journal of cancer 2009, 125(9):2159-2165.
O'Rorke MA, Ellison MV, Murray LJ, Moran M, James J, Anderson LA: Human papillomavirus related head and neck cancer survival: a systematic review and meta-analysis. Oral oncology 2012, 48(12):1191-1201.
Dikic I, Elazar Z: Mechanism and medical implications of mammalian autophagy. Nature reviews Molecular cell biology 2018, 19(6):349-364.
Galluzzi L, Green DR: Autophagy-Independent Functions of the Autophagy Machinery. Cell 2019, 177(7):1682-1699.
Levine B, Klionsky DJ: Development by self-digestion: molecular mechanisms and biological functions of autophagy. Developmental cell 2004, 6(4):463-477.
Doherty J, Baehrecke EH: Life, death and autophagy. Nature cell biology 2018, 20(10):1110-1117.
Kaushik S, Cuervo AM: The coming of age of chaperone-mediated autophagy. Nature reviews Molecular cell biology 2018, 19(6):365-381.
Mizushima N: A brief history of autophagy from cell biology to physiology and disease. Nature cell biology 2018, 20(5):521-527.
Yu L, Chen Y, Tooze SA: Autophagy pathway: Cellular and molecular mechanisms. Autophagy 2018, 14(2):207-215.
Folkerts H, Hilgendorf S, Vellenga E, Bremer E, Wiersma VR: The multifaceted role of autophagy in cancer and the microenvironment. Medicinal research reviews 2019, 39(2):517-560.
Tsuchihara K, Fujii S, Esumi H: Autophagy and cancer: dynamism of the metabolism of tumor cells and tissues. Cancer Lett 2009, 278(2):130-138.
Monkkonen T, Debnath J: Inflammatory signaling cascades and autophagy in cancer. Autophagy 2018, 14(2):190-198.
Cheng Y, Ren X, Hait WN, Yang JM: Therapeutic targeting of autophagy in disease: biology and pharmacology. Pharmacological reviews 2013, 65(4):1162-1197.
White E, DiPaola RS: The double-edged sword of autophagy modulation in cancer. Clin Cancer Res 2009, 15(17):5308-5316.
Liu JL, Chen FF, Lung J, Lo CH, Lee FH, Lu YC, Hung CH: Prognostic significance of p62/SQSTM1 subcellular localization and LC3B in oral squamous cell carcinoma. Br J Cancer 2014, 111(5):944-954.
Tang JY, Hsi E, Huang YC, Hsu NC, Chu PY, Chai CY: High LC3 expression correlates with poor survival in patients with oral squamous cell carcinoma. Hum Pathol 2013, 44(11):2558-2562.
Zhou ZR, Zhu XD, Zhao W, Qu S, Su F, Huang ST, Ma JL, Li XY: Poly(ADP-ribose) polymerase-1 regulates the mechanism of irradiation-induced CNE-2 human nasopharyngeal carcinoma cell autophagy and inhibition of autophagy contributes to the radiation sensitization of CNE-2 cells. Oncology reports 2013, 29(6):2498-2506.
Delgado MA, Elmaoued RA, Davis AS, Kyei G, Deretic V: Toll-like receptors control autophagy. The EMBO journal 2008, 27(7):1110-1121.
Sorbara MT, Ellison LK, Ramjeet M, Travassos LH, Jones NL, Girardin SE, Philpott DJ: The protein ATG16L1 suppresses inflammatory cytokines induced by the intracellular sensors Nod1 and Nod2 in an autophagy-independent manner. Immunity 2013, 39(5):858-873.
Harris J: Autophagy and cytokines. Cytokine 2011, 56(2):140-144.
Jiang GM, Tan Y, Wang H, Peng L, Chen HT, Meng XJ, Li LL, Liu Y, Li WF, Shan H: The relationship between autophagy and the immune system and its applications for tumor immunotherapy. Molecular cancer 2019, 18(1):17.
Wen Y, Graybill WS, Previs RA, Hu W, Ivan C, Mangala LS, Zand B, Nick AM, Jennings NB, Dalton HJ et al: Immunotherapy targeting folate receptor induces cell death associated with autophagy in ovarian cancer. Clin Cancer Res 2015, 21(2):448-459.
Gautier L, Cope L, Bolstad BM, Irizarry RA: affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 2004, 20(3):307-315.
Johnson WE, Li C, Rabinovic A: Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 2007, 8(1):118-127.
Wagner GP, Kin K, Lynch VJ: Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory in biosciences = Theorie in den Biowissenschaften 2012, 131(4):281-285.
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I et al: TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016, 44(8):e71.
Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, Schrock A, Campbell B, Shlien A, Chmielecki J et al: Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome medicine 2017, 9(1):34.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE, III, Koeppen H, Astarita JL, Cubas R et al: TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018, 554(7693):544-548.
Nidheesh N, Abdul Nazeer KA, Ameer PM: An enhanced deterministic K-Means clustering algorithm for cancer subtype prediction from gene expression data. Computers in biology and medicine 2017, 91:213-221.
Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26(1):139-140.
Hanzelmann S, Castelo R, Guinney J: GSVA: gene set variation analysis for microarray and RNA-seq data. BMC bioinformatics 2013, 14:7.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES et al: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(43):15545-15550.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z: Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell reports 2017, 18(1):248-262.
Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A et al: Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013, 39(4):782-795.
Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, Dawson N, O'Donnell PH, Balmanoukian A, Loriot Y et al: Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet 2016, 387(10031):1909-1920.
Senbabaoglu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, de Velasco G, Miao D, Ostrovnaya I, Drill E, Luna A et al: Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol 2016, 17(1):231.
Friedman J, Hastie T, Tibshirani R: Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of statistical software 2010, 33(1):1-22.
Iasonos A, Schrag D, Raj GV, Panageas KS: How to build and interpret a nomogram for cancer prognosis. Journal of clinical oncology : official journal of the American Society of Clinical Oncology 2008, 26(8):1364-1370.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP: Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome research 2018, 28(11):1747-1756.
Gu Z, Eils R, Schlesner M: Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016, 32(18):2847-2849.
Chen DS, Mellman I: Elements of cancer immunity and the cancer-immune set point. Nature 2017, 541(7637):321-330.
Zhou G, Liu Z, Myers JN: TP53 Mutations in Head and Neck Squamous Cell Carcinoma and Their Impact on Disease Progression and Treatment Response. Journal of cellular biochemistry 2016, 117(12):2682-2692.
Cancer Genome Atlas N: Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature 2015, 517(7536):576-582.
Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, Li J, Li F, Tan HB: Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett 2020, 470:126-133.
Wu T, Dai Y: Tumor microenvironment and therapeutic response. Cancer Lett 2017, 387:61-68.
Feng Y, He D, Yao Z, Klionsky DJ: The machinery of macroautophagy. Cell research 2014, 24(1):24-41.
Viry E, Paggetti J, Baginska J, Mgrditchian T, Berchem G, Moussay E, Janji B: Autophagy: an adaptive metabolic response to stress shaping the antitumor immunity. Biochemical pharmacology 2014, 92(1):31-42.
Wang S, Xia P, Rehm M, Fan Z: Autophagy and cell reprogramming. Cellular and molecular life sciences : CMLS 2015, 72(9):1699-1713.
Moskovitz J, Moy J, Ferris RL: Immunotherapy for Head and Neck Squamous Cell Carcinoma. Current oncology reports 2018, 20(2):22.
Alvaro T, Lejeune M, Salvado MT, Bosch R, Garcia JF, Jaen J, Banham AH, Roncador G, Montalban C, Piris MA: Outcome in Hodgkin's lymphoma can be predicted from the presence of accompanying cytotoxic and regulatory T cells. Clin Cancer Res 2005, 11(4):1467-1473.
Sato E, Olson SH, Ahn J, Bundy B, Nishikawa H, Qian F, Jungbluth AA, Frosina D, Gnjatic S, Ambrosone C et al: Intraepithelial CD8+ tumor-infiltrating lymphocytes and a high CD8+/regulatory T cell ratio are associated with favorable prognosis in ovarian cancer. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(51):18538-18543.
Salmon H, Franciszkiewicz K, Damotte D, Dieu-Nosjean MC, Validire P, Trautmann A, Mami-Chouaib F, Donnadieu E: Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. The Journal of clinical investigation 2012, 122(3):899-910.
Kim JM, Chen DS: Immune escape to PD-L1/PD-1 blockade: seven steps to success (or failure). Annals of oncology : official journal of the European Society for Medical Oncology 2016, 27(8):1492-1504.
Gajewski TF: The Next Hurdle in Cancer Immunotherapy: Overcoming the Non-T-Cell-Inflamed Tumor Microenvironment. Seminars in oncology 2015, 42(4):663-671.
Lowe KL, Cole D, Kenefeck R, I OK, Lepore M, Jakobsen BK: Novel TCR-based biologics: mobilising T cells to warm 'cold' tumours. Cancer treatment reviews 2019, 77:35-43.
Gajewski TF, Woo SR, Zha Y, Spaapen R, Zheng Y, Corrales L, Spranger S: Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Current opinion in immunology 2013, 25(2):268-276.
Gibney GT, Weiner LM, Atkins MB: Predictive biomarkers for checkpoint inhibitor-based immunotherapy. Lancet Oncol 2016, 17(12):e542-e551.
Fancello L, Gandini S, Pelicci PG, Mazzarella L: Tumor mutational burden quantification from targeted gene panels: major advancements and challenges. J Immunother Cancer 2019, 7(1):183.
Klinakis A, Rampias T: TP53 mutational landscape of metastatic head and neck cancer reveals patterns of mutation selection. EBioMedicine 2020, 58:102905.

Download PDF

Version 1

posted

You are reading this latest preprint version

Autophagy regulation patterns characterization identifies immune phenotypes and immunotherapy response in head and neck squamous cell carcinoma (HNSCC)

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1