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 (Supplementary 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 Fig. S1A-B and Supplementary 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 (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 . 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 Supplementary 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 Supplementary Table S7). The results from TIME landscape and pathway enrichment indicated that ATGclusterC/D were more likely to be immune-desert phenotype (Fig. 2B and Supplementary 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 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 . 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 Supplementary 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 (Fig. S2A-B and Supplementary 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 Supplementary Table S11). The expression distribution of immune phenotypes related markers displayed a similar trend to function annotation among five ATGclusters (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 (Supplementary 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 (Supplementary 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 . 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 Supplementary 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 (Fig. S3A-B and Supplementary 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 Supplementary Table S15). Then the expression of immune phenotype marker genes showed similar trend with results of ssGSEA (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 Supplementary 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 (Fig. S5A-B and Supplementary 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 (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 Supplementary 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, Fig. S7A-E, and Supplementary 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 (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 (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 (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 (Fig. S10A-B and Supplementary 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 (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 (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 (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 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 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. S7E). 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,Fig. 11D; Kruskal-Wallis test, p = 2.6e-5,Fig. 11F). 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, Fig. 11E; Fisher’s exact test, p = 0.116173, Fig. 11G). 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, 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, 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.