Identification and validation of a tumor mutation burden-related signature combined with immune microenvironment infiltration in adrenocortical carcinoma

1.   Introduction

Adrenocortical carcinoma (ACC) is a rare endocrine malignancy. It afflicts one in every million people each year, and the median overall survival is merely 3–4 years ^[1]. For individuals with local or locally progressive illness, radical resection is presently the sole curative option ^[2]. However, the cumulative recurrence rate is still high even after surgery ^[3]. The most widely used TNM (tumor, lymph node, and metastasis) classification was not satisfactory owing to the lack of predominant genomic and molecular characteristics ^[4,5]. Therefore, identifying pivotal genomic determinants to enhance the predictive accuracy is important for ACC treatment and survival analysis.

In recent years, multiple acknowledged biomarkers for immune responsiveness, including microsatellite instability (MSI), tumor-infiltrating lymphocytes (TILs), and tumor mutation burden (TMB), especially the immune microenvironment infiltration and TMB, have shown great potential in the prediction of advanced or aggressive cancers ^[6,7,8]. TILs constitute the most crucial part of immunity since they can mediate the response of the immune system to chemotherapy, and they have revolutionized the treatments for many malignancies ^[7,9]. In addition, TILs have been confirmed to have a considerable impact on the development of tumors and clinical outcomes in various cancers, including lung cancer, urothelial carcinoma, and colorectal cancer ^[10,11,12]. It has been found that high mast cell infiltration indicates a better survival rate in ACC patients ^[13]. TMB represents the total number of somatic missense mutations in one megabase of genomic regions and has been determined as an emerging biomarker accompanied by immune infiltrates in various malignancies ^[14,15,16]. Notably, previous studies revealed that the TMB level could predict immunotherapy effect and survival outcomes across most cancer types ^[17,18,19]. High TMB in kidney renal clear cell carcinoma patients indicated an awfully poor survival outcome and inhibited immune cell infiltration ^[20]. Yan et al. constructed a prognostic signature by combining TMB and immune cell infiltrates to predict survival outcomes in cutaneous melanoma ^[21].

Increasing evidence revealed that polygenic mutation was related to the carcinogenesis and aggressive progression in ACC, indicating the predictive potential of TMB ^[1,22]. Mutations were transcribed and translated into novel antigens, which could be recognized and targeted by the tumor-immune system ^[23]. More mutations contribute to more antigens, making tumors more immunogenic and responsive to the immune system ^[17]. Nevertheless, merely about 20% of cancer patients could benefit from the immune strategy, which may be due to the involvement of TILs and the status of the tumor immune microenvironment (TIME) ^[24]. What is worse is that no biosignatures for evaluating the status of immune microenvironment infiltration in ACC have been found based on TMB level. Accordingly, it is of critical necessity to investigate the underlying molecular mechanism of immune infiltrates and the role of TMB based on an effective model containing multiple biomarkers for ACC.

In this research, we intended to explore the prognostic role of TMB with the combination of the characteristics of immune infiltrates in an attempt to provide a distinctive perspective for the further development of ACC treatment strategies.

2.   Materials and methods

2.1.   Data acquisition and analysis

First of all, the somatic mutation data of 92 ACC patients were extracted from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). After that, the "Masked Somatic Mutation" data was selected and processed by VarScan software. The Mutation Annotation Format of somatic mutation data was prepared and implemented by the "maftools" R package, which provided a wide range of analysis modules to execute a feature‐rich customizable visualization ^[25]. Then, the gene expression data of 92 ACC samples in HTSeq-FPKM format were obtained. Moreover, we retrieved the corresponding clinical data of all samples. In order to facilitate downstream analysis, all Ensembl gene IDs were converted to gene symbols using an annotation GTF file obtained from GENCODE. Meanwhile, we downloaded transcriptome expression profiles and clinical information of GSE76019, GSE33371, and GSE10927 from the Gene Expression Omnibus (GEO) for validation. The probe matrix of the GSE76019 cohort, including 34 patients, was converted to a gene matrix using the GPL13158 platform. Additionally, the GSE33371 and GSE10927 datasets, containing 23 and 24 samples, respectively, were generated using the GPL570 platform.

2.2.   TMB value calculation and prognostic evaluation

TMB refers to the total number of somatic missense mutations in a megabase of the genomic region, comprising base substitutions, insertions, and deletions. The Perl scripts were developed using the JAVA platform to specifically calculate the mutation frequency of all samples. The average length of the human exons is 38 megabase (Mb). Accordingly, the TMB estimate is equal to the total number of variants/38 for each sample. The calculation of TMB for 92 ACC patients was shown in Table S1. The ACC samples were divided into the low and high TMB groups. Then, the Kaplan-Meier analysis was conducted to compare the survival differences between two groups using the "survival" R package. We further assessed the relationship between TMB levels and clinical variables via Wilcoxon rank-sum test.

2.3.   Differentially expressed gene (DEG) and functional enrichment analyses

The "Limma" package was selected for differential gene expression analysis without normalization between the two TMB groups with the screening criteria of | Fold change (FC) | > 1 ^[26]. The heatmap of all DEGs was analyzed and visualized utilizing the "Pheatmap" R package. Then, we used "org.Hs.eg.db" R package to convert all gene symbols into Entrez IDs for each DEG and implemented the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using "enrichplot", "ggplot2", and "clusterProfiler" packages ^[27,28,29]. Additionally, gene set enrichment analysis (GSEA) was conducted by JAVA software. The "c2.cp.kegg.v7.2 symbols.gmt gene sets" obtained from the MSigDB database were chosen as the reference gene set ^[30]. False Discovery Rate (FDR) < 0.05 was considered as a threshold. Homoplastically, the immune-related genes retrieved from the Immport were conducted by the "VennDiagram" package to select immune-related DEGs between two groups ^[31].

2.4.   Survival analysis

A total of 48 immune-related DEGs were screened out and the prognostic values in the two TMB groups were further evaluated. Ultimately, six core genes associated with survival outcomes were identified by the Kaplan-Meier analysis. Additionally, the Cox regression analysis was applied to verify the prognostic potential of these core genes.

2.5.   Analysis of the immune infiltrates of the core immune-related genes in ACC

The CIBERSORT method was utilized to estimate the immune infractions, which was a newly developed tool that can convert gene expression data of each patient into immune infractions ^[32]. Detailed differential distributions of immune cells were analyzed and shown by the heatmap. Additionally, multiple algorithms containing EPIC, XCELL, MCPCOUNTER, QUANTISEQ, CIBERSORT-ABS, CIBERSORT, and TIMER were commanded to determine the differential abundances of immune infiltrates of core immune-related genes in ACC.

2.6.   Establishment of TMB prognostic signature (TMBPS)

We established a novel TMBPS containing six core immune-related genes and assessed its predictive accuracy for all ACC patients. The formula of TMBPS was as follows: TMBPS = Ʃ (βi × Expi) (I = 6). Moreover, the ROC curve was utilized to evaluate the predictive value of multiple clinical parameters in ACC.

2.7.   Statistical analysis

The Wilcoxon rank-sum test was mainly used for comparisons between two groups based on the non-parametric hypothesis test. Kruskal-Wallis test was applied to analyze two or more categories. All statistical analyses were implemented using the R software (Version 4.1.1), and a P value less than 0.05 was considered to indicate statistically significant differences.

3.   Results

3.1.   The landscape of mutation profiles in ACC

The somatic mutation profiles in ACC were analyzed for a comprehensive landscape of mutation profiles. The percentages (≥ 5%) of the top 34 mutated genes and mutation types marked in different colors were shown in the waterfall plot (Figure 1). On the whole, missense mutation, comprising single nucleotide polymorphisms (SNP) and C > T, was the predominant mutation type (Figure 2A–C). The median number of variants in each sample was 21.5 (Figure 2D), and the variant classifications were represented by different colors (Figure 2E). The simultaneous and exclusive correlation of mutated genes was depicted in Figure 2F. In order to reveal the mutation difference in TMB levels, the landscapes in two TMB groups were compared, as illustrated in Figure 3. Moreover, the gene expression data obtained from TCGA, consisting of 92 ACC patients (32 males and 60 females), along with their clinicopathological features, were summarized in Table 1. The average age of these patients was 47.16 ± 16.30.

3.2.   The correlation of TMB with clinical prognosis

The distributed patterns of clinical features of ACC in two TMB groups were depicted by the heatmap. TMB levels were closely associated with survival status and tumor stage (Figure 4A). Kaplan-Meier survival analysis suggested that ACC patients in the high TMB group tended to have a significantly worse survival outcome (Figure 4B, C). This finding is contrary to the result of previous studies ^[21,33]. Similarly, we assessed the correlation between clinical features and TMB values and found that a higher TMB level was correlated to older age and advanced tumor stage and AJCC-T stage (Figure 4D–F). Nevertheless, no significant correlation was observed between TMB level and gender, AJCC-N stage or AJCC-M stage (Figure 4G–I).

3.3.   Identification of DEGs correlated to TMB

The heatmap showed that the levels of DEGs were generally lower in the low TMB group (Figure 5A). A total of 859 DEGs were determined by differential analysis for the following investigation (Table S2). To elucidate the potential biological functionality and pathways of DEGs, GO and KEGG enrichment analyses were performed. Nuclear division, DNA helicase activity, and microtubule binding were enriched in the GO category (Figure 5C; Table S3). Additionally, KEGG pathway enrichment analysis and the GSEA revealed that cell cycle, DNA replication, p53 signaling pathway, and pathways in cancer were enriched (Figure 5D, E; Tables S4 and S5). Owing to the fact that TMB was associated with the immune microenvironment, 48 immune-related genes were determined for the next analysis (Figure 5B; Table 2).

3.4.   Identification of core genes and their relation to the immune microenvironment

Survival analysis and Cox regression analysis were performed to identify core genes that were significantly correlated with survival outcomes. Six core genes were identified and the detailed workflow of screening was shown in Figure 6. We found ACC patients with higher expression levels of TMSB15A, MMP9, BIRC5, and LTBP1 had worse survival outcomes (Figure 7A–D), while patients with high expression levels of CCL14 and PTGFR had better prognosis (Figure 7E, F). The results of Cox regression analysis also indicated that the six core prognostic genes were significantly related to the survival outcomes in ACC patients (Figure 7G, H). To further assess the underlying relation of these core genes with immune infiltrates in ACC, multiple software was applied. It was found that the expression of the six core genes was robustly associated with the abundance of immune cell subtypes (Figure 8).

3.5.   Relationship between TMB and immune cells in ACC

We evaluated the proportions of 22 immune cells in all ACC samples based on the CIBERSORT method with P < 0.05 (Table S6), which were visualized in the heatmap and box plot (Figure 9A, B). We could find memory resting CD4+ T cells, CD8+ T cells, and M2 macrophage are the predominant immune cell types, while activated dendritic cells and M0 macrophage showed high abundance in the high TMB group. The violin diagram showed that infiltration levels of naive B cells and activated NK cells were relatively higher in the low TMB group, while activated memory CD4+ T cells showed a higher expression level in the high TMB group (Figure 9C).

3.6.   Construction and assessment of TMBPS for ACC patients

In view of the relationship between variants of TMB-related immune genes with poor prognosis and higher immune infiltration, we established a TMBPS comprising the above six core genes and assessed its predictive efficiency for ACC. Based on the Cox regression model, TMBPS was calculated as follows: TMBPS = (0.343290 × TMSB15A + 0.026519 × MMP9 + 0.075941 × BIRC5 + 0.109634 × LTBP1 − 0.798708 × CCL14 − 0.142566 × PTGFR) (Table S7). The relations among survival status, prognostic scores, and core gene expression based on the TMBPS were shown in Figure 10A. ROC curve was introduced and the robust predictive accuracy of TMBPS was illustrated, with AUC = 0.897, higher than that of any other clinical parameters (Figure 10B). In addition, survival analysis revealed that ACC patients with high TMBPS suffered from an unfavorable prognosis (Figure 10C). Similarly, the Cox analysis of the above clinical factors was employed to verify that TMBPS could be an independent and robust prognostic predictor for ACC (Figure 10D, E). In order to validate our results, GSE76019, GSE33371, and GSE10927 derived from the GEO database were applied to assess the effectiveness of ACC prediction by TMBPS (Figure 11).

4.   Discussion

Tumorigenesis is the result of the accumulation of genetic alterations in the DNA interacting with immune infiltrates ^[34]. It has been found that TIME paves a novel way for tumor progression and immune response. TILs have been determined as an emerging indicator that affects the prognosis and therapeutic response in various human cancers ^[35,36,37]. Although patients with locally progressive ACC generally have a high recurrence rate after radical resection, PD-L1 inhibitors could reactivate dormant TILs and a high PD-L1 expression level indicated a longer postoperative survival, which represented a promising strategy for ACC ^[38]. Nevertheless, the therapeutic strategies of ACC are still limited and effective biomarkers for immune responses are lacking.

The effective biomarkers can help to identify patients whose immune system would respond so as to avoid waste of money and severe toxicities for non-responders. TMB has proven to be a novel biomarker to predict immune responses in various malignancies, such as prostate adenocarcinoma, urothelial carcinoma of the bladder, and lung cancer ^[15,39,40]. Luo et al. found that a higher TMB level indicated worse BCR-free survival and TMB was associated with the immune infiltrates in prostate cancer ^[15]. Jiang et al. showed that a combination of TMB, immune infiltrates, and PD-L1 expression is feasible for the prediction of early-stage lung squamous cell carcinoma ^[41]. Nevertheless, there were few discussions about the role of TMB and its underlying connection to immune infiltrates in ACC.

In this research, we comprehensively analyzed and visualized the landscape of mutation profiles of ACC patients. It was found that 75% of ACC patients showed various mutation forms, with missense mutations comprised of SNP and C > T mutations accounting for the most. The two most common mutated genes were TP53 and CTNNB1. It was shown that TP53 was a tumor suppressor protein responding to different cellular stresses through regulating the expression of target genes, inducing cell cycle capture, apoptosis, senescence, and metabolism alterations ^[42,43]. CTNNB1 was a protein regulating cell growth and adhesion between cells ^[44,45]. Analysis of the association between TMB level and survival outcomes indicated that ACC patients with high TMB suffered from a worse prognosis. The results demonstrated that higher TMB levels were closely associated with older age and advanced tumor stage and AJCC-T stage. Accordingly, TMB is an effective predictor that could provide valuable information for immunotherapy in various kinds of cancers, including ACC ^[46,47].

To further elucidate the potential biological functionality and mechanisms of 859 DEGs between two TMB groups, we conducted differential analysis to single out TMB-related DEGs. GO enrichment analysis showed that DEGs were primarily involved in cell mitosis. KEGG pathway enrichment analysis and GSEA revealed that DEGs were mainly correlated with cancer progression and immune cell response to tumors, such as DNA replication, cell mitosis, and cell cycle. DNA replication and cell mitosis were the fundamental biological processes in which dysregulation could cause genome instability ^[48]. The accumulative errors of DNA replication and cell mitosis could cause tumorigenesis, including ACC ^[49]. These functions were also correlated with the occurrence of cancer, the influence of the tumor microenvironment, and the differentiation and activation of immune cells. Given that the physiological process of cell mitosis demands the homeostasis of the cell cycle, its dysregulation would bring about the disorder of cell growth and the occurrence of cancer ^[50].

Survival analysis was performed and six core prognostic TMB-related immune genes in ACC were identified. CCL14 was a chemokine inducing the activation of immune cells and a potential prognostic biomarker and tumor suppressor via regulating the cell cycle and promoting apoptosis in hepatocellular carcinoma ^[51]. BIRC5, as an immune-related gene, was greatly related to multiple immune cell infiltrates in diverse cancers and could inhibit apoptosis and facilitate cell proliferation ^[52]. We found BIRC5 was greatly related to abundant immune cell subtypes. On the whole, expression levels of these core immune genes were correlated to the abundance of immune infiltrates, including neutrophil cells, macrophage, cancer-associated fibroblast, B cells, and T cells.

We also analyzed the correlation between TMB level and immune infiltrates to reflect the status of the TIME in ACC. In this study, infiltration levels of naive B cells and activated NK cells in the low TMB group were higher, while activated memory CD4+ T cells showed a higher infiltration level in the high TMB group. The possible reason is that the increased amount of neoantigens caused by genomic mutation promoted the immune activation and recognition of memory CD4+ T cells. It has been shown that CD4+ TILs were extraordinarily associated with antigen processing, and the infiltration of memory T cells was involved in the prognosis of multiple malignancies ^[53,54]. Activated memory CD4+ T cells, which were stimulated by the proliferation of inactive ones, were able to release inflammatory cytokines, thereby promoting tumor growth and accelerating tumor metastasis ^[55]. These findings verified that immune cells played an extremely important role in antitumor immunity in ACC patients. The identification of effective immunological biomarkers can help to avoid immunotherapy resistance and improve the therapeutic effect, and would become a novel promising therapeutic strategy for ACC patients.

Finally, a new risk score signature containing six core genes was established and its predictive value for ACC was assessed with AUC. It was found that ACC patients with high TMBPS had a poor prognosis than those with low TMBPS. Our results were validated in the other three independent ACC patient cohorts retrieved from the GEO database. Therefore, we confirmed the superiority and effectiveness of our risk signature for the diagnosis and treatment of ACC, and it was expected to be applied in clinical practice in the future.

5.   Conclusions

In summary, we revealed a systematic landscape of TMB and identified a TMBPS combined with immune microenvironment infiltration in ACC. This paper will provide a reference for the development of ACC treatment strategies.

Acknowledgments

We are grateful to Guangdong Medical Science and Technology Research Fund Project (B2022332, B2020127) granted to Qingbiao Chen and Foshan Science and Technology Innovation Project (2020001005794) granted to Jingbo Lin.

Conflict of interest

The authors declare no competing interest.