Multi-modal characterization of metabolic and immune gene clusters in adrenocortical carcinoma treatment

Establishment of the TIME subtype

According to recent literature, knowledge of the TIME should be the main emphasis of further ACC research⁸. Therefore, we analyzed 24 microenvironmental cell subpopulations using GSEA and classified ACC patients (TCGA + GSE76019 + GSE76021) into three clusters based on ssGSEA (Fig. 1a–c). The results were clear: TIME cluster A had the lowest ssGSEA scores, while TIME cluster C had the highest. In 23 microenvironmental cell subpopulations (excluding plasma cells), the three clusters exhibited statistically significant differences in scores. Patients in TIME cluster B had the best prognosis (Fig. 1d). Furthermore, four key immune checkpoints (PD-1, PD-L1, PD-L2, CTLA4) exhibited the lowest expression in TIME cluster A (Fig. 1e–h).

Establishment of subtypes based on DEGs

Our analysis of the differential genes between the TIME clusters identified 18 common DEGs (Fig. 2a). These 18 genes underwent unsupervised hierarchical clustering analysis to find that the clustering stability was optimum at k = 2 (Fig. 2b, Supplementary Fig. 2a–c). We validated the classification by repeating the clustering analysis on the original dataset and two independent GEO datasets (GSE33371 and GSE10927). The results align with the previous classification pattern (Fig. 2c and Supplementary Fig. 2d–f). We named these groups Gene clusters A and B and plotted the expression heatmap of the 18 differentially expressed genes (Fig. 2e). Single-cell RNA sequencing data showed that these 18 genes were highly expressed in specific cell subsets and associated with immunity and metabolism (Fig. 2f). Patients in the Gene cluster B group had a better prognosis (Fig. 2g). This group showed high infiltration abundance in 19 microenvironment cell subsets and significantly high expression of four key immune checkpoints (Fig. 2h, i). Finally, we calculated the Steroid-related Immune Score (SIS) by PCA. We separated the 142 ACC patients into groups with high and low SIS based on the optimal cutoff value (1.041494).

a 18 common DEGs were selected from the TIME clusters. b Consensus clustering matrix heatmap with DEG-related molecular pattern in TCGA + GSE76019 + GSE76021 when k = 2. c Consensus clustering matrix heatmap with DEG-related molecular pattern in TCGA + GSE76019 + GSE76021 + GSE33371 + GSE10927 when k = 2. d The principal component analysis plot showing the distribution of Gene clusters. e The heatmap showing the expression distribution of DEGs. f The heatmap showing the expression distribution of DEGs in each cell type in the single-cell RNA sequencing data. g Kaplan-Meier curves for predicting the overall survival of ACC patients in Gene clusters (log-rank test, P = 0.003). h Differential infiltration of 24 immune cells in Gene clusters. * P < 0.05, ** P < 0.01, *** P < 0.001, ns, no significance. i Differential expression of CTLA4, PD-L1, PD-L2, and PD-1 in ACC patients in Gene clusters. * P < 0.05, ** P < 0.01, *** P < 0.001, ns, no significance.

AI validation of SIS and associated pathological features

Deep neural networks are an indispensable tool for medical image analysis. They can identify features that are difficult to detect by the naked eye, such as microsatellite instability (MSI), tumor mutation burden (TMB), and gene expression status^9,10. This study used deep learning technology to analyze whole slide images (WSIs), verified the effectiveness of SIS grouping, and explored its pathological characteristics to understand better the uniqueness of ACC (Fig. 3a–d). We used two mainstream deep learning networks, ResNet50 and Vision Transformer-B16, for validation and comparison (Table 1). The five-fold cross-validation results demonstrated that the AUC of the SIS classification reached 0.8 ± 0.01, with ResNet50 performing best (AUC = 0.8214, accuracy = 0.71) (Fig. 3e, f). Furthermore, ResNet50 excelled in the binary classification C1A/B model prediction of ACC (AUC = 0.848, accuracy = 0.74) (Supplementary Fig. 3b, c), confirming the efficacy of SIS grouping in pathology. We created a heatmap in the original slice (Fig. 3g) to more clearly visualize the SIS grouping and its classification probability. We used Class Activation Maps (CAMs) technology to visualize the pathological features associated with the SIS subgroups. The findings indicated that the SIS correlated with sinusoidal invasion and necrosis in the Weiss score, and surprise, with lymphocytic infiltration (Fig. 3h and Supplementary Fig. 3d). To enhance the verification of the model’s clinical application, we picked the ResNet50 model exhibiting optimal performance. We validated ACC patients from the First Affiliated Hospital of Dalian Medical University and the Second Affiliated Hospital of Dalian Medical University. The findings indicated that patients in the high SIS group exhibited significant lymphocytic infiltration. Conversely, patients in the low SIS group had minimal lymphocytic infiltration, with both groups marked by sinusoidal invasion and necrosis (Fig. 3i, Supplementary Fig. 3e). This suggests that the model exhibited consistent performance in the external validation set and possesses potential therapeutic applications. Subsequent prognostic follow-up assessments indicated that patients in the high SIS group had significantly prolonged survival periods compared to those in the low SIS group. However, owing to the restricted follow-up duration and patient population, the p-value failed to achieve statistical significance (p > 0.05). Survival studies, after integrating these patients with those from the TCGA cohort, indicated a decrease in the p-value, implying that these patients align with the prognostic characteristics of the SIS subgroup (Supplementary Fig. 3f, g, Supplementary Fig. 4a).

a Segmentation and background subtraction of pathological images. b Color normalization of patches. c Training of a deep learning model using patches from two subgroups of patients in the training set. d Testing of the trained model on all patches of patients in the test set and statistical classification of the patients. e Five-fold cross-validation AUC plot of ResNet50. f Five-fold cross-validation AUC plot of Vision Transformer. g WSIs from high SIS and low SIS tumor patients in the TCGA test set reveal spatial patterns associated with SIS prediction. h Lymphocyte infiltration in WSIs was found to correlate with SIS subtypes using CAMs. i The extent of lymphocyte infiltration found in WSIs in the external validation set.

The new ACC subgroup SIS closely related to clinical outcomes

Several subtypes have been identified in the ACC dataset, including the COC and C1A/C1B subtypes(Fig. 4a)^11,12. The high SIS group was more likely to correlate with the positive prognosis COC1 and the inactive C1B subtype. Conversely, the low SIS group was significantly correlated with the unfavorable prognosis of COC2 and COC3, as well as the aggressive C1A subtype. There are significant differences in SIS between the COC subtype, C1A/C1B subtype, expression subtype, and methylation subtype (Fig. 4b, c and Supplementary Fig. 4d, e). The TCGA immune subtype analysis shows that ACC patients are mainly distributed in the C3 group with the best prognosis and the C4 group with the worst prognosis. The C3 group had the highest number of high SIS patients, while the C4 group had the highest number of low SIS patients. Patients in the C3 group had a SIS level much greater than those in the C4 group (Fig. 4d). Patients with high SIS have a favorable prognosis in the TCGA and GEO datasets (Fig. 4e, f, and Supplementary Fig. 4a-c). Univariate and multivariate Cox regression analysis clearly identified SIS as an independent prognostic factor for ACC patients (Fig. 4h, i). The low SIS group was composed mainly of patients with the following clinical parameters: female gender, T4, N1, M1, Stage III, and Stage IV (Fig. 4j). The Weiss score is crucial in the pathological diagnosis of ACC. The low SIS group accounted for over 50% of the Weiss score-related indicators, including necrosis, mitotic rate > 5/50 HPF, venous invasion, and sinusoidal (lymphatic) invasion (Fig. 4g). Epithelial-to-mesenchymal transition (EMT), tumor mutational burden (TMB), mRNA expression-based stemness index (mRNAsi), and Ki-67 are common biomarkers in oncology. ACC patients with high EMT, high TMB, high mRNAsi, and high Ki-67 generally have a poor prognosis, while patients in the low SIS group within these stratifications have the worst survival outcomes and vice versa (Fig. 4k-n). We identified mutations in known ACC genes (e.g., CTNNB1, ZNRF3)¹² and also found mutations in some new genes (e.g., TTN, HLTF, ADAMTS16, NSD1, and PARP8) (Fig. 4o). The only exception was HLTF, which exhibited a high mutation rate in the high SIS group, while the other genes were almost exclusively mutated in the low SIS group.

a The distribution of various ACC subtypes and pathological parameters in SIS subgroups. b Distribution of SIS in COC subtypes. c Distribution of SIS in C1A/C1B subtypes. d Distribution of SIS in Immune subtypes. e The Kaplan-Meier curve was used to predict the overall survival of the SIS (TCGA + GSE76019 + GSE76021)(log-rank test, P = 2.281e-04). f The Kaplan-Meier curve was used to predict the overall survival of the SIS (TCGA + GSE76019 + GSE76021 + GSE33371 + GSE10927) (log-rank test, P = 2.498e-04). g The distribution of SIS subgroups in the Weiss Scoring System is shown in a Likert plot. * P < 0.05, ** P < 0.01, *** P < 0.001. h, i Univariate and multivariate analysis of clinical characteristics. j Distribution of SIS subgroups across clinical characteristics. k-n Kaplan-Meier curves for the prediction of overall survival in ACC patients with Ki-67, EMT, TMB, and mRNAsi and their stratification in SIS. o Genomic alterations in SIS subgroups.

The low SIS group closely related to steroid synthesis

According to KEGG GSEA and GSVA analysis, the high SIS group exhibited enrichment in immune-related pathways, whereas the low SIS group had enrichment in steroid biosynthesis pathways. The HALLMARK GSEA and GSVA analyses indicated that the high SIS group had enrichment in immunological pathways, while the low SIS group demonstrated enrichment in cholesterol homeostasis pathways (Fig. 5a, b, Supplementary Fig. 6a, b). We speculate that patients with high SIS may be related to immune response, while patients with low SIS may be related to adrenal function. SIS was substantially negatively correlated with the adrenal cortical differentiation index (ADS) (r = −0.62) (Fig. 5e). Most patients with high SIS were predicted to have low ADS (71.4%) and did not have abnormal hormone (61.5%) or cortisol (80.8%) secretion. In contrast, most patients with low SIS were predicted to have high ADS (66%) and had abnormal hormone (78.7%) and cortisol (57.4%) secretion (Fig. 5d). In the cortisol present group, 84% of patients exhibited low SIS or high ADS. In the hormone present group, 79% of patients had low SIS, which is greater than the 70% of patients with high ADS. Therefore, low SIS is a more effective assessment of adrenal function than high ADS (Fig. 5f, g). We screened nine core genes related to steroid hormone synthesis by HALLMARK_CHOLESTEROL_HOMEOSTASIS and KEGG_STEROID_BIOSYNTHESIS analysis (Fig. 5c). Mitotane is the main drug to inhibit ACC steroid hormone synthesis, which mainly inhibits SOAT1, which is related to cholesterol storage, and CYP11A1 and CYP11B1, which are associated with the conversion of cholesterol to cortisol and aldosterone (Fig. 5m). Among the 12 core genes, DHCR7, CYP51A1, SOAT1, and CYP11A1 were highly correlated with SIS and ADS. DHCR7 was significantly and stably differentially expressed among them in the hormone-related subgroup, while the other genes were more variable (Fig. 5h). Except for EBP, CYP11A1, and CYP11B1, genomic alterations were more common in the low SIS group compared to the high SIS group, particularly for DHCR7, SOAT1, and FDFT1, which showed no alterations in the high SIS group but ≥10% alterations in the low SIS group (Supplementary Fig. 5a). Among these 12 genes, only patients with DHCR7 genomic alterations showed poor prognosis in overall survival (OS), progression-free survival (PFS), and disease-specific survival (DSS) (Supplementary Fig. 5b–d). In contrast, patients with SOAT1 genomic alterations showed poor prognosis only in PFS (Supplementary Fig. 5e–g). The metabolic genes most associated with SIS features, DHCR7 and SC5D, were further investigated by Lasso regression, Random Forest(RF), and Support Vector Machine – Recursive Feature Elimination(SVM-RFE) machine learning methods (Fig. 5i). In the pan-cancer expression analysis, DHCR7 showed the highest expression in ACC, while SC5D exhibited an intermediate expression level. Nevertheless, high expression of both was associated with poorer patient prognosis (Fig. 5j and Supplementary Fig. 5h–j). After extensive analysis, we concluded that DHCR7 is a potential molecular marker for ACC. Single-cell RNA sequencing results showed that DHCR7 was specifically expressed in cancer cells (Fig. 5l).IHC results showed elevated DHCR7 protein expression in tumor tissues (Fig. 5k). CCK-8 assay results showed that NCI-H295R cells (IC50 = 4.684 μM) were more sensitive to Mitotane than SW-13 cells (IC50 = 6.918 μM) (Fig. 5l, m). Mitotane sensitivity was significantly increased in both DHCR7 knockdown cells (Fig. 5n, o).

a, b KEGG and HALLMARK GSEA in SIS subgroups. c Common genes screened from the core genes of the two enrichment pathways. d The distribution of the ensemble of ADS, hormone, and cortisol levels among SIS subgroups. e Correlation between SIS and ADS. f, g SIS and ADS are distributed over hormones and cortisol, respectively. h Correlation of steroid hormone-related genes and Mitotane primary target genes with SIS and ADS expression. Expression values were compared with low SIS, high ADS, hormone present, and cortisol present. i Lasso regression, RF, and SVM-RFE machine learning jointly screened for genes most associated with SIS. j Kaplan-Meier curve was used to predict DHCR7 overall survival. k Protein expression of DHCR7 in ACC tissues and normal adrenal tissues by IHC. l Distribution of DHCR7 expression in single-cell RNA sequencing. m Mechanistic pathways associated with DHCR7 and Mitotane. n Cell viability was assessed after treatment with different concentrations of Mitotane in SW-13 cells. Control or DHCR7 siRNA was transfected, incubated for 48 h, and then collected for RT-PCR analysis of the DHCR7 gene. Cell viability was assessed after treatment of SW-13 cells with control or DHCR7 siRNA transfected with Mitotane (6.9 μM) for 48 h. o Cell viability was evaluated in NCI-H295R cells treated with varying concentrations of Mitotane, as well as in cells transfected with control or DHCR7 siRNA and incubated for 48 h, followed by RT-PCR analysis of DHCR7 expression. Cell viability was assessed after treatment of NCI-H295R cells with control, or DHCR7 siRNA transfected with Mitotane (4.7 μM) for 48 hours. * P < 0.05, ** P < 0.01, *** P < 0.001.

The high SIS group closely related to immunity

The findings from GSEA and GSVA indicated that patients in the high SIS group exhibited more pronounced tumor immune features Supplementary Fig. 6a, b). Therefore, we investigated the relationship of SIS with the ESTIMATE score and the abundance of immune cell infiltration more closely. The results showed that the high SIS group had a much stronger correlation with the ESTIMATE score than the low SIS group did within the same dataset. In fact, the high SIS group had a correlation that was above 0.8 in both datasets (Fig. 6a). Regarding the abundance of immune cell infiltration, both groups showed correlations with increased monocytes, decreased neutrophils, decreased basophils, decreased natural killer cells, decreased naive CD8 T cells, and increased cytotoxic cells. The abundance of CD8 T cells and Tfh cells in the high SIS group increased with increasing SIS. In contrast, in the low SIS group, the abundance of macrophages, macrophages M1, macrophages M2, and activated NK cells increased with lower SIS (Fig. 6b). Through machine learning analysis, pan-cancer immune studies revealed four immune-related factors—MHC, CP, EC, and SC¹³. ACC analysis revealed higher MHC and EC levels, and lower CP and SC levels in the high SIS group (Fig. 6c). Meanwhile, the high SIS group exhibited a significantly higher immunophenoscore (IPS) (Fig. 6d). TIP analysis revealed higher activity scores for CD8 T cells, macrophages, NK cells, and infiltration of immune cells into tumors (step 5) in the high SIS group (Fig. 6e). And the overall activity score strongly correlated with SIS (r = 0.69) (Fig. 6f). Analysis of the 5 immune expression signatures revealed higher, statistically significant scores in the high SIS group (Fig. 6g). In addition, immune infiltration is associated with DNA damage¹⁴. Our analysis revealed lower SNV neoantigens, nonsilent mutation rate, copy number variation (CNV) burden (number of segments), and homologous recombination deficiency (HRD) in the high SIS group, in addition to the group having a low proliferation rate (Supplementary Fig. 6f). Immunity quantitative trait loci (immunQTLs) were applied to assess the impact of genetic variants on immune infiltration. Survival-associated immunQTLs (FDR < 0.05) GWAS analyses revealed the highest number of QTLs for regulatory T cells (Tregs), followed by CD8 T cells, and the lowest number for resting CD4 memory T cells (Supplementary Fig. 6g). Among these GWAS disorders, the highest number of calcium level-related QTLs were all from Tregs. Schizophrenia involves T-cell regulatory (Tregs), T-cell CD8, macrophages M1, and T-cell CD4 memory resting, with the widest coverage. TIDE is an important indicator for predicting tumor response to immune checkpoint blockade (ICB). Both TCGA and GEO datasets showed that SIS was significantly negatively correlated with TIDE (Supplementary Fig. 6c, d). Over 70% of patients in the high SIS group responded to ICB, while more than 65% in the low SIS group were resistant. Notably, the low SIS group represented over 80% of ICB-resistant patients (Fig. 6h, i). TIDE’s two main immune escape mechanisms are associated with T cell dysfunction and T cell exclusion, respectively¹⁵. Our data indicate that the immune escape mechanism in ACC is primarily associated with T cell dysfunction. Cancer-associated fibroblasts (CAFs), myeloid-derived suppressor cells (MDSCs), and the M2 subtype of tumor-associated macrophages (TAMs) are the three main cell types that suppress T-cell infiltration¹⁶. A significant negative correlation was observed between the high SIS group and TAM M2. In addition, the immune profile of ACC was strongly correlated with interferon-gamma (IFNG) response and T-cell inflammatory phenotype (Merck18)(Supplementary Fig. 6e)¹⁷. We further assessed the expression of three molecules associated with tumor escape mechanisms in the low SIS group—Most MHC molecules were underexpressed in the low SIS group, thus avoiding T cell recognition; immunosuppressive factors (e.g., TGFBR1) might be upregulated for tumor escape; and immunostimulatory factors (e.g., RAET1E) might be downregulated to avoid immune attack (Fig. 6j).

We selected a low-SIS patient (with abnormal hormone secretion) for single-cell RNA sequencing analysis to explore the interaction between the tumor and its microenvironment. The analysis revealed the cell subpopulations of this patient were mainly concentrated in cancer cells (78.95%), macrophages (15.44%), progenitor cells (3.42%), mesenchymal cells (1.17%), T cells (0.65%), and endothelial cells (0.38%) (Fig. 6k, Supplementary Fig. 7). We screened 79 exosome-associated genes specifically expressed in tumor cells (avg_log2FC > 0.4), and KEGG showed that these genes were enriched in hormone synthesis and secretion pathways, in addition to metabolic reprogramming and energy metabolism pathways, immune evasion and TIME pathways, and drug metabolism and resistance pathways (Fig. 6l). Using HitPredict and CellphoneDB analyses, we found that these exosomes interacted with ligand receptors of immune cells, with more than 65% of the interactions involving T cells and macrophages. Most of these exosomes were highly expressed in the TCGA and GEO datasets in low SIS patients (Fig. 6m). These genes were strongly associated with tumor, extracellular matrix, and immune processes (Fig. 6n). Unexpectedly, some genes (e.g., AXL, HGF, PDGFB, and PDGFRB) were associated with EGFR tyrosine kinase inhibitor resistance.

a Correlation analysis between SIS subgroups and tumor microenvironment ESTIMATE score. b Heatmap of immune cell infiltration in different SIS groups. c Relationship between SIS and MHC, EC, SC, and CP. d Distribution and comparison of IPS in SIS subgroups. e Heatmap of immune activity scores of seven-step cancer-immunity cycle under SIS subgroup. f Correlation between SIS and Overall activity scores. g Heatmap of five characteristic immune expression signature scores under SIS subgroups, * P < 0.05, ** P < 0.01, *** P < 0.001. h Correlation analysis of SIS with TIDE and immunotherapy response in the TCGA cohort. i Correlation analysis of SIS with TIDE and immunotherapy response in the GEO cohort. j Expression patterns of immune-related genes in SIS subgroups. k UMAP analysis of cell type distribution. l KEGG enrichment pathway analysis of ACC-associated exosomes. m Interaction relationship between ACC exosomes and immune-related cell ligand receptors and their expression. n KEGG enrichment pathway analysis of exosomes and their associated ligand receptors.

Drug sensitivity prediction based on SIS subgroups

We successfully predicted drug sensitivity for SIS subgroups using the GDSC1 and GDSC2 datasets. The results showed that patients in the high SIS group were significantly more responsive to drugs targeting the PI3K/Akt/mTOR pathway (90%). Additionally, drugs targeting the following pathways showed greater sensitivity in this group: 54.5% on the Autophagy pathway, 40.9% on PI3K/Akt/mTOR pathway, and 18.2% on the Protein Tyrosine Kinase/RTK pathway. Furthermore, the GDSC1 and GDSC2 datasets showed that BI-2536 (PLK1 inhibitor) was more effective in the low SIS group based on the TCGA and GEO datasets (Fig. 7a). SPIED3 also predicted increased sensitivity to the MTOR inhibitor in the high SIS group, consistent with the GDSC results (Fig. 7b). CMap analysis, on the other hand, confidently identified several drugs with increased sensitivity in the low SIS group, including calmodulin antagonists^18,19, dopamine receptor antagonists²⁰, oxidosqualene cyclase inhibitors²¹, serotonin receptor antagonists²², and sterol demethylase inhibitors²³. All of these drugs inhibited steroid hormone synthesis, supporting our finding that the low SIS group exhibited adrenal function characteristics (Fig. 7c). SPIED3 and CMap together predicted four classes of drugs to be effective for patients in the low SIS group: acetylcholine receptor antagonist (mebeverine), dopamine receptor antagonist, BCR-ABL kinase inhibitor (imatinib), and norepinephrine reuptake inhibitor (maprotiline) (Fig. 7d). Notably, all dopamine receptor antagonists are used for the treatment of schizophrenia. Furthermore, single-cell RNA sequencing revealed that the tumor cell-specific genes identified were significantly enriched in the CALCIUM ION BINDING and IRON ION BINDING pathways, both of which are associated with schizophrenia, in the GO-MF enrichment analysis(Fig. 7e)^24,25. The CHP1 gene plays a role in calcium metabolism and facilitates ferroptosis. The PCLO gene, linked to calcium metabolism, shows a significant association with schizophrenia(Fig. 7f)²⁶.

a Heatmap showing GDSC1 and GDSC2 drugs with statistically significant differences in the SIS subgroups in the TCGA and GEO cohorts, and the pathways affected by these drugs are included. b Screening of CMap drugs consistently associated with TCGA and GEO from the top 200 SPIED3 drugs related to SIS. c Screening of Score < -90 or Score>90 in TCGA and GEO consistently related drugs from CMap. d SPIED3 and CMap identified four drugs that are potentially effective for the low SIS group. e GO-MF enrichment analysis of genes specifically expressed in tumor cells. f Heatmap of genes associated with iron metabolism & ferroptosis, calcium metabolism, and schizophrenia in single-cell sequencing results.