Construction of a risk signature for adrenocortical carcinoma using immune-related genes.

BACKGROUND: Adrenocortical carcinoma (ACC) is considered a rare tumor with a dismal prognosis. Expression of immune-related genes (IRGs) in ACC and correlations between IRGs and ACC prognosis were assessed using The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases.
METHODS: To preliminarily predict immune cell infiltration, an immune score was calculated using ESTIMATE. Differentially expressed IRGs were screened, and potential biological functions were investigated. We then performed univariate Cox regression to identify IRGs associated with survival, and the regulatory mechanisms of IRGs associated with survival were predicted. We built a risk signature through multivariate Cox regression to evaluate patient overall survival (OS).
RESULTS: A high immune score predicted a good prognosis and an early clinical stage in ACC. We identified 30 IRGs associated with survival and then predicted associated regulatory mechanisms via protein-protein interaction (PPI) and transcription factor (TF) regulatory networks. The risk signature established by multivariate Cox regression correlated significantly with prognosis in ACC.
CONCLUSIONS: The vital roles of IRGs in ACC were assessed, and the risk signature obtained based on IRGs associated with survival independently predicted ACC prognosis. 2020 Translational Andrology and Urology. All rights reserved.

Introduction

Adrenocortical carcinoma (ACC), which originates from the cortex of the adrenal gland, is considered a rare tumor with a low incidence rate of 2 per million people (1). However, among endocrine carcinomas, the aggressiveness of ACC is only lower than that of anaplastic thyroid cancer, and patients with ACC at an advanced clinical stage have a dismal prognosis. The median overall survival (OS) time and 5-year survival rate of ACC patients are 3.21 years and 15% to 44%, respectively (2-5). At present, the only effective treatment for ACC is radical resection, which is recommended only when the size of the tumor exceeds 5 cm or increased circulating adrenal hormone levels are confirmed (6). Adjuvant therapies, such as mitoxantrone combined with radiotherapy, administered after surgery can delay ACC recurrence (7).

Recent studies have indicated that a disrupted immune system, the main factor that allows tumor cells to evade the immune response, is involved in the development and progression of tumors (8). As a promising antitumor strategy, immunotherapy aims to induce the immune system to recognize cancer antigens as foreign antigens and suppress the proliferation and metastasis of tumor cells through active and passive immune responses, thus accelerating the development of personalized medicine (9,10). Many studies have suggested a relationship between ACC and immunity. For example, one study showed that among the three ACC subsets integrated through cluster analysis, immune-mediated pathways were significantly upregulated in the subset with the best prognosis (11). Loss of major histocompatibility complex, class II, DR beta 1 (HLA-DRB1), a major histocompatibility complex (MHC) class II allele, and altered expression of Fas/Fas ligand are possible mechanisms of immune escape in ACC (6). A novel adjuvant for immunotherapy, autologous dendritic cells can induce tumor-specific responses in ACC patients, even though final outcomes have not been satisfactory (12). High levels of CD8+ T lymphocytes can independently improve prognosis in childhood ACC (13). Furthermore, upregulation of the immune markers neutrophil-lymphocyte ratio (NLR) and platelet-lymphocyte ratio (PLR) correlates significantly with a poor prognosis in ACC patients after resection (14). Although these findings show that immunity is important in ACC, the associated molecular mechanisms, particularly those underlying immunogenomic effects, remain unclear. Hence, we assessed expression of immune-related genes (IRGs) in ACC and investigated correlations between IRGs and ACC prognosis using public datasets. Gene expression and clinical data for ACC patients were downloaded from The Cancer Genome Atlas (TCGA), and gene expression data for the normal adrenal gland were obtained from the Genotype-Tissue Expression (GTEx) database. The authors present the REMARK reporting checklist (available at http://dx.doi.org/10.21037/tau-20-485).

Methods

Characteristics of the included ACC patients

Our study was a retrospective analysis. Transcriptomic data and corresponding clinical data for ACC patients were downloaded from GTEx and TCGA in the UCSC Xena database (http://xena.ucsc.edu/). We normalized and merged the transcriptome data from 79 ACC tissues and 127 normal adrenal tissues from the two databases into a single dataset to identify differentially expressed genes (DEGs). Seventy-seven ACC patients with both gene expression data and corresponding clinical data were included for additional analysis. OS was used as a prognostic indicator, and the median follow-up time was 1,171 days. Data for the 77 ACC patients are shown in . IRG information was selected from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org/) (15). To further explore regulatory mechanisms, we also obtained transcription factor (TF) information from the Cistrome Cancer database (http://cistrome.org/) (16). All procedures performed in this study were in accordance with the Declaration of Helsinki (as revised in 2013) and no ethical approval was required because the data we used were obtained from public databases. Because of the retrospective nature of the research, the requirement for informed consent was waived.

Table 1

Characteristics of the included ACC patients obtained from the TCGA database

Basic information	Total (n=77)	%
Age	48 (median)
Gender
Female	48	62.3
Male	29	37.7
Stage
I	9	11.6
II	37	48.1
III	16	20.8
IV	15	19.5
T classification
T1	9	11.7
T2	42	54.5
T3	8	10.4
T4	18	23.4
N classification
N0	68	88.3
N1	9	11.7
M classification
M0	62	80.5
M1	15	19.5

ACC, adrenocortical carcinoma.

Immune score

An immune score was obtained using the transcriptome data of ACC tissues in TCGA via Estimation of Stromal and Immune cells in Malignant Tumor tissues with Expression data (ESTIMATE), an algorithm that calculates immune and stromal scores based on expression of immune cells and stromal cell-related genes in tumors to predict infiltration of these cells. The current study only focused on the immune score (17,18). Then, we generated Kaplan-Meier survival curves to illustrate the relationship between ACC patient OS and the immune score (the median value of the immune score was used as the cutoff value).

DEGs

We compared expression of genes in 79 ACC tissue samples with that in 127 normal tissue samples via the edgeR package in R software (19). |Log2 fold change (FC)| >1.0 and a false discovery rate (FDR) adjusted to P<0.05 were defined as the cutoff criteria. Differentially expressed IRGs and TFs were identified from the DEGs. We carried out Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to estimate the biological functions of the differentially expressed IRGs initially via the clusterProfiler package in R software (20) and visualized these results using the GOplot package in R software (21).

Molecular mechanisms of IRGs associated with survival

We performed univariate Cox regression to extract IRGs associated with survival, setting a P value <0.01 as the cutoff. Protein-protein interaction (PPI) network analysis with the String database (score >0.4) (https://stringdb.org/) was conducted to assess relationships among the IRGs associated with survival in ACC patients. In the PPI network, each node represents the protein product of the IRGs associated with survival and each edge the interaction between the proteins (22). Genes with the largest number of edges, also called hub genes, were distinguished in our PPI network. Moreover, as regulatory mechanisms were the focus of the current study, we constructed a TF regulatory network (P value <0.05 and correlation coefficient >0.4 defined as the cutoff criteria) using Cytoscape software version 3.7.2 for visualization.

Risk signature

IRGs associated with survival were subjected to multivariate Cox regression to determine the coefficient of each selected IRG, and the risk score was calculated based on the formula below. (expression level of IRGs associated with survival) (23).

We divided patients into two groups (the median value of the risk score was used as the cutoff value): high-risk and low-risk groups. A receiver operating characteristic (ROC) curve was built to evaluate the accuracy of the risk score in predicting the prognosis of ACC patients (24). Fisher’s exact tests were applied to explore differences in clinical variables between the groups. Kaplan-Meier survival curves and Cox regression were utilized to confirm that the risk score was able to independently predict OS. The clinical variables included for evaluating the effectiveness of the signature were age, gender, TNM classification and clinical stage.

Statistical analysis

The Wilcoxon rank-sum test was used for comparisons between two groups, whereas the Kruskal-Wallis test with Bonferroni’s post hoc test was performed for multigroup comparisons. The survival curve represented by the Kaplan-Meier survival curve was analyzed for the difference in OS between the two groups via the log-rank test. Cox regression was conducted by using the survival package in R software.

Results

We initially found that compared with a low score, a high immune score predicted a better prognosis and early clinical-stage disease in ACC, indicating that immunotherapy may be effective for ACC treatment (). Relatively high immune scores were also observed in male patients and patients with advanced TNM-stage disease, but there were no significant differences (). These results showed that a high level of immune cell infiltration can affect the prognosis of ACC patients and is related to some early clinical variables.

Figure 1

Relationships between the immune score and overall survival (A), clinical stage (B), gender (C), T classification (D), N classification (E) and M classification (F). Genes differentially expressed between ACC tissue samples and normal adrenal gland samples are illustrated using a heatmap (G) and volcano plot (H). We extracted differentially expressed IRGs using a heatmap (I) and volcano plot (J). We screened differentially expressed transcription factors (TFs) via a heatmap (K) and volcano plot (L). GO analysis of differentially expressed IRGs was performed and visualized via a bar plot (M) and dot plot (N). KEGG analysis of differentially expressed IRGs was performed and visualized via a bar plot (O) and dot plot (P). IRGs associated with survival were determined with a P value <0.01 set as the cutoff value in univariate Cox regression (Q). A PPI network (score >0.4) was established to assess relationships among IRGs associated with survival in ACC patients (R). The PPI network showed that IL-6 contained the largest number of edges; thus, it was identified as the hub gene of IRGs associated with survival in ACC. We visualized the number of gene edges in the PPI network using a bar plot. (S). A TF regulatory network was constructed based on 33 differentially expressed TFs and 30 IRGs associated with survival (T). In H, J and L, the red and green dots represent upregulated and downregulated genes, respectively. In X, the red and green circles represent high-risk and low-risk IRGs, respectively. Purple triangles represent TFs. Red indicates upregulation, and green indicates downregulation. The red and green edges represent positive and negative regulation, respectively. IRGs, immune-related genes; ACC, adrenocortical carcinoma; GO, Gene Ontology; BP, biological process; CC, cell component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; TFs, transcription factors.

DEGs in ACC

We identified 5,045 DEGs from all genes in the transcriptomic dataset obtained from UCSC Xena; 2,646 of the DEGs were upregulated and 2,399 downregulated (). Of the 5,045 DEGs, 227 differentially expressed IRGs were screened, with 93 exhibiting upregulated expression and 134 exhibiting downregulated expression (). Finally, we identified 33 differentially expressed TFs, which consisted of 22 with upregulated expression and 11 with downregulated expression (). The results of GO and KEGG analyses for the differentially expressed IRGs are illustrated in . According to the GO results, receptor ligand activity was the most frequent GO biological process category. The “cytokine-cytokine receptor interaction” was identified as the most enriched pathway in KEGG analyses.

IRGs associated with survival

We extracted 30 IRGs significantly associated with survival via univariate Cox regression (all P values <0.01) (). A PPI network that included all 30 IRGs associated with survival was established, with scores of >0.4 in String (). In the PPI network, interleukin 6 (IL-6) contained the largest number of edges; thus, it was identified as the hub gene of the IRGs associated with survival in ACC (). To further investigate associated regulatory mechanisms, a TF regulatory network was constructed based on the 33 differentially expressed TFs and 30 IRGs associated with survival (). In addition, we evaluated correlations between clinicopathological factors and IRGs associated with survival ().

Figure 2

Increased expression of AMH was significantly related to an advanced clinical stage (Aa) and a high T classification (Ab). Increased expression of BIRC5 was significantly related to an advanced clinical stage (Ba), a high T classification (Bb) and a high M classification (Bc). Decreased expression of CCL14 was significantly related to an advanced clinical stage (Ca), a high T classification (Cb), a high N classification (Cc) and a high M classification (Cd). Overexpression of CCL22 was significantly related to female (D). Decreased expression of CXCL1 (E), FPR2 (F), or HLA-G (G) was significantly related to high N classification. Downregulation of KCNH2 expression was significantly related to a high T classification (H). Downregulation of QRFP expression was significantly related to an advanced clinical stage (Ia) and high T (Ib), N (Ic), and M (Id) classifications. Increased expression of SEMG1 was significantly related to an advanced clinical stage (Ja) and high T (Jb), N (Jc), and M (Jd) classifications. Decreased expression of SLPI was significantly related to an advanced clinical stage (Ka) and high T (Kb) and M (Kc) classifications. Increased expression of TMSB15A was significantly related to an advanced clinical stage (La) and high T (Lb) and M (Lc) classifications. Downregulation of VIPR1 expression was significantly related to an advanced clinical stage (Ma) and high T (Mb), N (Mc), and M (Md) classifications. IRGs, immune-related genes; AMH, anti-Mullerian hormone; BIRC5, baculoviral IAP repeat containing 5; CCL14, C-C motif chemokine ligand 14; CCL22, C-C motif chemokine ligand 22; CXCL1, C-X-C motif chemokine ligand 1; FPR2, formyl peptide receptor 2; HLA-G, major histocompatibility complex, class I, G; KCNH2, potassium voltage-gated channel subfamily H member 2; QRFP, pyroglutamylated RFamide peptide; SEMG1, semenogelin 1; SLPI, secretory leukocyte peptidase inhibitor; TMSB15A, thymosin beta 15a; VIPR1, vasoactive intestinal peptide receptor 1.

High risk scores indicated a poor prognosis in ACC

The risk score was calculated from coefficients determined by multivariate Cox regression. A risk signature was then constructed, dividing the patients into two groups (the median value of the risk score was used as the cutoff value), namely, a high-risk group and a low-risk group (). The coefficients of the included genes obtained from multivariate Cox regression are presented in .

Figure 3

Expression of included IRGs in different groups (A). The distribution of ACC patients into different groups (B). Survival status of the patients in the high-risk group or the low-risk group (C). The AUC was 0.970 (D). The results of Fisher’s exact tests and Kaplan-Meier survival analysis revealed that a high risk score was significantly related to a poor OS (E). Univariate Cox regression demonstrated that clinical stage, T classification, M classification and the risk score could predict OS (F). Multivariate Cox regression showed that T classification and the risk score were independent prognostic factors (G). ACC, adrenocortical carcinoma; IRGs, immune-related genes; AUC, area under curve; ROC, receiver operating characteristic; OS, overall survival.

Table 2

The coefficients of included genes obtained from multivariate Cox regression

Gene	Coefficient	P
HLA−G	−2.85909	0.0002
ULBP1	−5.99312	0.0002
SLPI	2.41088	0.0001
CXCL1	6.04584	0.0010
TMSB15A	1.14829	0.0723
MMP9	3.45316	0.0000
SEMG1	−3.33718	0.0004
WNT5A	2.76633	0.0011
KCNH2	0.88324	0.0473
PTGS2	−1.12506	0.1308
BIRC5	−2.96852	0.0014
CCL14	−1.60270	0.0235
CCL22	1.96742	0.0004
SEMA3D	2.83077	0.1287
SEMA7A	3.04236	0.0001
FPR2	−22.91990	0.0017
AMH	3.13594	0.0002
NPPA	12.59732	0.0012
NUDT6	−6.26884	0.0003
QRFP	6.26674	0.0032
SCT	0.80979	0.0014
UTS2	−2.64917	0.0044
NR1H3	−2.36130	0.0831
NR6A1	2.95381	0.0004
VIPR1	−6.43599	0.0982

The area under the ROC curve (AUC) was 0.970, suggesting that the risk signature had the potential to accurately predict the prognosis of ACC patients (). The results of Kaplan-Meier survival analyses and Fisher’s exact tests revealed that a high risk score correlated significantly with an advanced clinical stage (P=0.011), a high T classification (P=0.004), a high M classification (P=0.010) and poor OS (P<0.001) ( and ). Univariate Cox regression demonstrated that T classification [hazard ratio (HR) =10.286; 95% confidence interval (CI): 3.976–26.608; P<0.001), M classification (HR =6.150; 95% CI: 2.710–13.959; P<0.001), clinical stage (HR =6.476; 95% CI: 2.706–15.498; P<0.001) and the risk score (HR =2.707; 95% CI: 2.070–3.539; P<0.001) could predict prognosis (, ). We then carried out multivariate Cox regression, which showed that the risk score (HR =2.819; 95% CI: 1.973–4.030; P<0.001) was associated with prognosis (). Therefore, we considered the risk score to be an independent prognostic factor for ACC.

Table 3

Differences in the characteristics of ACC patients between the high risk and low risk

Basic information	Low risk	High risk	P value
Total	39	38
Age			0.257
≤48	22	16
>48	17	22
Gender			1.000
Female	24	24
Male	15	14
Stage			0.011
I&II	29	17
III&IV	10	21
T classification			0.004
T1&T2	32	19
T3&T4	7	19
N classification			0.310
N0	36	32
N1	3	6
M classification			0.010
M0	36	26
M1	3	12

ACC, adrenocortical carcinoma.

Discussion

Currently, immunotherapy, which exploits the immune system to fight tumors, is a potential treatment for various types of malignancies. Nevertheless, new findings suggest that tumor cells can escape the immune response by using immune checkpoints, such as programmed death-1 (PD-1), programmed death ligand-1 (PD-L1) and cytotoxic T lymphocyte antigen-4 (CTLA-4), in the tumor microenvironment (25). Therefore, research on the application of immune checkpoint inhibitors to prevent immune escape is receiving much attention at present (26-28). In ACC, PD-L1 expression is detected on the membrane of tumor cells, and silencing of mismatch repair genes, which can improve sensitivity to anti-PD-1 therapy, has been demonstrated in 30% of ACC patients, suggesting that immunotherapy may be an effective treatment for ACC (29-31). A recent study showed that IFNγ activates immune cell infiltration in anti-ACC therapy but that it might increase expression of PD-L1, suppressing the immune response. Although it ultimately improved the postoperative disease-free survival (DFS) of ACC patients, the effect of the immune response was not satisfactory. Therefore, the application of anti-PD-L1 therapy might alleviate its suppression of the immune response (32). In addition to expression of immune checkpoint molecules, reduced signaling and expression of the innate immune receptor Toll-like receptor 4 (TLR4) was recently reported in ACC, providing a novel marker for immunotherapy in ACC (33). Although immunotherapy provides hope for the treatment of ACC, there are still some challenges. For instance, immune cytolytic activity in ACC is lower than that in other cancers, and immune responses are inhibited by the secretion of steroid hormones, an important feature of ACC; both of these are considered major barriers to effective immunotherapy in ACC (25,34). Unfortunately, research on immunology and ACC is limited, and the scarcity of case samples caused by the low incidence of ACC may be the main reason. Therefore, further analyses of IRGs using transcriptomic data from public databases will be important for the prediction of ACC prognosis and the development of immunotherapeutic strategies.

In this study, we initially identified that compared with a low score, a high immune score predicted a better prognosis and early clinical stage disease in ACC, indicating that immunotherapy may be effective for ACC treatment. Differentially expressed IRGs and TFs were then extracted, and KEGG functional enrichment analysis showed that the five most significant signaling pathways for activation of the differentially expressed IRGs were the PI3K-Akt signaling pathway, JAK-STAT signaling pathway, chemokine signaling pathway, Ras signaling pathway and MAPK signaling pathway. Previous studies have shown that ACC cell proliferation is partially activated by PI3K-Akt and MAPK signaling (35-37). Mutations in Ras-related genes have also been observed in ACC, and we found that the JAK-STAT signaling promoted ACC cell aldosterone secretion (38-42). Studies on the specific mechanisms underlying chemokine signaling and the relationships between these five signaling pathways and immunotherapy in ACC have rarely been reported. Thus, we propose for the first time that these pathways are associated with immunotherapy efficacy in ACC.

IRGs associated with survival were screened by univariate Cox regression. IL-6 was confirmed to be the hub gene in the PPI network. Among the IRGs associated with survival in ACC patients, a higher level of IL-6 was detected in ACC samples than in normal samples; a high IL-6 level can stimulate the secretion of glucocorticoids, and the serum level of matrix metallopeptidase 9 (MMP-9) decreases significantly after radical resection (43-45). Then, we constructed a TF regulatory network based on IRGs associated with survival and differentially expressed TFs. Among all differentially expressed TFs in ACC patients, centromere protein A (CENPA) (46), E2F transcription factor 1 (E2F1) (47) and forkhead box M1 (FOMX1) (48) showed upregulated expression that was involved in ACC progression and predicted a poor prognosis. In contrast, downregulation of transcription factor 21 (TCF21) expression resulted in the accumulation of secreted glucocorticoids and accelerated proliferation of ACC cells, indicating that TCF21 is a potential prognostic marker in ACC (49-51). Additionally, a previous study reported that Activin could induce x-zone apoptosis to inhibit ACC cell growth induced by SMAD family member 2 (Smad2) (52). However, no studies on other IRGs associated with survival or TFs in ACC are available, and previous publications only offer limited information about the regulatory mechanisms involved. Hence, we constructed a TF regulatory network for further analysis. Conclusively, a risk signature was built from the coefficients of selected IRGs determined by multivariate Cox regression, and it was validated that the risk score can serve as an independent prognostic factor for ACC.

Nonetheless, the deficiencies in our study should be acknowledged. First, we merged transcriptomic data from two databases due to the small number of ACC samples in TCGA and the lack of paired normal tissue samples. Second, our analysis of IRG functions was not verified by in vitro or in vivo experiments. Third, the infiltration of each specific immune cell type was not estimated clearly.

In conclusion, we assessed the roles of IRGs in ACC and investigated correlations between IRGs and ACC prognosis. The signature based on IRGs may be used as a tool to predict prognosis in patients with ACC, and our results provide preliminary evidence for the application of immunotherapy in ACC. Although our findings offer a novel perspective for immunotherapy in ACC, further studies are needed.

The article’s supplementary files as