The gene expression data and clinical information of 79 patients with ACC were obtained from The Cancer Genome Atlas (TCGA) on June 13, 2021 (https://portal.gdc.cancer.gov/repository). The gene expression data of 258 normal adrenal gland samples were obtained from the Genotype-Tissue Expression (GTEx) database (https://toil.xenahubs.net/download/gtex_RSEM_gene_tpm.gz). The GSE19750 and GSE10927 datasets from the Gene Expression Omnibus (GEO) database were used as validation cohorts (https://www.ncbi.nlm.nih.gov/geo/). The “removeBatchEffect” function of the R package “limma” (version 4.1.0) was used to remove the batch effects of GSE19750 and GSE10927. The detailed flow chart of our study was displayed in Figure 1.

Seventy-nine patients with ACC were clustered into three groups based on hypoxia-related gene expression profiles using the R software package “ConsensusClusterPlus” (version 1.54.0). The maximum number of clusters is 6, and 80% of the total sample is drawn 100 times, clusterAlg = “hc,” innerLinkage = “ward.D2.” Clustering heatmaps were created using the R software package “pheatmap” (version 1.0.12). The Kaplan–Meier survival analysis was then used to compare the survival difference among three groups using the R software package “survival” (version 3.2-13) and “survminer” (version 0.4.9). The expression distribution of ferroptosis-related genes in three hypoxia-related subtypes was analyzed by the Wilcoxon test. The box plot was implemented by the R software package “ggplot2.” Ferroptosis-related gene expression heat map was implemented by the R software package “pheatmap.”

The R software package “Limma” (version 4.1.0) was used for differential gene analysis by comparing the gene expression pattern in ACC tissues to that in normal tissues. The adjusted P-value was analyzed to correct for false positive results in TCGA or GTEx. “Adjusted P < 0.05 and Log (Fold Change) >1 or Log (Fold Change) < −1” were defined as the thresholds for identifying the differentially expressed genes (DEGs). The volcano plot of DEGs was implemented by the R software package “ggplot2.” The hierarchical clustering heat map was implemented by the R software package “pheatmap.”

Seventy-five hypoxia-related genes (Supplementary Table S1) and 24 ferroptosis-related genes (Supplementary Table S2) were derived from published literature (17, 20). Subsequently, 34 hypoxia-related DEGs and 12 ferroptosis-related DEGs were extracted. For Kaplan–Meier curves, P values and hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests and univariate Cox proportional hazards regression using the R software package “survival” (version 3.2-13) and “survminer” (version 0.4.9). Prognosis-related genes with P < 0.05 were screened for the establishment of a hypoxia and ferroptosis-related gene signature.

A protein–protein interaction (PPI) network was constructed using the web tool called “STRING” (https://string-db.org/) to explore the interactions among the above key genes that affected overall survival (OS).

The key genes that affected OS were further screened with the least absolute shrinkage and selection operation (LASSO) Cox regression model, using the R software package “glmnet” (version 4.1-1). The risk score was calculated by combining the regression coefficients and the corresponding gene expression levels. Based on the median risk score, patients with ACC were divided into high- and low-risk groups. The Kaplan–Meier survival curves were generated using the R software package “survival” (version 3.2-13) and “survminer” (version 0.4.9). Receiver operator characteristic (ROC) curve analysis was performed to test the sensitivity and specificity of the predictive model using the R software package “timeROC” (version 0.4). Multivariate Cox regression analysis was performed to evaluate whether the risk score was an independent prognostic factor for OS. Age, gender, stage, and the risk score were used to construct a nomogram. The GSE19750 and the GSE10927 dataset were utilized as the validation cohort. The Kaplan–Meier survival analysis and the ROC curve analysis were repeated in the same way as before.

Tumor mutation burden (TMB) is the total amount of substitutions, insertions, or deletions per mega bases in the exon-coding regions in tumors and is a potential biomarker for predicting response to immune checkpoint inhibitor therapy (21). The expression distribution of 5 hub genes in ACC and normal adrenal gland tissues was analyzed by the student's t-test. The Sankey diagram was drawn with each column representing a characteristic variable, different colors representing different types, and lines representing the distribution of the same sample in different characteristic variables. Then, Spearman's correlation analysis was used to explore the correlation of the gene expression and TMB.

Based on the median risk score, patients with ACC were divided into high- and low-risk groups. Differential gene analysis was performed between the two groups. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the R software package “clusterProfiler” (version 4.0.5). The R software package “GSVA” was used to analyze the correlation between genes and pathways (version 1.42.0).

The gene expression data of 33 kinds of tumors were obtained from the TCGA database. Spearman's correlation analysis was used to explore the correlation of the gene expression and TMB in 33 kinds of tumors.

Statistical analyses were conducted through R language (version 4.0.3). Differences between the two groups were analyzed via the student's t-test. Differences among the three groups were analyzed via the Wilcoxon test. The statistical methods related to bioinformatic analysis are detailed in the corresponding sections. P < 0.05 was considered statistically significant.

The clinicopathologic features of 79 patients in the TCGA-ACC cohort are shown in Table 1. Our study classified these 79 ACC samples into molecular subtypes based on the expression profiles of hypoxia-related genes using the ConsensusClusterPlus package. The classification was reliable and stable at k = 3; therefore, the samples were divided into C1, C2, and C3 subtypes (Figures 2A–D). Patients with the C3 subtype had shorter survivals than those with the C1 and C2 subtypes (p = 0.0013; Figure 2E). Therefore, hypoxia-related gene expression can be used to characterize three ACC subtypes with distinct clinical outcomes.

We analyzed the associations between the three hypoxia-related subtypes and the expression of ferroptosis-related genes in ACC. We found that 17 of 25 ferroptosis related genes, including acyl-CoA synthetase long-chain family member 4 (ACSL4), atlastin GTPase 1 (ATL1), ATP synthase membrane subunit C locus 3 (ATP5MC3), cysteinyl tRNA synthetase 1 (CARS1), CDGSH iron sulfur domain 1 (CISD1), citrate synthase (CS), dipeptidyl-dipeptidase-4 (DPP4), Fanconi anemia complementation group D2 (FANCD2), farnesyl-diphosphate farnesyltransferase 1 (FDFT1), heat shock protein family A member 5 (HSPA5), heat shock protein beta 1 (HSPB1), lysophosphatidylcholine acyltransferase 3 (LPCAT3), nuclear receptor coactivator 4 (NCOA4), nuclear factor, erythroid 2 like 2 (NFE2L2), solute carrier family 1 member 5 (SLC1A5), solute carrier family 7 member 11 (SLC7A11), and transferrin receptor (TFRC), were significantly differentially expressed in the three hypoxia-related subtypes (Figure 3A). Among these ferroptosis-related genes, the expressions of FDFT1 and ATP5MC3 were significantly higher in C3 samples than in C1 and C2 samples. By contrast, the expression of HSPB1 was significantly lower in C3 samples than in C1 and C2 samples. The heatmap showed similar results (Figure 3B). Therefore, these three hypoxia-related subtypes had a strong association with ferroptosis in ACC.

Using the limma package for analysis of differential gene expression in 79 ACC samples and 258 normal adrenal gland samples (Adjusted P < 0.05; |log2FoldChange| > 1), 4,743 DEGs were obtained, including 4,029 upregulated and 714 downregulated genes (Figures 4A,B). Next, 34 hypoxia-related DEGs and 12 ferroptosis-related DEGs were separately extracted from 75 previously identified hypoxia-related genes and 24 previously identified ferroptosis-related genes. After further screening these key genes by log-rank tests and univariate Cox proportional hazards regression, we found that 12 genes were significantly correlated with OS in patients with ACC (Figure 4C). A PPI network of 12 genes was then retrieved from String to explore the interaction among these genes (Figure 4D). We then constructed a LASSO Cox regression model (Figures 5A,B). The risk score was calculated based on the results of the LASSO Cox regression model, as follows: (−0.3843) ^* ACSL4 + (0.9418) ^* FANCD2 + (−0.1097) ^* HIF3A + (0.3538) ^* HSPA5 + (0.5395) ^* PSMB7 (Figure 5C). Seventy-nine patients with ACC were then separated into high- and low-risk groups based on the median risk score. Our results showed that patients in the high-risk group exhibited shorter OS than those in the low-risk group (P < 0.001; Figure 5D). We further verified the predictive efficacy of this model using ROC curve analysis. The AUC of the model for predicting OS at 3 years was 0.925 (Figure 5E), suggesting that the genetic model had accurate and robust prognostic prediction performance. To further explored the five hub genes (ACSL4, FANCD2, HIF3A, HSPA5, and PSMB7), we found that the expression of five hub genes in ACC was significantly different from that in the normal tissues (Figure 6A). Figure 6B represented the relationship among the expression of hub genes, the tumor stages, and survival status. We further found that each gene was significantly associated with TMB (Figure 6C).

Univariate Cox regression analysis showed that the risk score calculated based on our predictive model was significantly associated with poor prognosis in patients with ACC (P < 0.001, HR: 3.57, 95% CI: 2.39–5.32; Figure 7A). Multivariate Cox regression analysis suggested that the risk score was an independent risk factor for poor prognosis in patients with ACC (p < 0.001, HR: 4.21, 95% CI: 2.45–7.24; Figure 7B). A nomogram was established to predict the OS of ACC patients. Age, sex, stage, and the risk score were used to construct the nomogram. The risk score contributed the most to the prediction of OS (Figure 7C). To verify the generalizability of the predictive model, its predictive performance was externally validated using pooled data from the GSE19750 and GSE10927 datasets. Consistent with the results of the TCGA dataset, the risk score was significantly associated with OS in these datasets (P < 0.001; Figure 7D). ROC curve analysis confirmed the reliability and stability of the predictive model for predicting OS (AUC = 0.827; Figure 7E). The heatmap displays the difference in the expression profiles of the five prognostic genes in the low-risk and high-risk groups (Figure 7F).

To explore the biological functions and pathways that were related to the risk score, 1,860 DEGs (1,379 up-regulated genes and 481 down-regulated genes) between the high- and low-risk groups were obtained (Figures 8A,B). GO enrichment and KEGG pathway analyses were performed. Figure 8C displayed the significantly enriched KEGG pathways, including cell cycle, DNA replication, basal cell carcinoma, oocyte meiosis, and cellular senescence for up-regulated DEGs, and chemokine signaling pathway, drug metabolism-cytochrome P450, hematopoietic cell lineage, cytokine-cytokine receptor interaction, and viral protein interaction with cytokine and cytokine receptor for down-regulated DEGs. Figure 8D displayed the significantly enriched GO items in up-regulated DEGs, including organelle fission, nuclear division, and chromosome segregation for biological process (BP), and chromosomal region, chromosome, centromeric region, and condensed chromosome for cellular component (CC), and microtubule binding, catalytic activity, and DNA helicase activity for molecular function (MF). Figure 8E displayed the significantly enriched GO items in down-regulated DEGs, including T cell activation, mononuclear cell differentiation, and cellular calcium ion homeostasis for BP, and external side of plasma membrane, secretory granule membrane, and specific granule for CC, and immune receptor activity, cytokine receptor activity, and C-C chemokine receptor activity for MF. Figure 8F represented the interactive relationship between enriched BP pathways. GSVA enrichment analysis strongly implicated cell cycle checkpoint, mitotic sister chromatic segregation, and nuclear division as the main biological process (Figure 8G). The chord diagram in Figure 8H displayed GO biological terms for the top 100 genes with the largest fold change. GO terms including pattern specification process, regionalization, and cell fate specification were enriched.

We further explored the effect of 5 hub genes (ACSL4, FANCD2, HIF3A, HSPA5, and PSMB7) in 33 kinds of cancer in the TCGA database (Figure 9). Expressions of five genes were all significantly correlated with TMB in ACC. Among the 33 cancer types, ACSL4 expression was most correlated with TMB in ACC, whereas FANCD2, HIF3A, HSPA5, and PSMB7 were most correlated with TMB in thymoma (THYM).