Identification of an immune-related gene signature as a prognostic target and the immune microenvironment for adrenocortical carcinoma
Chengdang Xu1 | Caipeng Qin2 | Jingang Jian3 | Yun Peng2 | Xinan Wang1 | Xi Chen1 | Denglong Wu1 | Yuxuan Song2,4
1Department of Urology, Tongji Hospital, School of Medicine, Tongji University, Shanghai, China
2Department of Urology, Peking University People’s Hospital, Beijing, China
3Department of Urology, The First Affiliated Hospital of Soochow University, Dushu Lake Hospital Affiliated to Soochow University, Suzhou Medical College of Soochow University, Suzhou, China
4Department of Urology, Tianjin Medical University General Hospital, Tianjin, China
Correspondence
Denglong Wu, Tongji Hospital, 389 Xincun Rd, Putuo District, 200065 Shanghai, China. Email: wudenglong2009@tongji.edu.cn
Yuxuan Song, Peking University People’s Hospital, 11 Xizhimen South Dajie, 100044 Beijing, China. Email: yuxuan_song2013@163.com
Funding information
Science and Technology Commission of Shanghai Municipality, Grant/Award Numbers: 21ZR1458300, 19ZR1448700, 22ZR1456800, 20Y11904400; Shanghai Hospital Development Center,
Grant/Award Numbers: SHDC2020CR3074B, SHDC12019112
Abstract
Background: Adrenocortical carcinoma (ACC) is a rare endocrine malignancy. Even with complete tumor resection and adjuvant therapies, the prognosis of patients with ACC remains unsatisfactory. In the microtumor environment, the impact of a disordered immune system and abnormal immune responses is enormous. To improve treatment, novel prognostic predictors and treatment targets for ACC need to be identified. Hence, credible prognostic biomarkers of immune- associated genes (IRGs) should be explored and developed.
Material and methods: We downloaded RNA-sequencing data and clinical data from The Cancer Genome Atlas (TCGA) data set, Genotype-Tissue Expression data set, and Gene Expression Omnibus data set. Gene set enrichment analysis (GSEA) was applied to reveal the potential functions of differentially expressed genes.
Results: GSEA indicated an association between ACC and immune-related functions. We obtained 332 IRGs and constructed a prognostic signature on the strength of 3 IRGs (INHBA, HELLS, and HDAC4) in the training cohort. The high- risk group had significantly poorer overall survival than the low-risk group (p <. 001). Multivariate Cox regression was performed with the signature as an independent prognostic indicator for ACC. The testing cohort and the entire TCGA ACC cohort were utilized to validate these findings. Moreover, external validation was conducted in the GSE10927 and GSE19750 cohorts. The tumor-infiltrating immune cells analysis indicated that the quantity of T cells, natural killer cells, macrophage cells, myeloid dendritic cells, and mast cells in the immune microenvironment differed between the low-risk and high-risk groups.
Conclusion: Our three-IRG prognostic signature and the three IRGs can be used as prognostic indicators and potential immunotherapeutic targets for ACC.
@ 2022 The Authors. Immunity, Inflammation and Disease published by John Wiley & Sons Ltd.
Inhibitors of the three novel IRGs might activate immune cells and play a synergistic role in combination therapy with immunotherapy for ACC in the future.
KEYWORDS
adrenocortical carcinoma, immunology, prognosis, signatures, tumor microenvironment
1 INTRODUCTION |
Adrenocortical carcinoma (ACC) is a rare endocrine malignancy worldwide.1,2 The incidence is reported to be one to two per million per year. The prognosis is poor with a median overall survival (OS) of 3-4 years.3 Even with complete tumor resection and adjuvant therapies, the prognosis of patients with ACC following treatment remains unsatisfactory.4 The risk of recurrence and metastasis ranged from at least 15% to 40% among patients with ACC after initial treatment.5,6 Conse- quently, novel prognostic predictors and treatment targets for ACC are needed and credible prognostic biomarkers of immune-associated genes (IRGs) should be explored.
In the microtumor environment, the impact of a disordered immune system and abnormal immune responses is enormous.7 Immune cells may modify the signaling molecule to influence the immunological response, which significantly impacts the risk of recur- rence, change, and diffusion of the neoplasm.8,9 Previous research studies have reported that some tumor cells may evade immune surveillance and immune elimination, resulting in tumor infiltration and metastasis.1º An abnormal immune state is thought to be associated with glioblastoma progression.11 Martinez-Bosch et al.1º sug- gested that immunosuppressive macrophages and myeloid-derived suppressive cells may significantly suppress the immune microenvironment in pancreatic cancer. Peng et al.12 identified two immune-related genes could serve as potential biomarkers for immunotherapy in ACC, which provide better insights into ACC microenvironment. In addition, N6-methyladenosine methylation regulators had impacts on ACC prognosis through regulating immune-related functions.13 Bio- informatics analyses based on The Cancer Genome Atlas (TCGA) constructed a competitive endogenous RNA network associated with tumor-infiltrating immune cells and determined that the immune system is implicated in the microenvironment of ACC.14
Aberrantly expressed IRGs may be potential therapeutic targets and prognostic biomarkers for patients with ACC. For this purpose, we explored IRGs and constructed an IRG prognostic signature. This study was designed using the gene
expression profiles from the TCGA ACC and Genotype- Tissue Expression (GTEx) data sets. To investigate the IRG signature’s feasibility, we validated the predictive power in two Gene Expression Omnibus (GEO) data sets. We also evaluated the correlation between the IRG signature and tumor-infiltrating immune cells.
2 MATERIAL AND METHODS |
2.1 TCGA ACC and GTEx data sets |
The TCGA ACC data set provided ACC specimens (n =79), whereas the GTEx data set provided normal adrenal specimens (n =258). We acquired the RNA- sequencing (RNA-seq) data and clinical samples from the GTEx (http://www.gtexportal.org/home/index.html) and TCGA (http://portal.gdc.cancer.gov/) databases. Differ- entially expressed genes (DEGs) were aberrated tran- scripts identified under different biological context through RNA-seq analysis.15 DEGs makes RNA-seq analysis possible to get a deep understanding of pathogenesis-related molecular mechanisms and biologi- cal functions. Therefore, differential analysis has been regarded as valuable for diagnostics, prognostics, and therapeutics of tumors. To identify DEGs from the ACC and normal adrenal samples, we used the edgeR R package.16 The cutoff criteria were |log Fold Change| ≥ 1 and an adjusted p (adj. p) <. 05.
2.2 .2 | Gene set enrichment analysis (GSEA) and identification of immune-related genes
GSEA (http://www.broadinstitute.org/gsea/index.jsp) was applied to reveal the potential functions of DEGs. Gene set permutations were performed 1000 times for each analysis. p <. 05 was established as the cutoff point for Gene Ontology (GO) enrichments in the GSEA. We search for IRGs from the GSEA v4.2.3 (http://www.gsea-msigdb.org/gsea/downloads.jsp: Immune system process M13664, Immune response M19817)17 and identified 332 IRGs.
2.3 Identification of an immune-related gene prognostic signature |
We construct the IRG signature by using the survival R package. The identification and verification criteria for IRG signatures in the ACC specimens were as follows: (1) intact information regarding the IRG expression and clinical features (age, gender, tumor TNM stage, and survival time); and (2) specimens with total survival <30 days were excluded (some nonneoplastic factors such as severe infection or hemorrhage may have destroyed these specimens). From these, we selected 77 ACC samples for the IRG signature construction. The 77 specimens were randomly divided into the training (n =39) and the testing (n = 38) cohorts. The training cohort was utilized for the IRG signature construction. The other cohort was utilized for validation of the IRG signature.
Next, we screened for prognosis-related IRGs (p <. 01), regardless of whether the IRGs were protective or deleterious according to hazard ratio (HR) as well as 95% confidence interval (CI). The 332 IRGs were analyzed by univariate Cox regression analysis in the training cohort. After that, we selected an optimal model from the prognosis-related IRGs using the Akaike information criteria (AIC) method. The three IRGs with minimum AIC values were finally selected for the prediction model construction.
We analyzed the prognosis-related IRGs using multivariate Cox regression analyses to construct a prognostic IRG signature and identify the coeffi- cients.18-20 The following in-house formula was used to calculate the risk score of the prognostic IGR signal for each specimen: Risk score = Expression IRG1 X CoefficientIRG1 + Expression IRG2 × CoefficientIRG2 + … + ExpressionIRGn X CoefficientIRGn. The risk score of the prognostic IRG signature may be determined based on a linear combination of the IRG expression level weighted by the regression coefficients. We acquired the coefficient by log-transformed HR, derived from the multivariate Cox regression analy- sis.21,22 We divided the specimens into the low-risk and high-risk groups based on the median value of the risk score. Univariate and multivariate Cox regression analyses, time-dependent receiver operating charac- teristic (ROC) curve, and Kaplan-Meier (KM) survival curve were conducted to evaluate the predictive value of the risk score in the training cohort.
2.4 4 | | Internal validation and external validation of the IRG signature
To further validate the predictive ability of the IRG signature, we used the testing cohort and the entire
Genotype-Tissue Expression dataset 258 normal adrenal samples
The Cancer Genome Atlas Adrenocortical Carcinoma dataset 79 adrenocortical carcinoma samples
Gene Expression Omnibus database
5609 Differentially Expressed Genes (DEGs)
332 Immune-related Genes (IRGs)
GSE19027 cohort (n = 24) and GSE19750 cohort (n = 22)
3519 Up-regulated DEGs
2090 Down-regulated DEGs
Identification of a 3-IRG signature in training cohort
Gene set enrichment analysis (GSEA)
The 3-IRG signature associated with survival in training cohort
External Validation in GSE19027 cohort and GSE19750 cohort
Immune-related functions
Validation in testing and entire TCGA ACC cohorts
Cox regression
Stratification analyses
Functional enrichment analysis
Survival analysis and Correlation analysis of the 3 IRGs (INHBA, HELLS and HDAC4)
Tumor-Infiltrating Immune Cells Analysis
TCGA ACC cohort to perform internal validation, and used GSE10927 and GSE19750 cohorts from GEO to perform external validation. The validation was con- ducted by following the same analysis methods of the training cohort.
2.5 Functional enrichment analysis between different risk groups and correlation analysis |
Principal component analysis (PCA) was carried out to profile the expression patterns of the low-risk and
high-risk groups according to the IRG signature by using scatterplot3d R package. We identified DEGs in both groups and adopted the clusterProfiler and GOplot R packages for GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses23 to uncover the potential capacities in the two risk groups and the TCGA ACC cohort. Furthermore, the adj. p <. 05 was used as the cutoff value. We investigated the pairwise gene correlation of three IRGs (Inhibin subunit BA [INHBA], Helicase, lymphoid specific [HELLS], and Histone deacetylase 4 [HDAC4]) by adopting the Pearson correlation analysis. The R was determined. p <. 05 was established as the cutoff value.
(A)
(B)
humoral immune response (GO:0006959)
0.4 -
Volcano
Running Enrichment Score
NES=1.461
0.3
P=0.027
10
0.2
5
0.1
logFC
0
0.0
?
-10
Ranked list metric
9
6
0
20
40
60
80
100
3
-log10(FDR)
0
(C)
(D)
1000
Rank in Ordered Dataset
2000
3000
immune effector process (GO:0002252)
immune response (GO:0006955)
0,0
0.0
Running Enrichment Score
NES =- 1.364
Running Enrichment Score
NES =- 1.270
-0.1
P-value=0.026
-0.1
P-value=0.038
-0.2
-0.2
-0.3-
-0.3
0
0
Ranked list metric
Ranked list metric
-5
-5
-10
-10
-15
-15
500
1000
2000
500
1000
Rank in Ordered Dataset
1500
Rank in Ordered Dataset
1500
2000
2.6 Tumor-infiltrating immune cells of the three-IRG signature |
We retrieved 22 types of tumor-infiltrating immune cells from CIBERSORT, a web tool for analyzing tumor-infil trating immune cells (http://cibersort.stanford.edu/).24,25 To identify the differences between the low-risk and high-risk groups, we measured the quantity of each immune cell. p <. 05 was set as the cutoff value.
2.7 Statistical analysis |
Statistical analysis was carried out by GraphPad Prism 7.0 and R software (v3.5.3: http://www.r-project.org). The RNA-seq data were log2-transformed. To assess the influences of the IRG prognostic signature on OS and other clinical features, Cox regression analyses, log-rank test, and KM method were used. The x test and Student’s t test were adopted to measure qualitative variables and quantitative variables, respectively. According to the IRG signature and other indexes, we evaluated the predictive value using time-dependent ROC. The cutoff value of the two-sided p was set at .05 (two-sided p <. 05).
3 RESULTS |
Figure 1 shows the study flowchart.
3.1 Identification of DEGs and functional enrichment analysis |
All genes from the TCGA ACC and GTEx data sets were analyzed via gene differential expression analysis. From these, we selected 79 ACC specimens and 258 normal adrenal specimens. We used edgeR R package to identify DEGs between ACC and normal adrenal samples to further explore the immune-related DEGs and immune- related functions. In total, we identified 5609 DEGs from these specimens, including 3519 upregulated and 2090 downregulated DEGs (Figure 2A). The biological capaci- ties of DEGs determined via GO enrichment analysis was identified by GSEA. The results indicated that upregu- lated genes were involved in the humoral immune response (GO:0006959; p = . 027; Figure 2B). Moreover, downregulated genes correlated with the effector process of the immune system (GO:0002252; p =. 026) and immune response (GO:0006955; p = . 038; Figure 2C,D), which showed that the tumorigenesis and progress of ACC are related to the immune system and immune responses.
3.2 Identification of immune-related genes and the clinical characteristics of TCGA patients with ACC |
We identified and obtained 332 IRGs from the Molecular Signatures Database v4.0. To construct and validate the IRG signature, 77 patients with ACC were included. These patients were randomly divided into the training (n = 39) and the testing (n = 38) cohorts. Table 1 was the clinical features of 77 ACC patients.
3.3 Identification of the IRG prognostic signature in the training cohort |
We identified prognostic IRGs from the 262 IRGs in the training cohort using univariate Cox regression analysis. There were seven prognosis-related IRGs. To construct an IRG signature, a multivariate Cox regression model with the smallest AIC was applied to the seven prognostic IRGs. Ultimately, three IRGs were selected and we constructed the three-IRG prognostic signature. Table 2 summarizes the details relating to the univariate Cox regression analysis, descriptions, coefficients, Ensembl IDs,
| Characteristics | Entire TCGA ACC cohort (N = 77) | Detailed data | a pª | |
|---|---|---|---|---|
| Training cohort (N= 39) | Testing cohort | |||
| (N=38) | ||||
| Age at diagnosis (years) | .737 | |||
| ≤60 | 60 (77.9%) | 31 (79.5%) | 29 (76.3%) | |
| >60 | 17 (22.1%) | 8 (20.5%) | 9 (23.7%) | |
| Gender | .206 | |||
| Male | 29 (37.7%) | 12 (30.8%) | 17 (44.7%) | |
| Female | 48 (62.3%) | 27 (69.2%) | 21 (55.3%) | |
| Tumor stage | .809 | |||
| Stage I | 9 (11.7%) | 6 (15.4%) | 3 (7.8%) | |
| Stage II | 37 (48.1%) | 18 (46.2%) | 19 (50.0%) | |
| Stage III | 16 (20.8%) | 8 (20.5%) | 8 (21.1%) | |
| Stage IV | 15 (19.4%) | 7 (17.9%) | 8 (21.1%) | |
Abbreviations: ACC, adrenocortical carcinoma; TCGA, The Cancer Genome Atlas.
ax2 test.
| Gene symbol | Ensembl ID | Description | Coefficient | Univariate Cox regression analysis | ||
|---|---|---|---|---|---|---|
| HR | 95% CI | p | ||||
| INHBA | ENSG00000122641 | Inhibin subunit BA | 1.647 | 5.193 | 1.111-24.262 | 0.006 |
| HELLS | ENSG00000119969 | Helicase, lymphoid specific | 2.583 | 9.882 | 1.956-49.925 | <0.001 |
| HDAC4 | ENSG00000068024 | Histone deacetylase 4 | 2.610 | 13.242 | 2.919-60.073 | 0.002 |
Abbreviations: CI, confidence interval; HR, hazard ratio.
(A)
(D)
Risk
High risk
Low risk
Survival time (years)
~
.
Dead
1.00
.
Alive
20
AD
+
N
Survival probability
0.75
0
0
10
20
30
40
(B)
0.50
Patients (increasing risk socre)
+
High risk
0.25
·
p=3.666e-05
Risk score
~
low Risk
0
0.00
N
0
1
2
3
4
5
8
7
8
9
10
11
12
Time(years)
0
10
20
30
40
(C)
Patients (increasing risk socre)
Risk
High risk
19
17
9
7
5
3
2
2
2
2
2
1
1
type
type
Low risk
20
20
17
12
9
7
5
4
2
2
2
1
1
1.4
high
INHBA
1.2
0
1
2
3
4
5
6
7
8
9
10
11
12
low
(E)
Time(years)
1
HELLS
0.8
0.6
HDAC4
0,4
0.2
1.0
(F)
pvalue
Hazard ratio
I
0.8
Age
0.751
1.008(0.962-1.056)
True positive rate
1
0.6
Gender
0.801
1.178(0.330-4.202)
0.4
Stage
0.752
1.118(0.561-2.229)
0.2
risk score (AUC=0.923)
Age (AUC=0.534)
Gender (AUC=0.504)
riskScore
<0.001
2.585(1.569-4.258)
0.0
Stage (AUC=0.811)
0
1
2
3
4
0.0
0.2
0.4
0.6
0.8
1.0
Hazard ratio
False positive rate
and gene symbols. All the three IRGs (INHBA, HELLS, and HDAC4) were deleterious IRGs. For all three IRGs, the univariate Cox HR was greater than 1, which suggested that patients with high expression of the three IRGs may have a shorter survival time. The following formula was established to calculate the risk score of the three- IRG signature based on the expression of the three IRGs for OS prediction: Risk score = (1.647 x Expression INHBA) + (2.583 × ExpressionHELLS) + (2.61 x Expression HDAC4).
3.4 Assessing the predictive capacity of the three-IRG prognostic signature for patients with ACC in the training cohort |
From the training cohort, 39 patients with ACC patients were divided into two groups based on the median value of the three-IRG signature risk score, namely, the high-risk (n = 19) and low-risk (n = 20) groups. Figure 3A,B summarizes the details of each patient in the training cohort, including survival time, survival status, and the risk score. Figure 3C
| Variables | Univariate analysis | Multivariate analysis | ||||
|---|---|---|---|---|---|---|
| HR | 95% CI | p | HR | 95% CI | p | |
| Training cohort | ||||||
| Three-IRG signature risk score | 1.003 | 0.968-1.040 | .866 | 1.008 | 0.962-1.056 | .751 |
| Age at diagnosis | 1.505 | 0.496-4.563 | .470 | 1.178 | 0.330-4.202 | .801 |
| Gender | 2.579 | 1.441-4.617 | .001 | 1.118 | 0.561-2.229 | .752 |
| Tumor stage | 2.718 | 1.752-4.216 | <. 001 | 2.585 | 1.569-4.258 | <. 001 |
| Testing cohort | ||||||
| Three-IRG signature risk score | 1.029 | 0.984-1.076 | .209 | 1.067 | 1.002-1.137 | .044 |
| Age at diagnosis | 0.907 | 0.288-2.862 | .868 | 0.847 | 0.1750-4.103 | .837 |
| Gender | 3.499 | 1.700-7.202 | .001 | 3.434 | 1.264-9.329 | .016 |
| Tumor stage | 1.932 | 1.372-2.720 | <. 001 | 2.168 | 1.281-3.668 | .004 |
| Entire TCGA ACC cohort | ||||||
| Three-IRG signature risk score | 1.010 | 0.985-1.036 | .430 | 1.019 | 0.987-1.052 | .236 |
| Age at diagnosis | 1.043 | 0.473-2.299 | .917 | 0.863 | 0.353-2.112 | .747 |
| Gender | 2.903 | 1.844-4.569 | <. 001 | 1.707 | 1.049-2.780 | .031 |
| Tumor stage | 1.920 | 1.562-2.360 | <. 001 | 1.842 | 1.471-2.307 | <. 001 |
| GSE10927 cohort | ||||||
| Three-IRG signature risk score | 1.012 | 0.977-1.048 | .519 | 0.992 | 0.9589-1.026 | .632 |
| Age at diagnosis | 1.471 | 0.510-4.242 | .475 | 1.251 | 0.326-4.802 | .745 |
| Gender | 1.818 | 1.095-3.019 | .021 | 2.330 | 1.298-4.184 | .005 |
| Tumor stage | 1.863 | 1.143-3.036 | .013 | 2.660 | 1.345-5.261 | .005 |
| GSE19750 cohort | ||||||
| Three-IRG signature risk score | 1.035 | 0.994-1.079 | .097 | 1.019 | 0.978-1.063 | .368 |
| Age at diagnosis | 0.294 | 0.091-0.956 | .042 | 0.905 | 0.210-3.900 | .894 |
| Gender | 1.143 | 0.700-1.865 | .593 | 1.284 | 0.737-2.234 | .377 |
| Tumor stage | 1.008 | 1.003-1.013 | <. 001 | 1.008 | 1.003-1.013 | .001 |
Note: Bold values indicat statistically significant results.
Abbreviations: ACC, adrenocortical carcinoma; CI, confidence interval; HR, hazard ratio; IRG, immune-associated gene; TCGA, The Cancer Genome Atlas.
illustrates the expression of the three IRGs. The KM survival curve of all the patients with ACC in the training cohort is shown in Figure 3D. We analyzed the survival time and found that patients in the low- risk group had a longer survival time than those in the high-risk group (p <. 001). As shown in Figure 3E, by plotting time-dependent ROC curves of the three-IRG prognostic signature and other clinical features in the training cohort, we identified that the area under the ROC curve (AUC) value was 0.923. We also evaluated the age, three-IRG prognostic signature, gender, and tumor stage using multivariate Cox regression analy- ses, as shown in Figure 3F. From the risk score
(HR = 2.585, 95% CI = 1.569-4.258, p <. 001) of the three-IRG signature, a unique prognostic target of overall poor survival was identified (Table 3).
3.5 Authentication of the three-IRG prognostic signature in the entire TCGA ACC cohort and the testing cohort |
We utilize the three-IRG signature to evaluate the OS of patients with ACC in the testing cohort to prove the stability and predictive capacity of the three-IRG signature. Risk scores of the three-IRG signature in 38 patients with
(A)
Survival time (years)
(D)
Risk
High risk
Low risk
1
00
.
Dead
.
Alive
1.00
OD
4
~
Survival probability
0.75
0
10
20
30
(B)
Patients (increasing risk socre)
0.50
High risk
0.25
¥
.
Risk score
.
low Risk
p=1.5450-03
S
0.00
N
0
1
2
3
4
5
6
7
8
9
10
Time(years)
0
10
20
30
(C)
Patients (increasing risk socre)
Risk
High risk
19
17
11
8
2
2
0
0
0
0
0
type
type
Low risk
19
19
17
13
11
10
6
4
3
1
0
1.5
high
INHBA
low
0
1
2
3
4
5
6
7
8
9
10
Time(years)
HELLS
1
(E)
HDAC4
0.5
(F)
pvalue
Hazard ratio
1.0
0.8
Age
0.044
1.067(1.002-1.137)
True positive rate
0.6
Gender
0.837
0.847(0.175-4.103)
+
0.4
Stage
0.016
3.434(1.264-9.329)
0.2
risk score (AUC=0.867)
Age (AUC=0.534)
riskScore
0.004
2.168(1.281-3.668)
0.0
Gender (AUC=0.504)
Stage (AUC=0.811)
0
2
4
6
8
Hazard ratio
0.0
0.2
0.4
0.6
0.8
1.0
False positive rate
ACC from the testing cohort were computed using the previously mentioned formula. The 38 ACC patients were divided into the high-risk and low-risk groups (n = 19 per group) based on the median value. Figure 4A-B provides the clinical features of the testing cohort, the survival time, survival status, and the risk score. The expression of the three IRGs is shown in Figure 4C.
Figure 4D displays the KM survival curve of the testing cohort. Our results found that patients of the low- risk group subsisted for a longer time compared with the high-risk group (p = . 002). Figure 4E shows the plotting time-dependent ROC curves of the three-IRG prognostic signature and other clinical features in the training cohort, which demonstrated that the AUC of the three- IRG prognostic signature was distinctly better than other clinical features (risk value was 0.867). The results between the training cohort and the testing cohort were the same. From the multivariate Cox regression analysis, the risk score (HR = 2.168, 95% CI = 1.281-3.668, p =. 004) confirmed the three-IRG signature as a unique prognostic target (Figure 4F and Table 3).
All the patients in the TCGA ACC cohort were analyzed using time-dependent ROC curves, KM survival curve, and multivariate Cox regression analysis. The results are shown in Figure 5 and Table 3.
3.6 External validation using the GSE10927 and GSE19750 cohorts |
To externally validate the three-IRG signature in the GSE10927 (n=24) and GSE19750 (n=22) cohorts, we applied the time-dependent ROC curve (Table 4), univariate and multivariate Cox regression analyses, and KM survival curve (Table 3). The KM survival curve and time-dependent ROC curve for the GSE10927 cohort are shown in Figure 6A,B. The results confirmed that patients of the low-risk group subsisted for a longer time compared to those of the high-risk group (p = . 005). The AUC of three-IRG prognostic signature was distinctly better than other clinical features (risk value was 0.861). Similar results were found in the GSE19750 cohort (Figure 6C,D).
(A)
(D)
Risk
High risk
Low risk
Survival time (years)
~
·
Dead
.
Alive
1.00
8
N 2 4 6
0.75
0
Survival probability
0
20
40
60
80
(B)
Patients (increasing risk socre)
0.50
00
·
High risk
Risk score
+
-
low Risk
0.25
p=1.385e-07
2
0
0.00
0
1
2
3
4
5
6
7
8
9
10
11
12
0
20
40
60
80
Time(years)
(C)
Patients (increasing risk socre)
Risk
High risk
38
34
20
16
7
5
2
2
2
2
2
1
1
type
Low risk
39
39
34
24
20
17
11
8
5
3
2
1
1
type
1.5
high
0
1
2
3
4
5
6
7
8
9
10
11
12
INHBA
low
(E)
Time(years)
HELLS
1
HDAC4
0.5
1.0
(F)
pvalue
Hazard ratio
0.8
Age
0.236
1.019(0.987-1.052)
True positive rate
0.6
Gender
0.747
0.863(0.353-2.112)
0.4
Stage
0.031
1.707(1.049-2.780)
0.2
risk score (AUC=0.914)
Age (AUC=0.534)
0.0
Gender (AUC=0.504)
Stage (AUC=0.811)
riskScore
<0.001
1.842(1.471-2.307)
0.0
0.5
1.0
1.5
2.0
2.5
0.0
0.2
0.4
0.6
0.8
1.0
Hazard ratio
False positive rate
| Variables | Entire TCGA ACC cohort | GSE10927cohort | GSE19750 cohort |
|---|---|---|---|
| Three-IRG signature risk score | |||
| 1-Year | 0.883 | 0.853 | 0.792 |
| 3-Year | 0.923 | 0.956 | 0.897 |
| 5-Year | 0.914 | 0.861 | 0.765 |
| Age at diagnosis | |||
| 1-Year | 0.607 | 0.537 | 0.607 |
| 3-Year | 0.534 | 0.536 | 0.643 |
| 5-Year | 0.534 | 0.512 | 0.717 |
| Gender | |||
| 1-Year | 0.512 | 0.694 | 0.512 |
| 3-Year | 0.504 | 0.548 | 0.480 |
| 5-Year | 0.504 | 0.611 | 0.333 |
| Tumor stage | |||
| 1-Year | 0.611 | 0.616 | 0.611 |
| 3-Year | 0.811 | 0.791 | 0.676 |
| 5-Year | 0.811 | 0.438 | 0.488 |
Abbreviations: ACC, adrenocortical carcinoma; AUC, area under the ROC curve; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.
(A)
GSE19027 cohort
(C)
GSE19750 cohort
Risk
High
Low
Risk
+
High
+
Low
1.00
P=0.005
1.00
P=0.035
Survival probability
0.75
Survival probability
0.75
0.50
0.50
0.25
0.25
0.00
0.00
0
2
4
6
8
10
12
0
2
4
6
8
10
12
14
16
18
20
Time(years)
Time(years)
Risk
High Low
12
12
2
2
0
0
0
0
Risk
High
11
11
8
1
3
3
1
1
3
1
N
1
2
0
1
0
1
LOW
5
3
O
0
2
4
6
8
10
12
0
2
4
6
8
10
12
14
16
18
20
(B)
Time(years)
(D)
Time(years)
1.0
1.0
0.8
0.8
True positive rate
0.6
True positive rate
0.6
0.4
0.4
0.2
risk score (AUC=0.861)
0.2
Age (AUC=0.512)
risk score (AUC=0.765)
Gender (AUC=0.611)
Age (AUC=0.717)
0.0
Stage (AUC=0.438)
0.0
Gender (AUC=0.333)
Stage (AUC=0.488)
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
False positive rate
False positive rate
3.7 Survival, clinical features, and the three-IRG prognostic signature |
We analyzed the entire TCGA ACC cohort by stratified survival analysis, to further verify the prognostic capacity and explore the extensive feasibility of the three-IRG signature (Figure 7A-G). The samples with Stage I-II (p =. 011), Stage III (p =. 008), and Stage IV (p =. 031) were assigned to the low-risk group, which had better survival than the high-risk group. Similar results were found in diverse ages and different gender.
Clinical features were distributed along with the three-IRG prognostic signature. The high-risk group had a higher number of mortalities, which was reflective of poorer survival in these patients who had a high-risk value (p <. 05) (Figure 8B-D). Furthermore, there were different risk scores for patients in different disease
stages. Patients at the inchoate stages (Stage I and Stage II) have lower risk values than those at more advanced stages (Stage III and Stage IV; p <. 05; Figure 8E-G), which verified that risk scores of the IRG prognostic signature were distinctly related to the development of ACC.
3.8 PCA and functional enrichment analysis |
To determine the disparate distribution patterns between the high-risk and low-risk groups based on the three IRGs, we performed a PCA analysis. The two groups were obviously separated in two different directions based on the three IRGs. It was evident that the samples in the high-risk group were significantly different from
(A)
Stage I and Stage II
(B)
Stage III
(C)
Stage IV
2
High risk
10
High risk
:
High risk
Low risk
Low risk
Low risk
0
g
Survival rate
0
P=0.011
Survival rate
P=0.008
Survival rate
P=0.031
0.6
9%
0.6
0.4
0.4
0.4
2
2
2
:
:
:
0
1000
2000
3000
4000
0
1000
2000
3000
4000
0
500
1000
1500
Time (day)
Time (day)
Time (day)
(D)
Age =< 60
(E)
Age>60
(F)
Male
(G)
Female
2
High risk
2
High risk
9
2
Low risk
Low risk
High risk
High risk
Low risk
Low risk
8
:
:
:
Survival rate
Survival rate
P=0.033
P<0.001
Survival rate
Survival rate
6
0.6
0.6
P<0.001
0.6
P<0.001
0.4
0.4
3
3
=
8
2
2
2
8
8
8
0
1000
2000
3000
4000
0
1000
2000
3000
4000
0
1000
2000
3000
0
1000
2000
3000
4000
Time (day)
Time (day)
Time (day)
Time (day)
(A)
Risk
Risk
Risk Score
3
Survival Status ***
High
Stage ***
2
Low
Gender
Risk Score 6
Age
1
0
INHBA
-1
-2
Survival Status ***
-2
Alive
Dead
-3
Stage ***
HELLS
Stage I
Stage II
Stage III
Stage IV
Gender
Female Male
HDAC4
Age
=< 60
>60
(B)
(C)
(D)
Relative risk score of signature
Entire TCGA ACC cohort
GSE19027
10-
Relative risk score of signature
22-
Relative risk score of signature
GSE19750
P<0.001
P=0.005
15-
P=0.012
20-
5-
10-
18-
0
5-
16-
-5
Alive
14
Alive
0
Dead
Dead
Dead
Alive
(E)
(F)
(G)
Entire TCGA ACC cohort
Relative risk score of signature
10-
P=0.045
Relative risk score of signature
GSE19027
GSE19750
P=0.049
20-
Relative risk score of signature
P=0.004
12-
P<0.001
5-
18-
10-
8-
0-
16-
6-
-5
14
Stage I
Stage II
Stage III
Stage IV
Stage I
Stage II
Stage III
Stage IV
Stage I
Stage II
Stage III
Stage IV
Open Access
PC2
(A)
Entire TCGA ACC cohort
(B)
Training cohort
(C)
Testing cohort Low risk
Low risk
High risk
Low risk
High risk
High risk
2
·
19
10
12
0
0
2
0
8
8
PC3
:
PC3
8
PC3
8
0
-0.5
-0.5
2
PC2
-1.0
₪
PC2
N
:
3
1
2
-
40
0
0
10
1
-1
-1.5
1
0
A
-2
-2
-3
0
-1
-2
0
2
-3
0
-3
-2
-1
1
3
-3
-2
=1
0
1
2
3
S
-3
-2
-1
0
1
2
3
(D)
PC1
PC1
PC1
(E)
hsa04218
hsa04060
-log (adj p-value)
BP
CC
MF
hsa04110
hsa04933
20
0
hsa04512
hsa05205
15
0
008
10
hsa04672
hsa04510
5
o
0
GO:0140014
GO:0060485
GO:0030198
GO:0043062
GO:0048871
GO:0031012
GO:0062023
GO:0071556
GO:0098553
GO:0042611
GO:0005201
GO:0032395
GO:0042605
GO:0033218
GO:0042277
hsa04061
hsa04657
30
30
800
hsa04974
hsa04612
hsa04151
hsa04514
hsa04145
z-score
z-score
logFC
decreasing
increasing
downregulated
upregulated
decreasing
increasing
that of the low-risk group (Figure 9A-C). From the entire TCGA ACC cohort, 664 upregulated and 563 down- regulated DEGs were revealed between low-risk and high-risk groups. We identified their capacity using the KEGG and GO functional enrichment analyses. The top 15 terms of GO are listed in Figure 9D and Table 5. The top 15 terms of KEGG are shown in Figure 9E and Table 6. We found that they were significantly enriched in the extracellular matrix regulation, cell cycle, and phosphatidylinositol 3-kinase (PI3K)-Akt signaling pathway.
3.9 Survival analysis and correlation analysis of INHBA, HELLS, and HDAC4 |
We further explored the prognostic value of the three IRGs (INHBA, HELLS, and HDAC4) for construction of the
signature and identified that higher expression levels of INHBA, HELLS, and HDAC4 were related to worse OS time among patients with ACC (p <. 05; Figure 10A-C). Higher expression levels of HELLS and HDAC4 also indicated worse disease-free survival (Figure 10D-F). INHBA, HELLS, and HDAC4 in the TCGA ACC data set were analyzed by pairwise correlation analysis. The expression levels of HELLS and HDAC4 demonstrated an obvious positive correlation (p <. 001 and R =0.648). An increase in HELLS was related to an increase in HDAC4 (Figure 10G).
3.10 Correlation between tumor- infiltrating immune cells and the three-IRG signature |
We compared the tumor-infiltrating immune cells between low-risk and high-risk groups in the entire
| ID | Term | Category | Adjusted p |
|---|---|---|---|
| GO:0030198 | Extracellular matrix organization | Biological process | 2.92E - 15 |
| GO:0043062 | Extracellular structure organization | Biological Process | 2.92E -15 |
| GO:0060485 | Mesenchyme development | Biological process | 6.13E -09 |
| GO:0048871 | Multicellular organismal homeostasis | Biological process | 1.40E -07 |
| GO:0140014 | Mitotic nuclear division | Biological process | 4.74E -07 |
| GO:0031012 | Extracellular matrix | Cellular component | 6.31E - 22 |
| GO:0062023 | Collagen-containing extracellular matrix | Cellular component | 2.08E -20 |
| GO:0042611 | MHC protein complex | Cellular component | 1.56E -16 |
| GO:0071556 | Integral component of lumenal side of endoplasmic reticulum membrane | Cellular component | 1.09E -10 |
| GO:0098553 | Lumenal side of endoplasmic reticulum membrane | Cellular component | 1.09E -10 |
| GO:0005201 | Extracellular matrix structural constituent | Molecular function | 4.23E - 12 |
| GO:0042605 | Peptide antigen binding | Molecular function | 2.80E - 08 |
| GO:0033218 | Amide binding | Molecular function | 1.52E -07 |
| GO:0042277 | Peptide binding | Molecular function | 1.67E -07 |
| GO:0032395 | MHC class II receptor activity | Molecular function | 1.75E-06 |
Abbreviations: ACC, adrenocortical carcinoma; GO, Gene Ontology; MHC, major histocompatibility complex; TCGA, The Cancer Genome Atlas.
| ID | Term | Adjusted p |
|---|---|---|
| hsa04061 | Viral protein interaction with cytokine and cytokine receptor | 4.67E- 08 |
| hsa04672 | Intestinal immune network for IgA production | 3.63E - 07 |
| hsa04512 | ECM-receptor interaction | 3.63E - 07 |
| hsa04110 | Cell cycle | 3.77E - 07 |
| hsa04218 | Cellular senescence | 5.61E - 07 |
| hsa04060 | Cytokine-cytokine receptor interaction | 1.97E - 05 |
| hsa04933 | AGE-RAGE signaling pathway in diabetic complications | 1.97E - 05 |
| hsa05205 | Proteoglycans in cancer | 2.23E - 05 |
| hsa04151 | PI3K-Akt signaling pathway | 7.38E - 05 |
| hsa04974 | Protein digestion and absorption | 1.09E - 04 |
| hsa04612 | Antigen processing and presentation | 1.41E - 06 |
| hsa04514 | CAMs | 3.53E - 06 |
| hsa04145 | Phagosome | 6.89E-06 |
| hsa04510 | Focal adhesion | 1.65E -04 |
| hsa04657 | IL-17 signaling pathway | 4.21E - 04 |
Abbreviations: ACC, adrenocortical carcinoma; AGE-RAGE, advanced glycationend products-receptor of AGES; CAMS, Cell adhesion molecules; ECM, extracellular matrix; IL-17, interleukin-17; KEGG, Kyoto Encyclopedia of Genes and Genome; PI3K, phosphatidylinositol 3-kinase; TCGA, The Cancer Genome Atlas.
TCGA ACC cohort, to investigate the relationship between tumor immune microenvironment and the three-IRG prognostic signature. T cells, natural killer (NK) cells, macrophage cells, mast cells, and myeloid dendritic cells showed disparate abundance in the low- risk and high-risk groups (p <. 05; Figure 11). The correlation heat map was also plotted based on the distribution of immune cells in patients with ACC (Figure 12).
4 DISCUSSION |
Accumulating evidence indicated that aberrantly ex- pressed genes in tumors may be potential biomarkers for diagnosis, target therapy, and prognosis through bio- informatics analysis.26 In cancer treatment, immuno- therapy is gaining more and more attention. Currently, we can identify potential prognostic targets, analyze underlying mechanisms, and explore valuable IRGs through high-throughput sequencing data.27,28 We re- searched IRGs using the TCGA ACC data set and constructed a three-IRG prognostic signature that may be used as a prognostic indicator or immune therapeutic target in patients with ACC.
The immune system impacts the tumor micro- environment, particularly tumor immune escape.26 Immune components have quite a significant influence on the clinical outcomes and gene expression by tumor
(A)
Overall Survival
(B)
Overall Survival
(C)
Overall Survival
1.0
Low INHBA TPM
1.0
1.0
High INHBA TPM
Low HELL’S TPM
High HELLS TPM
Low. HDAC4 TPM
P=0.025
P<0.001
High HDAC4 TPM
P<0.001
0.8
0.8
0.8
Percent survival
Percent survival
Percent survival
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
0.0
0
50
100
150
0
50
100
150
0
50
100
150
Months
Months
Months
(D)
Disease Free Survival
(E)
Disease Free Survival
(F)
Disease Free Survival
0
Low INHBA TPM
a
High INHBA
Low HELLS TPM
1.0
High
Low HDAC4 TPM
P=0.970
P<0.001
High.HDAC4.T.P.M
P<0.001
0.8
0.8
0.8
Percent survival
Percent survival
Percent survival
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
0.0
0
50
100
150
0
50
100
150
0
50
100
150
(G)
Months
(H)
Months
(1)
Months
4
R = 0.648
R = 0.100
..
R = 0.165
HELLS Expression Level (log2 TPM)
P < 0.001
INHBA Expression Level (log2 TPM)
2.0
P = 0.381
2
INHBA Expression Level (log2 TPM)
P = 0.146
3
1.5
2
1.0
1
1
0.5
0
0
0.0
HDAC4 Expression Level (log2 TPM)
1
2
3
0
HELLS Expression Level (log2 TPM)
1
2
3
4
HDAC4 Expression Level (log2 TPM)
1
2
3
tissues, which exist in the tumor microenvironment.9,29,30 There were 3519 upregulated and 2090 downregulated DEGs found in ACC samples and normal adrenal samples in this study. GSEA suggested that these DEGs are involved in the humoral immune response (GO:0006959), and the immune effector process (GO:0002252) and immune response (GO:0006955), which indicated that the development of ACC is significantly correlated with immune responses and the immune system (Figure 2).
We extracted 332 IRGs and constructed an IRG signature based on three IRGs (INHBA, HELLS, and HDAC4) in the training cohort. Based on the median value of the three-IRG signature risk score, we divided the samples into the low-risk and high-risk groups. We discovered that the prognosis of samples could be distinguished by the three-IRG signature. The OS of the
high-risk group was worse than those of the low-risk group (Figure 3D). The multivariate Cox regression analysis suggested that the risk score may be regarded as a unique prognostic target of OS through multivariate Cox regression analyses and the three-IRG signature was superior to other clinicopathologic features according to the ROC curve. The above findings were verified the testing cohort, and the whole TCGA ACC, GSE10927, and GSE19750 cohorts.
We then explored the wide applicability of the three-IRG signature. Survival analyses were con- ducted in disparate subgroups, stratified by age, tumor stage, and gender. We found that the three- IRG signature was applicable in various subgroups (Figure 7). Moreover, we observed that the risk score of patients with early-stage ACC (Stage II and Stage I)
(A)
(B)
(C)
NK cell resting P=0.005
(D)
(E)
T cell CD4+ memory resting
T cell gamma delta
NK cell activated
P=0.049
Macrophage MO
0.5-
P=0.028
0.08-
P=0.037
0.15-
0.15-
0.6-
P<0.001
0.4-
0.06-
0.10-
0.10-
0.4-
0.3-
0.04-
0.05-
0.05-
0.2-
0.2-
0.02-
0.1-
0.00-
0.00-
0.00-
0.0-
…**
0.0
low risk (n=39)
high risk (n=38)
-0.02-
low risk (n=39)
high risk (n=38)
-0.05
low risk (n=39)
high risk (n=38)
-0.05
low risk (n=39)
high risk (n=38)
-0.2
low risk (n=39)
high risk (n=38)
(F)
Macrophage M1 P=0.024
(G)
Macrophage M2 P<0.001
(H)
(1)
(J)
Myeloid dendritic cell resting
Myeloid dendritic cell activated
Mast cell activated P<0.001
0.20-
0.8-
0.15-
P=0.002
0.5-
P<0.001
0.25-
0.15-
0.6-
0.10-
0.4-
0.20-
0.10-
0.3-
0.15-
0.4-
0.05-
0.2-
0.10-
0.05-
0.00-
0.2
0.00-
0.1-
0.05-
0.0-
0.00-
…
-0.05
low risk (n=39)
high risk (n=38)
0.0
low risk (n=39)
high risk (n=38)
-0.05
low risk (n=39)
high risk (n=38)
-0.1
low risk (n=39)
high risk (n=38)
-0.05
low risk (n=39)
high risk (n=38)
(K)
Type
8
Type
B cell naive
High risk
B cell memory
6
Low risk
B cell plasma
T cell CD8+
4
T cell CD4+ naive
T cell CD4+ memory resting
2
T cell CD4+ memory activated
T cell follicular helper
0
T cell regulatory (Tregs)
T cell gamma delta
-2
NK cell resting
NK cell activated
Monocyte
Macrophage MO
Macrophage M1
Macrophage M2
Myeloid dendritic cell resting
Myeloid dendritic cell activated
Mast cell activated
Mast cell resting
Eosinophil
Neutrophil
was lower than the advanced one (Stage IV and Stage III; Figure 8). Our findings demonstrated that the three-IRG prognostic signature is not only a potential predictive indicator of tumor progression but can be applied to all patients with ACC.
Functional enrichment analyses demonstrated that the DEGs between low-risk and high-risk groups were obviously enriched in cell cycle function and PI3K-Akt signaling pathway. Cell cycle regulation is known to play a vital role in the differentiation and growth of tumor cells.31 Subramanian et al.32 identified a series of genes involved in regulating pathways of the cell cycle and overexpression of these genes was correlated with poorer OS in patients with ACC. LncRNA-HOTAIR participated in ACC development by regulating the cell cycle and proliferating ACC cells.33 PI3K-Akt signaling pathway is a key driver in carcinogenesis and Akt overactivation has been verified in various endocrine gland neoplasms, including thyroid carcinoma subtypes, parathyroid
carcinoma, pituitary tumor, and pheochromocytoma.34 In vitro experiments revealed that low expression of SCTR could stimulate proliferation via the PI3K/AKT signaling cascade in ACC cells.35
Recently, with the advent of checkpoint inhibitors and targeted immunotherapy, a greater focus has been placed on tumor immunotherapy for patients with ACC.36,37 Our study evaluated the prognostic value of three IRGs for the construction of the IRG signature. Higher expressions of INHBA, HELLS, and HDAC4 were correlated with worse OS time for patients with ACC. Hence, the three IRGs may be promising therapeutic targets for ACC.
INHBA encodes a member of the transforming growth factor-@ super-family of proteins.38 INHBA participates in the regulation of the immune system, especially in skin suffering from repetitive ultraviolet-B irradiation.39 Mesenchymal stromal cells inhibited monocytes through upregulating INHBA in patients
Macrophage M1
T cell CD4+ memory activated
B cell plasma
T cell gamma delta
T cell CD8+
T cell follicular helper
T cell regulatory (Tregs)
Myeloid dendritic cell resting
T cell CD4+ memory resting
B cell naive
T cell CD4+ naive
Monocyte
NK cell activated
Mast cell resting
Macrophage M2
Mast cell activated
Myeloid dendritic cell activated
B cell memory
Eosinophil
Neutrophil
NK cell resting
Macrophage MO
1
B cell plasma
1
0.21
0.31
0.18
0.2
0.26
0.03
0.02
0.02
-0.12
-0.01
0
-0.33
-0.11
-0.06
-0.06
-0.13
-0.09
0.02
-0.01
-0.18
0.02
T cell gamma delta
0.21
1
0.2
0.17
0.05
0.05
-0.01
-0.1
0.2
-0.11
0.27
0.01
-0.16
0.09
0.23
-0.1
-0.24
-0.09
-0.06
-0.01
-0.27
-0.17
T cell CD8+
0.31
0.2
1
0.5
0.21
0.28
-0.16
-0.23
0.09
-0.05
-0.01
0.05
-0.41
-0.07
0
0.05
-0.38
-0.15
0.3
0.03 -0.23
-0.06
0.8
Macrophage M1
0.18
0.17
0.5
1
0.2
0.21
-0.16
-0.1 0.03
-0.28
-0.15
-0.1
-0.25
0.18
0.01
0
-0.4
-0.18
0.55
0.11
-0.12
-0.11
T cell follicular helper
0.2
0.05
0.21
0.2
1
0.36
0.21
-0.13
-0.06
-0.14
0.04 0.15
-0.37
-0.32
-0.26
0.01
0.32
0.14
0.2
-0.01
-0.04
0.05
0.6
T cell regulatory (Tregs)
0.26
0.05
0.28
0.21
0.36
1
-0.06
-0.15
-0.05
-0.09 -0.22
0.08
-0.51
-0.35
-0.02
-0.12
0.16
0.01
0.05
0.03
0.2
0.38
B cell naive
0.03
-0.01
-0.16
-0.16
0.21
-0.06
1
0.05
-0.12
-0.03
0.13
-0.08
-0.03
-0.26
-0.27
-0.22
0.22
-0.12
-0.14
-0.03
-0.07
0.01
0.4
T cell CD4+ naive
0.02
-0.1
-0.23
-0.1
-0.13
-0.15
0.05
1
-0.09
-0.11
-0.06 -0.1
-0.06
-0.06
0.05
-0.01
0.08
0 -0.06
-0.05
0.02
0.03
Myeloid dendritic cell resting
0.02
0.2
0.09
0.03
-0.06
-0.05
-0.12
-0.09 1
0.06
0.04
0.1
-0.14
-0.01
0.1
-0.04
-0.03
-0.12
0.04
-0.05
-0.08
-0.03
0.2
Monocyte
-0.12
-0.11
-0.05
-0.28
-0.14
-0.09
-0.03
-0.11 0.06
1
0.08
0.19
0.03
-0.13
-0.29
-0.11
0.02
0.1
-0.15
-0.13
0.09
-0.2
NK cell activated
-0.01
0.27
-0.01
-0.15
0.04
-0.22
0.13
-0.06
0.04
0.08
1
0.29
0
-0.04
0.03
-0.13
-0.04
0.06
-0.15
0.04
-0.51
-0.23
0
Mast cell resting
0
0.01
0.05
-0.1
0.15
0.08
-0.08
-0.1
0.1
0.19
0.29
1
-0.15
-0.21
-0.37
0.03
0.08
0.27
0.05
-0.02
-0.13
0.06
T cell CD4+ memory resting
-0.33
-0.16
-0.41
-0.25
-0.37
-0.51
-0.03
-0.06
-0.14
0.03
0
-0.15
1
0
0.04
0.09
-0.1
-0.09
-0.3
-0.16
0.07
-0.27
Macrophage M2
-0.11
0.09
-0.07
0.18
-0.32
-0.35
-0.26
-0.06
-0.01
-0.13
-0.04
-0.21
0
1
0.28 -0.12
-0.4
-0.11
0.08
0.03
-0.35
-0.44
-0.2
Mast cell activated
-0.06
0.23
0
0.01
-0.26
-0.02
-0.27
0.05
0.1
-0.29
0.03
-0.37
0.04
0.28
1
-0.05
-0.26
-0.1
-0.2
-0.13
-0.16
-0.18
B cell memory
-0.06
-0.1
0.05
0
0.01
-0.12
-0.22
-0.01
-0.04
-0.11
-0.13
0.03
0.09
-0.12
-0.05
1
0.13
0.12
0.2
-0.03
0.25
0.05
-0.4
Myeloid dendritic cell activated
-0.13
-0.24
-0.38
-0.4
0.32
0.16
0.22
0.08
-0.03
0.02
-0.04
0.08
-0.1
-0.4
-0.26
0.13
1
0.25 -0.17
-0.01
0.16
0.05
Eosinophil
-0.09
-0.09
-0.15
-0.18
0.14
0.01
-0.12
0
-0.12
0.1
0.06
0.27
-0.09
-0.11
-0.1
0.12
0.25
1
-0.1
-0.04
-0.15
0
-0.6
T cell CD4+ memory activated
0.02
-0.06
0.3
0.55
0.2
0.05
-0.14
-0.06
0.04
-0.15
-0.15
0.05
-0.3
0.08
-0.2
0.2
-0.17
-0.1
1
0.56
0.15
0.05
Neutrophil
-0.01
-0.01
0.03
0.11
-0.01
0.03
-0.03
-0.05
-0.05
-0.13
0.04
-0.02
-0.16
0.03
-0.13
-0.03
-0.01
-0.04
0.56
1
-0.03
0.14
-0.8
NK cell resting
-0.18
-0.27
-0.23
-0.12
-0.04
0.2
-0.07
0.02
-0.08
0.09
-0.51
-0.13
0.07
-0.35
-0.16
0.25
0.16
-0.15
0.15
-0.03
1
0.47
Macrophage MO
0.02
-0.17
-0.06
-0.11
0.05
0.38
0.01
0.03
-0.03
-0.2
-0.23
0.06
-0.27
-0.44
-0.18
0.05
0.05
0
0.05
0.14
0.47
1
-1
with myelodysplastic syndrome.4º IDO inhibitor induced the upregulation of IHNBA in regulating the functions of NK cells among sarcoma patients.41
HELLS belongs to the SNF2-family, such as chromatin remodeling proteins.42 HELLS mutation leads to immune deficiency syndrome in children.43 A study based on adult mice suggested that although HELLS is not indispensable for lymphoid development, it is still necessary for the proliferation of peripheral T lymphocytes.44 In addition, high-dose-rate y-irradiation could impair immune signaling pathways by inducing the downregulation of HELLS.45
HDAC4 is one of the histone-modifying enzymes, which are the main epigenetic regulators in the control of inflammatory processes. Furthermore, overexpres- sion of HDAC4 may suppress the production of Type I
interferons induced by pattern-recognition receptors in innate immunity.46 Aberrant epithelial responses may arise via a self-defense mechanism generated by miR-22 through suppressing HDAC4 expression in epithelial cells.47 Moreover, the HDAC4/PGRN and HDAC4/ nuclear factor-KB axes have been involved in the regulation of inflammatory cytokines in rheumatoid arthritis. 48
To identify gene signatures and further create clinical predictive models, high-throughput sequencing data has been utilized in large-scale research studies. Xu et al. stated that a signature could forecast the prognosis of ACC according to survival-associated alternative splicing events.49 We constructed an immune-related signature based on 3 IRGs in the TCGA ACC data set and further validated its application in two GEO data
sets. Furthermore, the predictive capacity of the signa- ture was further confirmed. Besides, we also explored the association between tumor-infiltrating immune cells and the three-IRG signature and we found that T cells, NK cells, macrophage cells, myeloid dendritic cells and mast cells showed diverse abundance between the low- risk and high-risk groups, which further verified the immune correlation of this signature.14
Although the three IRGs were regarded as therapeutic targets for the treatment of ACC, the underlining mechanisms have not been clarified, which was one limitation of this study. Therefore, more in vitro and in vivo experiments based on a larger sample size are necessary for validating these findings and applying our results into clinical immunotherapy practice in the future.
5 CONCLUSION |
We identified an immune-related signature based on 3 IRGs in ACC. The three-IRG signature has important clinical significance, not only as a good classifier in distinct subgroups of ACC, but also as a unique prognostic target for ACC. Inhibitors of the three novel IRGs (INHBA, HELLS, and HDAC4) might activate immune cells in ACC immune microenvironment and play a synergistic role in combination therapy with immunotherapy for ACC in the future.
AUTHOR CONTRIBUTIONS
Chengdang Xu, Caipeng Qin, Yun Peng, Jingang Jian, Chengdang Xu, and Xinan Wang performed the research. Chengdang Xu, Jingang Jian, and Yuxuan Song designed the research study. Chengdang Xu, Yuxuan Song, Jingang Jian, Xi Chen, and Xinan Wang analyzed the data. Chengdang Xu and Denglong Wu contributed essential tools. Chengdang Xu, Yun Peng, Jingang Jian, Yuxuan Song, Xi Chen, Xinan Wang, and Denglong Wu wrote the paper.
ACKNOWLEDGMENT
This study was supported by grants from New Frontier Technology Joint Research Project of Shanghai Municipal Hospital ( No. SHDC12019112 ) , Clinical Research Plan of SHDC (No. SHDC2020CR3074B), and Natural Science Foundation of Shanghai Municipal Science and Technology Committee (Nos. 21ZR1458300, 19ZR1448700, 22ZR1456800, and 20Y11904400).
CONFLICT OF INTEREST
The authors declare no conflict of interest.
DATA AVAILABILITY STATEMENT
All data were extracted from public databases: TCGA ACC data set (http://portal.gdc.cancer.gov/), GTEx data- base (http://www.gtexportal.org/home/index.html), and GEO database (http://www.ncbi.nlm.nih.gov/geo/).
ETHICS STATEMENT
This study does not involve relevant ethical issues.
ORCID
Chengdang Xu [ http://orcid.org/0000-0001-9253-2458
REFERENCES
1. Else T, Kim AC, Sabolch A, et al. Adrenocortical carcinoma. Endocr Rev. 2014;35:282-326.
2. Fassnacht M, Dekkers OM, Else T, et al. European Society of Endocrinology Clinical Practice Guidelines on the manage- ment of adrenocortical carcinoma in adults, in collaboration with the European Network for the study of adrenal tumors. Eur J Endocrinol. 2018;179:G1-G46.
3. Fassnacht M, Libé R, Kroiss M, Allolio B. Adrenocortical carcinoma: a clinician’s update. Nat Rev Endocrinol. 2011;7: 323-335.
4. Puglisi S, Perotti P, Cosentini D, et al. Decision-making for adrenocortical carcinoma: surgical, systemic, and endocrine management options. Expert Rev Anticancer Ther. 2018;18: 1125-1133.
5. Bellantone R, Ferrante A, Boscherini M, et al. Role of reoperation in recurrence of adrenal cortical carcinoma: results from 188 cases collected in the Italian National Registry for adrenal cortical carcinoma. Surgery. 1997;122:1212-1218.
6. Gonzalez RJ, Tamm EP, Ng C, et al. Response to mitotane predicts outcome in patients with recurrent adrenal cortical carcinoma. Surgery. 2007;142:867-875.
7. Yu WD, Wang H, He QF, Xu Y, Wang XC. Long noncoding RNAs in cancer-immunity cycle. J Cell Physiol. 2018;233: 6518-6523.
8. Pagès F, Berger A, Camus M, et al. Effector memory T cells, early metastasis, and survival in colorectal cancer. N Engl J Med. 2006;353:2654-2666.
9. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
10. Martinez-Bosch N, Vinaixa J, Navarro P. Immune evasion in pancreatic cancer: from mechanisms to therapy. Cancers. 2018;10:6.
11. Kim R, Emi M, Tanabe K. Cancer immunoediting from immune surveillance to immune escape. Immunology. 2007; 121:1-14.
12. Peng Y, Song Y, Ding J, Li N, Zhang Z, Wang H. Identification of immune-related biomarkers in adrenocortical carcinoma: immune-related biomarkers for ACC. Int Immunopharmacol. 2020;88:106930.
13. Jin Y, Wang Z, He D, et al. Analysis of m6A-related signatures in the tumor immune microenvironment and identification of clinical prognostic regulators in adrenocortical carcinoma. Front Immunol. 2021;12:637933.
14. Huang R, Liu Z, Tian T, et al. The construction and analysis of tumor-infiltrating immune cells and ceRNA networks in metastatic adrenal cortical carcinoma. Biosci Rep. 2020; 40:BSR20200049.
15. Anders S, McCarthy DJ, Chen Y, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013;8:1765-1786.
16. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Biocon- ductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009;26:139-140.
17. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for inter- preting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005;102:15545-15550.
18. Bao ZS, Li MY, Wang JY, et al. Prognostic value of a nine-gene signature in glioma patients based on mRNA expression profiling. CNS Neurosci Ther. 2014;20:112-118.
19. Cai J, Zhang W, Yang P, et al. Identification of a 6-cytokine prognostic signature in patients with primary glioblastoma harboring M2 microglia/macrophage phenotype relevance. PLoS One. 2015;10:e126022.
20. Zhang CB, Zhu P, Yang P, et al. Identification of high risk anaplastic gliomas by a diagnostic and prognostic signature derived from mRNA expression profiling. Oncotarget. 2015;6: 36643-36651.
21. Lossos IS, Czerwinski DK, Alizadeh AA, et al. Prediction of survival in diffuse large-B-cell lymphoma based on the expression of six genes. N Engl J Med. 2004;350:1828-1837.
22. Zhang W, Zhang J, Yan W, et al. Whole-genome microRNA expression profiling identifies a 5-microRNA signature as a prognostic biomarker in Chinese patients with primary glioblastoma multiforme. Cancer. 2013;119:814-824.
23. Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package forcomparing biological themes among gene clusters. OMICS. 2012;16:284-287.
24. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453-457.
25. 5. Newman AM, Steen CB, Liu CL, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37:773-782.
26. Cao J, Yang X, Li J, et al. Screening and identifying immune- related cells and genes in the tumor microenvironment of bladder urothelial carcinoma: based on TCGA database and bioinformatics. Front Oncol. 2020;9:1533.
27. Wang W, Zhao Z, Yang F, et al. An immune-related IncRNA signature for patients with anaplastic gliomas. JNO. 2018;136: 263-271.
28. Wei C, Liang Q, Li X, et al. Bioinformatics profiling utilized a nine immune-related long noncoding RNA signature as a prognostic target for pancreatic cancer. J Cell Biochem. 2019;120:14916-14927.
29. Cooper LA, Gutman DA, Chisolm C, et al. The tumor microenvironment strongly impacts master transcriptional regulators and gene expression class of glioblastoma. Am J Pathol. 2012;180:2108-2119.
30. Şenbabaoğlu Y, Gejman RS, Winer AG, et al. Tumor immune microenvironment characterization in clear cell renal cell
carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17: 231.
31. Pereira SS, Monteiro MP, Bourdeau I, Lacroix A, Pignatelli D. Mechanisms of endocrinology: cell cycle regulation in adrenocortical carcinoma. Eur J Endocrinol. 2018;179:R95-R110.
32. Subramanian C, Cohen MS. Over expression of DNA damage and cell cycle dependent proteins are associated with poor survival in patients with adrenocortical carcinoma. Surgery. 2019;165:202-210.
33. Yan ZC, He L, Qiu JH, et al. LncRNA HOTAIR participates in the development and progression of adrenocortical carcinoma via regulating cell cycle. Eur Rev Med Pharmacol Sci. 2018;22: 6640-6649.
34. Robbins HL, Hague A. The PI3K/Akt pathway in tumors of endocrine tissues. Front Endocrinol. 2016;6:188.
35. Lee M, Waser B, Reubi JC, Pellegata NS. Secretin receptor promotes the proliferation of endocrine tumor cells via the PI3K/AKT pathway. Mol Endocrinol. 2012;26:1394-1405.
36. . Fiorentini C, Grisanti S, Cosentini D, et al. Molecular drivers of potential immunotherapy failure in adrenocortical carci- noma. J Oncol. 2019;2019:1-7.
37. Jasim S, Habra MA. Management of adrenocortical carci- noma. Curr Oncol Rep. 2019;21:20.
38. Gaddy-Kurten D, Tsuchida K, Vale W. Activins and the receptor serine kinase superfamily. Recent Prog Horm Res. 1995;50:109-129.
39. Chen W, Zhang J. Potential molecular characteristics in situ in response to repetitive UVB irradiation. Diagn Pathol. 2016;11: 129.
40. Sarhan D, Wang J, Sunil Arvindam U, et al. Mesenchymal stromal cells shape the MDS microenvironment by inducing suppressive monocytes that dampen NK cell function. JCI Insight. 2020;5:e130155.
41. Nafia I, Toulmonde M, Bortolotto D, et al. IDO targeting in sarcoma: biological and clinical implications. Front Immunol. 2020;11:274.
42. Geiman TM, Durum SK, Muegge K. Characterization of gene expression, genomic structure, and chromosomal localization of Hells (Lsh). Genomics. 1998;54:477-483.
43. Thijssen PE, Ito Y, Grillo G, et al. Mutations in CDCA7 and HELLS cause immunodeficiency-centromeric instability- facial anomalies syndrome. Nat Commun. 2015;6:7870.
44. Geiman TM, Muegge K. Lsh, an SNF2/helicase family member, is required for proliferation of mature T lympho- cytes. Proc Natl Acad Sci USA. 2000;97:4772-4777.
45. Shin SC, Lee KM, Kang YM, et al. Differential expression of immune-associated cancer regulatory genes in low- versus high-dose-rate irradiated AKR/J mice. Genomics. 2011;97: 358-363.
46. Yang Q, Tang J, Pei R, et al. Host HDAC4 regulates the antiviral response by inhibiting the phosphorylation of IRF3. J Mol Cell Biol. 2019;11:158-169.
47. Moheimani F, Koops J, Williams T, et al. Influenza A virus infection dysregulates the expression of microRNA-22 and its targets; CD147 and HDAC4, in epithelium of asthmatics. Respir Res. 2018;19:145.
48. Shao L, Hou C. miR-138 activates NF-KB signaling and PGRN to promote rheumatoid arthritis via regulating HDAC4. Biochem Biophys Res Commun. 2019;519: 166-171.
49. Xu N, Ke ZB, Lin XD, et al. Identification of survival- associated alternative splicing events and signatures in adrenocortical carcinoma based on TCGA SpliceSeq data. Aging. 2020;12:4996-5009.
How to cite this article: Xu C, Qin C, Jian J, et al. Identification of an immune-related gene signature as a prognostic target and the immune microenvironment for adrenocortical carcinoma. Immun Inflamm Dis. 2022;10:e680. doi:10.1002/iid3.680