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

cox regression to the autophagy regulation integrated

which has provided us targets for combination of immunotherapy [28]. However, de nite 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 pro les 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 in ltration, 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 le 1: Table S1). After ltering out the samples without complete prognosis information, we nally 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 ve 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 le 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 le 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 classi cation.
Differentially expressed genes (DEGs) between autophagy regulation patterns Based on the distinct expression of ATGs, patients with HNSCC were successfully classi ed into ve 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 signi cance criteria for determining DEGs were set as a false discovery rate (FDR) < 0.05 and |log2FC| > 1.0.

Estimation of in ltrating 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 in ltration 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 in ltration 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 in ltrating 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. In order to de ne 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 pvalue < 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 de ned by the relative expression of candidate prognostic DEGs within the model and their cox coe cient. The ATGscore = ∑_(i=1)^n (coe × Expri), where Expri is the relative expression of the DEG in the signature for patient i and coe is the LASSO cox coe cient of DEG i.

Statistical analyses
The statistical signi cance 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 identi ed 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 signi cance was set at p < 0.05.

Characterization of autophagy regulation patterns in HNSCC
A work ow 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 ve 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 signi cantly differential expressed among ATGclusterA-E in TCGA-HNSCC cohort (Fig. 1c). Moreover, Kaplan-Meier survival curves showed ATGclusterA-E displayed a complete different survival bene cial 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 in ltration and autophagy in HNSCC. The landscape of TIME was calculated via ssGSEA algorithm and showed via cluster heat map (Additional le 4: Table S4). The relative amount of TIME immune cells was strikingly different in distinct ATGclusters as follows (Fig. 2a): ATGclusterA exhibited high in ltration 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 in ltrated with immunosuppressive cells such as regulatory T cells (Treg cells), macrophages and mast cells; ATGclusterC/D displayed low in ltration with all immune cells; ATGclusterE was characterized by high in ltration in activated CD8 and cytotoxic cells but low in ltration 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 le 5: Fig. S1a-b and Additional le 6-7: Table S5-6, ve distinct ATGclusters were enriched in different HALLMARK and KEGG signalling pathways respectively: ATGclusterA was annotated with pathways associated with immune in ltration and stromal activation such as interferon a/g response, in ammatory 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 le 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 in ltrated but patients with this pattern showed a worse prognosis. Recent study have determined three immune phenotypes of tumours: desert, excluded, and in amed. Immune in amed phenotype was characterized as high in ltration with immune cells, while immune desert phenotype displayed the opposite situation. Moreover, immune excluded phenotype were featured with in ltration 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 speci c 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-broblast TGFb response signalling pathways (pan-F-TBRS) were strikingly activated, which could hamper the bene cial effect of highly immune cells in ltration ( Fig. 2b and Additional le 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-in amed phenotype ( Fig. 2b and Additional le 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 le 8: Table S7). Then major markers representing for immune phenotypes related signalling pathways were subsequent curated as follows: . 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 le 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 identi ed 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 in ltration 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 identi ed ve ATGclusters in meta-HNSCC cohort (Fig. 3a). The cluster heat map showed that the transcriptional pro le of ATGs were differential distributed among ve ATGclusters (Fig. 3b). Moreover, Kaplan-Meier survival curves demonstrated that the prognosis of ve 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 ve distinct ATGclusters of TCGA-HNSCC was determined in meta-HNSCC cohort as follows ( Fig. 3d and Additional le 9: Table S8): ATGclusterA/B/E was more in ltrated with immune cells, while ATGclusterC/D was less in ltrated. 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 ve ATGclusters. ATGclusterC/D displayed lowest in ltration of effector immune cells among ve 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 in ltration 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 le 10: Fig. S2a-b and Additional le 11-12: Table S9-10). Moreover, ssGSEA of speci c 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-in amed phenotype. And ATGclusterC/D was slightly annotated with DDR related signalling pathways but less enriched in immune in ltration signalling pathways, which was the characteristic of immune-desert phenotype ( Fig. 3e and Additional le 13: Table S11). The expression distribution of immune phenotypes related markers displayed a similar trend to function annotation among ve ATGclusters (Additional le 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.
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 ve 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 classi ed HNSCC into atypical, basal, classical, and mesenchymal molecular subtypes [49]. Then we surprised to nd 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, signi cantly accumulated in high ATGscore group (Fig. 4d-e). The correlation matrix showed that ATGscore and autophagy candidate genes was signi cantly correlated with immune cells which ATGscore was robustly negative associated with almost immune cells, indicating its role in TIME regulation ( Fig. 4f and Additional le 16: Table S14). Then ssGSEA of KEGG and HALLMARK gene sets demonstrated that low ATGscore group was enriched in immune in ltration related pathways, but high ATGscore group was induced in EMT/TGF-b signalling pathways (Additional le 17: Fig. S3a-b and Additional le 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 le 18: Table S15). Then the expression of immune phenotype marker genes showed similar trend with results of ssGSEA (Additional le 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 in ltration distribution among ve ATGclusters was similar to meta-HNSCC cohort ( Fig. 5a and Additional le 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 in ltration (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 in ltration and activation related signalling pathways enrichment (Additional le 21: Fig. S5a-b and Additional le 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 le 24: Fig. S6a-d). We also found that ATGscore was negative correlated with the expression of immune checkpoints (Fig. 5c). Moreover, ATGscore was signi cantly 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 le 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 le 27: Fig. S7ae, and Additional le 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 ve 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 nd 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 le 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 le 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. Strati cation survival analyses demonstrated that ATGscore could e ciently predict the prognosis of HNSCC patients in almost all the subgroups from the aforementioned clinical features (Additional le 28-29: Fig. S8-9). Furthermore, univariate and multivariate cox regression analyses were enrolled to gure 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 le 30: Fig. S10a-b and Additional le 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 le 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 le 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 ndings demonstrated that ATGscore could not only distinguish the prognosis of patients, but also determine TIME in ltration and immune phenotype, which indirectly con rmed 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 bene ts. Kaplan-Meier survival curves showed that patients with low ATGscore exhibited a signi cantly 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 classi cation system and TCGA II subtype TCGA classi cation 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 in ltrated with Treg, macrophages and mast cells in patients with GU and TCGA II subtypes (Additional le 32: Fig.  S11A). Moreover, ssGSEA of speci c 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 le 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 le 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 ndings (Fig. 7c). As IMvigor210 (mUC) cohort contained a complete information of immune phenotypes, we were amazed to nd that the patients with immunein amed 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-in amed 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 de nitely con rm our ndings and hypothesis that autophagy was able to shape the TIME in ltration and distinguish the immune phenotype.
Although ATGscore showed no correlation with TMB in TCGA-HNSCC cohort, here we were surprised to nd 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 bene t 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 nd 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 classi cation system respectively (Kruskal-Wallis test, p = 8e-8, Additional le 32: Fig. S11d; Kruskal-Wallis test, p = 2.6e-5, Additional le 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 le 32: Fig. S11e; Fisher's exact test, p = 0.116173, Additional le 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 nd 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 le 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 le 32: Fig.  S11k; Fisher's exact test, p = 0.006009, Fig. 8i). All our ndings implied that there exist distinct autophagy regulation patterns in tumours, which could shape the TIME in ltration and immune phenotypes, as well as predict the clinical outcome of ICIs immunotherapy.

Discussion
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, broblasts, 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 in ammatory, 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][53][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, 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 identi ed 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 identi ed with differential ATGs expression, survival bene t, tumour immune microenvironment (TIME) in ltration, and functional annotation. In TCGA-HNSCC cohort, ATGclusterA/B exhibited high TIME in ltration, while ATGclusterC/D/E was relatively less in ltrated 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 ndings in TIME in ltration. We noticed that

ATGclusterA was both in ltrated with effector immune and immunosuppressive cells, while ATGclusterB/E was highly in ltrated with activated CD8 T cells and cytotoxic cells and less in ltrated
with Treg cells, macrophages, and mast cells. So we inferred that highly in ltrated 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 in ltration have comprehensively studied. The ratio of CTL/Treg has nally been recognized as the independent prognostic factor in many tumour types [56,57]. Here in our study, we found that CTL/Treg ratio was signi cantly higher in ATGclusterB/E than that in ATGclusterA/C/D, which well explains the mismatch of immune in ltration and survival analysis. Moreover, function annotation demonstrated ATGclusterA was not only enriched in signalling pathways associated with in ammation, but also induced in angiogenesis, EMT, and pan-broblast 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 speci c 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 in ltrated 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-in amed 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 de ned the non-in amed or immunesuppressed tumours with immune-excluded and immune-desert phenotypes as "cold" tumour and tumours with immune-in amed 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 nd 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 in ltration, 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-in amed phenotype will obtained a durable responses and better overall survival when receiving with ICIs treatment [59]. But we also noticed that not all patients bene t 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 e cacy 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 identi ed to in uence the immune tolerance, such as hypo-in ltration 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 in ltration 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 ve distinct autophagy regulation patterns, where ATGclusterB showed lowest TMB and ATGclusterD showed highest TMB, which is not consistence with the results from TIME in ltration and function annotation. But next we were surprised to nd 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 signi cantly 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 ve distinct autophagy regulation patterns. Moreover, in TCGA-HNSCC, microarray-HNSCC, and IMvigor210 (mUC) cohorts, we were pleasantly to nd that low ATGscore group was signi cantly 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 nd that immune-in amed 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 in ltration 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 nd 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 ndings in GU molecular subtype with high mutation load. Moreover, Kaplan-Meier survival curves showed that combination of ATGscore and TMB could signi cantly 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 bene t from ICIs treatment and ICIs targeting immunotherapy responders exhibited a lower ATGscore. According to all above ndings, autophagy regulation patterns and ATGscore was found to be signi cantly 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 in uence the clinical outcome of ICIs immunotherapy.

Conclusion
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 in ltration, 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.

Declarations
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 ndings 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  The lines in the boxes represented median value and the black dots showed outliers. The statistical difference among ve distinct autophagy regulation patterns was tested by the Kruskal-Wallis test. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). ATGscore were more likely to be dead and patients with low ATGscore were more inclined to alive. Red indicated patients were dead and Blue indicated patients were alive. Dot-plot indicated that survival time of patients with high ATGscore were less than patients with low ATGSscore. Red indicated dead patients   was tested by Wilcoxon test. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001).   (mUC) cohort (Log-rank test, p < 0.001). (E-F) Immune activation relevant markers (e) and immunecheckpoints relevant markers (f) were differential expressed between high and low ATGscore groups in