MDPI
Article
Comprehensive Characterization of the Regulatory Landscape of Adrenocortical Carcinoma: Novel Transcription Factors and Targets Associated with Prognosis
João C. D. Muzzi 1,2,3D, Jéssica M. Magno 2,3, Jean S. Souza 3, Larissa M. Alvarenga 10D, Juliana F. de Moura 1D, Bonald C. Figueiredo 3,40 and Mauro A. A. Castro 2,*
1 Laboratório de Imunoquímica (LIMQ), Pós-Graduação em Microbiologia, Parasitologia e Patologia, Departamento de Patologia Básica, Universidade Federal do Paraná (UFPR), Curitiba 81530-990, Brazil
2 Laboratório de Bioinformática e Biologia de Sistemas, Pós-Graduação em Bioinformática, Universidade Federal do Paraná (UFPR), Curitiba 81520-260, Brazil
3 Oncology Division, Instituto de Pesquisa Pelé Pequeno Príncipe, Curitiba 80250-060, Brazil
4 Molecular Oncology Laboratory, Centro de Genética Molecular e Pesquisa do Câncer em Crianças (CEGEMPAC), Curitiba 80030-110, Brazil
* Correspondence: mauro.castro@ufpr.br
check for updates
Citation: Muzzi, J.C.D .; Magno, J.M .; Souza, J.S .; Alvarenga, L.M .; de Moura, J.F .; Figueiredo, B.C .; Castro, M.A.A. Comprehensive Characterization of the Regulatory Landscape of Adrenocortical Carcinoma: Novel Transcription Factors and Targets Associated with Prognosis. Cancers 2022, 14, 5279. https://doi.org/10.3390/ cancers14215279
Academic Editor: Peter Igaz
Received: 13 August 2022
Accepted: 25 October 2022 Published: 27 October 2022
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- iations.
CC
İ
BY
Copyright: @ 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).
Simple Summary: Adult adrenocortical carcinoma (ACC) is a rare and aggressive tumor in adults, usually associated with excessive steroid secretion. It is highly metastatic and has few therapeutic options and a poor prognosis. Here, we explore the hallmarks of ACC influenced by transcription factors and their target genes (regulons) to provide a prognostic overview of ACC biology. Using an in silico clinical data analysis approach, we assessed human transcriptomic data from publicly available datasets. We found four distinct clusters of regulons associated with good and worse prognoses associated with cell proliferation and/or immunologic activity. Some findings require further bench analyses, primarily focusing on worse prognostic regulons and their targets.
Abstract: We reconstructed a transcriptional regulatory network for adrenocortical carcinoma (ACC) using transcriptomic and clinical data from The Cancer Genome Atlas (TCGA)-ACC cohort. We inves- tigated the association of transcriptional regulatory units (regulons) with overall survival, molecular phenotypes, and immune signatures. We annotated the ACC regulons with cancer hallmarks and assessed single sample regulon activities in the European Network for the Study of Adrenal Tumors (ENSAT) cohort. We found 369 regulons associated with overall survival and subdivided them into four clusters: RC1 and RC2, associated with good prognosis, and RC3 and RC4, associated with worse outcomes. The RC1 and RC3 regulons were highly correlated with the ‘Steroid Phenotype,’ while the RC2 and RC4 regulons were highly correlated with a molecular proliferation signature. We selected two regulons, NR5A1 (steroidogenic factor 1, SF-1) and CENPA (Centromeric Protein A), that were consistently associated with overall survival for further downstream analyses. The CENPA regulon was the primary regulator of MKI-67 (a marker of proliferation KI-67), while the NR5A1 regulon is a well-described transcription factor (TF) in ACC tumorigenesis. We also found that the ZBTB4 (Zinc finger and BTB domain-containing protein 4) regulon, which is negatively associated with CENPA in our transcriptional regulatory network, is also a druggable anti-tumorigenic TF. We anticipate that the ACC regulons may be used as a reference for further investigations concerning the complex molecular interactions in ACC tumors.
Keywords: adrenocortical carcinoma; transcriptomics; regulatory network; transcription factor; immuno-oncology; KI-67; NR5A1; CENPA; ZBTB4; IZKF1
1. Introduction
Adrenocortical carcinoma (ACC) is a rare, aggressive endocrine malignancy with a bimodal age distribution and distinct characteristics between pediatric and adult tu-
mors [1,2]. ACC is characterized by a highly proliferative and immune-suppressed tumor microenvironment, high production of corticoids, TP53 mutation, and an upregulation of the WNT/ -Catenin pathway [3-6].
Zheng et al. (2016) classified The Cancer Genome Atlas (TCGA)-ACC patients into four groups based on mRNA K4 clustering: Steroid Phenotype High (HSP), HSP + Proliferation, Steroid Phenotype Low (LSP), and LSP + Proliferation [7]. The Steroid Phenotype is related to the activation of steroid biosynthesis pathways, while Proliferation was assessed by a proliferation score proposed by Wirapati et al. (2008) [8]. This classification presented a high overlap with the C1A/C1B molecular classification presented by Reynies et al. (2009) [9], being HSP and HSP + Proliferation related to C1A and worse outcomes, and LSP associated with C1B. In our previous study [6], we have shown that LSP and LSP + Proliferation pre- sented a significant presence of immune infiltration compared to the immune-suppressed microenvironment of HSP and HSP + Proliferation. Furthermore, Landwehr et al. (2020) showed that excessive glucocorticoid levels, present in nearly 60% of ACC patients, are related to T cell depletion in the tumor microenvironment [5].
The functional interplay between the tumor and infiltrating immune cells within the tumor microenvironment provides insights into genes associated with the anti-tumor immune response [10-12]. The levels and distribution of CD3+ and CD8+ T cell infiltra- tion distinguish four solid tumor phenotypes: hot (or inflamed), altered excluded, altered immunosuppressed, and cold (or non-inflamed) [11-13]. ACC is described as an immuno- logically cold tumor, presenting one of the lowest immune infiltrates among 30 solid cancer types from TCGA [7,14]. However, the amount and efficacy of the immune infiltrate depend on pre-existing low levels of intratumoral steroids [5,6]. This suggests that in order to boost the anti-tumor immune response, it is necessary to eliminate excessive glucocorticoid levels, and that may be the main reason for heterogeneous results in five clinical trials using four different immune checkpoint inhibitors (avelumab, nivolumab, pembrolizumab, and ipilimumab) [15-19].
2022 WHO Classification of Adrenal Cortical Tumors recommends that diagnosti- cians specify the mitotic count and the nuclear protein Ki-67 (Ki-67) labeling index in all ACCs [20]. A Ki-67 (or MKi-67) labeling index relates to proliferation and malignancy, besides being used in the risk stratification and the rationale of adjuvant mitotane ther- apy [20-22]. Steroidogenic Factor 1 (SF1) immunohistochemistry is the most reliable and specific biomarker to confirm the adrenal cortical origin [23]. This transcription factor (TF) is encoded by the Nuclear Receptor Subfamily 5 Group 1 (NR5A1) gene, whose over- expression is associated with increased steroid metabolism, proliferation, and a worse outcome [24,25].
Different combinations of regulators and molecular factors may be associated with cancer development. The inference of regulatory networks helps to understand how these factors may be related, converging on cellular mechanisms, which can add to the understanding of the biology of the disease or intervention strategies [26-29]. To reconstruct a regulatory network, gene expression data can be used to evaluate mutual information between a TF and potential target genes, generating regulatory hubs called regulons [26]. This reverse-engineering method has been successfully applied in other cancer types (e.g., [26,27,29]).
In the present study, we inferred a transcriptional regulatory network for ACC using publicly available transcriptomic and clinical data from the TCGA-ACC cohort [7]. Through multivariate Cox analysis, we found 369 regulons, composed of a TF and its direct and indirect targets, relating to overall survival. We investigated how these regulons correlate with molecular phenotypes and immune signatures [14]. In addition, we annotated these regulons with Molecular Signatures Database (MSigDb) Hallmarks, representing well- defined biological states or processes [30]. Finally, we tested the prognostic value of the regulon activity in the European Network for the Study of Adrenal Tumors (ENSAT) cohort [31].
2. Methods
2.1. The Cancer Genome Atlas-ACC Data
The TCGA-ACC RNA-seq and clinical data [7] available from the GDC repository were downloaded using the TCGABiolinks package v.2.20.1 in R [32-34]. Next, we assessed the curated survival data of TCGA-ACC participants using the Xena Browser [35]. Then, the gene expression matrix was filtered using the AnnotationHub package v.3.0.2 in R [36] for protein-coding genes. Finally, we normalized the raw counts with the variance stabilizing transformation (VST) method from the DESeq2 package v.1.32.0 in R [37], using the Steroid Phenotype classification in the experimental design (see Section 2.6).
2.2. Regulon Inference
The normalized gene expression matrix was used to call regulons with the RTN package v.2.16.0 in R [26]. First, we reconstructed regulons for 1605 TFs [38,39] using the ARACNe algorithm [40]. Then, we used the tni.alpha.adjust() function [41] to define a statistical threshold that controls the trade-offs between Type I and Type II errors at a similar level described by Castro et al. (2016) [26], using a Benjamini-Hochberg [42] adjusted p-value cutoff of 0.05.
2.3. Regulon Activity
The regulon activity was estimated using a two-tailed gene set enrichment analysis (GSEA2) [26] which produces a differential enrichment score (dES) for each sample. A positive dES represents activated regulons, while a negative dES represents suppressed regulon activity. Values near zero indicate inconclusive activity. We selected regulons with a minimum of 15 positive and 15 negative targets [30] to assess regulon activity, which is regarded as the minimum gene set size for downstream enrichment analyses [30].
2.4. Survival Analysis
We used the regulon activity for multivariate Cox analysis [43] relating to overall survival (OS) and Progression-Free Interval (PFI) (the period from the date of diagnosis until loco-regional or systemic recurrence, second malignancy, or death from any cause) [35]. For this analysis, we considered the tumor stage and age as covariates using the RTNSurvival package v.1.20.0 pipeline [44], generating the hazard ratio (HR) and 95% confidence interval (CI) for each regulon. Since the Cox analysis assessed the time-to-event association between steroid-related regulons and OS, steroid-secreting status was not included in the Cox model due to potential co-linearity with steroid-related regulons. Therefore, the regulons with an adjusted p-value below 0.05 for OS were selected for the subsequent analysis. For Kaplan-Meier survival curves [45], samples were divided into high, inconclusive, or low regulon activity, and the p-values were calculated using log-rank statistics [46,47].
2.5. Clustering
The regulons’ activity dES was used for the unsupervised consensus clustering using the ConsensusClusterPlus package v.1.56.0 in R [48]. We chose k equal to four.
2.6. Steroid and Proliferation Classification
Of the 92 samples in the TCGA-ACC cohort, 79 had RNA-seq data, and 78 were listed in the mRNA K4 molecular classification [7], which assigns a Steroid Phenotype to the samples.
The mRNA K4 classified the participants in “steroid-phenotype-high” (n =25), “steroid- phenotype-high + proliferation” (n = 22), “steroid-phenotype-low” (n =27), and “steroid- phenotype-low + proliferation” (n = 4).
We separated the Steroid Phenotype and the Proliferation profiles, resulting in two groups for the Steroid Phenotype, HSP (n = 47) and LSP (n = 31), and two groups for the higher and lower proliferation scores (n = 26 and 52, respectively). The proliferation score
used by Zheng et al. (2016) [7] was based on a proliferation gene set signature described by Wirapati et al. (2008) [8].
2.7. Differential Expression for Steroid and Proliferation Phenotypes
We called differentially expressed genes (DEGs) for the Steroid and Proliferation phenotypes using the DESeq2 package v.1.32.0 in R [37]. To avoid confounding effects, the two profiles were independently analyzed, generating DEGs relating to the Steroid independent of Proliferation (IP) classification and the Proliferation independent of Steroid (IS) classification. We considered significant DEGs with adjusted p-values below 0.05 in the Wald test [37].
2.8. Transcriptional Network Analysis (TNA)
The association between the activity of a regulon and the Steroid and Proliferation phenotypes was assessed using GSEA2 implemented in the TNA pipeline [38]. Here, we used the list of DEGs described in Section 2.7. Associations with adjusted p-values below 0.01 were considered significant.
2.9. Functional Annotation with MSigDb Hallmarks
The prognostic regulons were annotated with MSigDb Hallmarks [30] using the tni.annotate.regulons() function from the RTN package in R [26]. We used the Hallmark gene sets from the msigdbr package v.7.4.1 [49].
2.10. Immune Correlation
Values for leukocyte fraction, immune signatures, T-cell receptor (TCR), and B-Cell Receptor (BCR) metrics were retrieved from the master table for TCGA-ACC participants from Thorsson et al. (2018) [14]. We calculated the Spearman correlation between these values and the regulon activity.
2.11. Duals Inference
To search for co-regulatory associations between pairs of prognostic regulons, we used the RTNduals pipeline as described in Chagas et al. (2019) [50].
2.12. ENSAT Cohort Data
The clinical data and the normalized gene expression matrix were assessed from the ENSAT cohort using the GEOquery package v.2.60.0 [51]. This cohort comprises 44 ACC samples and is available in the Gene Expression Omnibus (GEO) portal under the accession number GSE49278 [31]. When more than one entry from the gene expression data referred to the same gene symbol, we selected the one with the higher coefficient of variation between the samples. The additional clinical data was obtained from the supplementary tables of Assié et al. (2014) [31].
2.13. Regulon Activity and Survival Analysis in the ENSAT Cohort
The regulatory network inferred in the TCGA-ACC cohort was used to calculate the regulon activity in the ENSAT cohort using GSEA2 [26]. We selected the regulons related to OS in the TCGA-ACC for multivariate Cox analysis in the ENSAT cohort. The HR and 95% CI for OS were inferred using tumor stage and age as covariates. For the Kaplan-Meier survival analysis, we followed the same protocol described in Section 2.4.
2.14. Statistics and Visualization
The R packages available in CRAN [52] and Bioconductor [53] repositories were used for statistical analyses. All p-values were corrected for multiple hypotheses using the Benjamini-Hochberg correction [42], and if not specified otherwise, we considered the result significant when below 0.05.
The two-sided Mann-Whitney-Wilcoxon test [54,55] was used for the boxplot com- parison when only two pairs were available. For general comparison between more than two groups, we used the Kruskal-Wallis rank-sum test [56] followed by Dunn’s test [57] for multiple pairwise comparisons.
For the construction of heatmaps, we used the ComplexHeatmap package v.2.8.0 [58] in R. To visualize the regulon network, we used the RedeR package v.2.0.0 [59] in R.
3. Results
3.1. The Identification of 369 Regulons with Prognostic Values Related to Molecular Phenotypes and Leukocyte Fractions
The TCGA-ACC RNA-seq data was used to call regulons using the RTN R pack- age [26] (Supplementary Table S1). Of the 1605 TFs annotated in the Lambert et al. (2018) collection [39], 611 had at least 15 positive and 15 negative targets in our transcriptional regulatory network (Supplementary Table S2, Supplementary Figure S1), which is regarded as the minimum gene set size for downstream enrichment analyses [30]. In Figure 1A, we present the general workflow used in this study. Figure 1B illustrates an example of a regulon and its targets’ distribution along the karyogram.
A
B
4
4
TCGA ACC cohort (78 samples)
A
4
R
N
7
K
V
gene expression matrix (24233 protein coding genes )
D
4
4
D
x transcription factors (1605 )
A
A
1
4
7
D
7
1605 regulons
Step 1
Regulator
Negative
Regulon reconstruction
Target
Positive
924 regulons
Step 2
Regulons with more than 15 targets
chr 01
611 regulons
Step 3
chr 02
Regulons with more than 15 positive and negative targets
chr 03
369 regulons
Step 4
chr 04
Prognosis Regulons in Cox multivariate Analysis
· Consensus clustering
Regulon activity
· Survival analysis
· GSEA for molecular classificartion
. Correlation with leukocyte fraction
chr 22
· Funcional annotation
. Prognostic validation in ENSAT cohort
chr X
0 Mb
50 Mb
100 Mb
150 Mb
200 Mb
Using multivariate Cox analysis, we used the OS data to assess the association with the regulons’ activity. Of the 611 regulons, we found 369 related to OS, of which 330 (89.4%) were also related to PFI (Supplementary Tables S3 and S4). Figure 2A shows the activity profile of the 369 prognostic regulons in the TCGA-ACC cohort (Supplementary Table S5). In the unsupervised clustering, the 188 good prognostic regulons showed higher activity in the LSP and C1B participants, and in participants with a higher leukocyte fraction. Conversely, the 181 poor prognosis regulons showed higher activity in the HSP and C1A participants, and participants with a lower leukocyte fraction. Consistently, the activity of high- and low-risk regulons presented opposite correlations with the leukocyte fraction (Figure 2B, Supplementary Table S4) and with the Steroid IP and Proliferation IS scores (Figure 2C and 2D, respectively, Supplementary Tables S4 and S6).
3.2. Consensus Clustering Resulted in Four Regulon Clusters with Different Functional and Molecular Characteristics
We used consensus clustering to look for subgroups within the low- and high-risk regulons (Supplementary Figure S2). The 188 regulons with good prognosis were divided into two clusters: regulon cluster (RC) 1 and RC2, with 62 and 126 regulons, respectively. The 181 regulons related to a worse prognosis were divided into RC3 and RC4, with 113 and 68 regulons, respectively (Figure 3A). RC1 activity showed the highest positive correlation with the presence of immune infiltrate and with the Steroid IP score, as opposed to RC3 (Figure 3B,C). RC2 had the lowest scores for the Proliferation IS classification, in contrast to RC4 (Figure 3D). The regulons’ correlation with leukocyte fraction was negatively associated with the Steroid IP (p = - 0.94, p-value < 2.2 x 10-16) (Figure 3E) and the Proliferation IS scores (p = - 0.52, p-value < 2.2 x 10-16) (Figure 3F).
Figure 4A shows the regulon enrichment in the MSigDb Hallmarks (Supplementary Figure S3, Supplementary Table S7). RC3 and RC4, in contrast to RC1 and RC2, were positively enriched in proliferation Hallmarks such as MYC targets v1 and v2, E2F targets, Mitotic Spindle, G2M checkpoint, and WNT/ß-Catenin Signaling (Figure 4B). In addition, RC1 was activated, while RC3 was repressed, in the immune Hallmarks (i.e., IL6 JAK STAT3 Signaling, Interferon (IFN) y and «, and Inflammatory responses-Figure 4C, Complement, Allograft rejection, and Coagulation), in addition to Apoptosis, P53 pathway, and some immune-related signaling pathways (i.e., IL2 STAT5 signaling and TNF-« signaling via NFK-B). On the other hand, RC2 had the lowest, while RC4 had the highest scores for PI3K Akt mTOR (Figure 4D) and Hedgehog signaling.
We also evaluated the RCs concerning the six immune signatures and the TCR metrics in- ferred by Thorsson et al. (2018) [14] for the TCGA-ACC participants (Supplementary Figure S4A, Supplementary Table S4). RC1 had the highest and RC3 the lowest correlations for Macrophage Regulation and Lymphocyte Infiltrate Signature Scores. Moreover, RC1 and RC2 negatively correlated with Proliferation and Wound Healing signatures, as opposed to the positive values found in RC3 and RC4. The four clusters showed a weak correlation with the IFN-y Response and TGF-ß Response signatures. Regarding the adaptive immune response (Supplementary Figure S4B, Supplementary Table S4), RC1 related positively to TCR Shannon and TCR Richness, while RC3 related negatively to these features. The RC activities did not present a significant correlation with TCR Evenness. Concerning the B cell response, the BCR metrics inferred by Thorsson et al. (2018) [14] were available for only five samples, making comparisons with sufficient statistical power unfeasible.
A
TCGA-ACC cohort (n=78)
B
1.0
7
Tumor Stage
Leukocyte Fraction Spearman Correlation
Steroid Phenotype
0.5
Proliferation
C1A/C1B
Purity
0.0
Leukocyte Fraction
-0.5
-1.0
C
2
1
369 prognostic regulons
Steroid IP
Score
0
-1
-2
D
2
Proliferation IS
1
Score
0
-1
HR
Steroid IP
Prolif. IS
LF rho
Hazard Ratio
Steroid IP Score
Proliferation IS Score
Leuk. Fraction Correlation
Regulon Activity
Regulon Cluster
Cluster
-2
Low Risk Regulons
High Risk Regulons
<1
RC1
RC3
1
>1
-2
0
2
-2
0
2
-1
0
1
-2
0
2
RC2
RC4
Tumor Stage
Steroid
Proliferation
C1A.C1B
Leukocyte Fraction
Prognostic regulons
Regulons with HR <1 (n=188)
1
3
NA
Steroid_High
Low Proliferation
C1A
Purity
2
☒ 4
Steroid_Low
High Proliferation
C1B
Regulons with HR >1 (n=181)
0
0.5
1
0
0.2
0.4
A
Kruskal-Wallis p < 2.2 x 10-16
B
Kruskal-Wallis p < 2.2 x 10-16
Overall Survival Hazard Ratio
Leukocyte Fraction Correlation
ns
ns
5.00
1.0
3.00
2.00
0.5
1.00
0.0
0.50
-0.5
0.25
RC1
RC2
RC3
RC4
-1.0
RC1
RC2
RC3
RC4
C
Kruskal-Wallis p < 2.2 x 10-16
D
Kruskal-Wallis p < 2.2 x 10-16
2
2
Steroid IP score
Proliferation IS
1
score
1
0
0
-1
-1
-2
-2
RC1
RC2
RC3
RC4
RC1
RC2
RC3
RC4
E
2
F
2
R= - 0.94,
0< 2.2 x 10-16
R= - 0.52,
p < 2.2 × 10-16
1
1
Steroid IP score
Proliferation IS
0
score
0
-1
-1
-2
-2
-0.5
0.0
0.5
-0.5
0.0
0.5
Leukocyte Fraction correlation
Leukocyte Fraction correlation
Regulon Cluster
RC1 (n=62)
RC2 (n=126)
RC3 (n=113)
RC4 (n=68)
A
369 prognostic regulons
Hallmark Category
Regulon Cluster
B
RC1
RC2
RC3
RC4
cellular component
metabolic
Kruskal-Wallis p < 2.2 x 10-16
Hazard Ratio
development
pathway
WNT Beta Catenin Signaling
Steroid IP score
Proliferation IS score
DNA damage
proliferation
Leuk.Fraction rho
immune
signaling
2
MYC_TARGETS_V2
G2M_CHECKPOINT
1
proliferation
E2F_TARGETS
MYC_TARGETS_V1
0
MITOTIC_SPINDLE
P53_PATHWAY
-1
DNA damage
UV_RESPONSE_DN
DNA_REPAIR
UV_RESPONSE_UP
-2
cellular component
APICAL_SURFACE PEROXISOME
RC1
RC2
RC3
RC4
APICAL_JUNCTION
SPERMATOGENESIS PANCREAS_BETA_CELLS
C
Kruskal-Wallis p < 2.2 x 10-16
development
MYOGENESIS
EPITHELIAL_MESENCHYMAL_TRANSI
ANGIOGENESIS ADIPOGENESIS
Inflammatory Response
WNT_BETA_CATENIN_SIGNALING
2
MTORC1_SIGNALING
TGF_BETA_SIGNALING
1
NOTCH_SIGNALING
PI3K_AKT_MTOR_SIGNALING
signaling
HEDGEHOG_SIGNALING
0
KRAS_SIGNALING_DN
ESTROGEN_RESPONSE_EARLY KRAS_SIGNALING_UP IL2_STAT5_SIGNALING
1
TNFA_SIGNALING_VIA_NFKB
-2
ANDROGEN_RESPONSE
ESTROGEN_RESPONSE_LATE
RC1
RC2
RC3
RC4
UNFOLDED_PROTEIN_RESPONSE
pathway
PROTEIN_SECRETION HYPOXIA
D
Kruskal-Wallis p < 2.2 x 10-16
APOPTOSIS
*
REACTIVE_OXYGEN_SPECIES_PATHW
PI3K Akt mTOR Signaling
CHOLESTEROL_HOMEOSTASIS GLYCOLYSIS
metabolic
BILE_ACID_METABOLISM
2
HEME_METABOLISM
XENOBIOTIC_METABOLISM
1
”
FATTY_ACID_METABOLISM
OXIDATIVE_PHOSPHORYLATION
IL6_JAK_STAT3_SIGNALING
0
-
INTERFERON_GAMMA_RESPONSE
immune
INTERFERON_ALPHA_RESPONSE
INFLAMMATORY_RESPONSE
-1
COMPLEMENT
ALLOGRAFT_REJECTION
2
COAGULATION
RC1
RC2
RC3
RC4
Regulon Cluster
Hazard Ratio
Steroid IP Score
Proliferation IS Score
Leuk. Fraction Correlation
Regulon Enrichment
Regulon cluster
RC1 (n=62)
RC3 (n=113)
RC1
RC3
<1
RC2 (n=126)
RC4 (n=68)
RC2
RC4
>1
-2
0
2
-2
0
2
-1
0
1
-2
0
2
3.3. The ENSAT Cohort Confirmed the Prognostic Value of 89.5% of the Regulons Related to Survival on TCGA-ACC
We investigated the regulon activity in the ENSAT cohort (44 ACC samples) (Supplementary Figure S5, Supplementary Table S5) and confirmed the association with survival for the 369 regulons. Of the 369 regulons, 361 were present in the ENSAT gene expression matrix, of which 323 (89.5% of the 361) had significant HRs in multivariate Cox analysis for OS (Supplementary Table S3). High-risk regulons (RC3 and RC4) showed greater activity in the C1A phenotype. They were suppressed in the C1B phenotype, contrary to what was observed for the low-risk regulons (RC1 and RC2), which agrees with the pattern observed in the TCGA-ACC cohort (Figure 5A,B).
A
TCGA-ACC cohort
B
ENSAT cohort
C1A
C1B
C1A
C1B
Single sample regulon activity (median of the cluster)
Kruskal-Wallis, p < 2.2 x 10-16
Kruskal-Wallis, p < 2.2 x 10-16
Single sample regulon activity (median of the cluster)
Kruskal-Wallis, p = 2.5x 10-10
Kruskal-Wallis, p = 4 x 10-11
2
*
ns
ns
ns
2
ns
ns
ns
ns
1
1
0
0
O
1
-1
=
RC1 RC2 RC3 RC4
RC1
RC2
RC3
RC4
RC1
RC2
RC3
RC4
RC1
RC2
RC3
RC4
(n=43 samples in each cluster) (n=35 samples in each cluster)
(n=26 samples in each cluster)
(n=18 samples in each cluster)
C
TCGA-ACC cohort
D
ENSAT cohort
Tumor Stages 1 and 2
Tumor Stages 3 and 4
Tumor Stages 1 and 2
Tumor Stages 3 and 4
Single sample regulon activity (median of the cluster)
Kruskal-Wallis, p = 5.9 x 10-5
Kruskal-Wallis, p = 1.2 x 10-6
Single sample regulon activity (median of the cluster)
Kruskal-Wallis, p = 4.7 x 10-3
Kruskal-Wallis, p = 1 x 104
**
*
*
*
*
**
*
**
2
ns
**
ns
ns
**
ns
2
ns
*
ns
**
ns
**
ns
1
1
O
.
.
O
O
0
o
8
0
SE
Q
.
O
&
.
-1
.
.
-1
I
O P
(n=46 samples in each cluster) (n=30 samples in each cluster)
RC1
RC2
RC3
RC4
RC1
RC2
RC3
RC4
RC1 RC2 RC3 RC4 (n=28 samples in each cluster)
RC1
RC2
RC3
RC4
(n=15 samples in each cluster)
Regulon Clustes
Vital Status
☐ RC1 ☐ RC3
Alive
☐ RC2
RC4
Dead
Concerning the tumor stages in the TCGA-ACC and ENSAT cohorts, RC1 and RC2 had stronger activity at low tumor stages, while RC3 and RC4 had stronger activity at higher stages (Figure 5C,D). We summarized the results for the 369 regulons in Supplementary Table S4.
3.4. NR5A1 Relates to Worse Outcomes in TCGA-ACC and ENSAT Cohorts. CENPA Has a Strong Association with Proliferation and Relates to a Bad Prognosis
NR5A1 is a well-described TF related to a worse prognosis, adrenal differentiation, and steroidogenesis in ACC. NR5A1 presented 248 targets in the regulatory network, with 176 negatives and 72 positives (Supplementary Figure S6, Supplementary Tables S1 and S2). In the Cox multivariate analysis, the NR5A1 regulon related to worse outcomes in PFI (HR = 2.15 [95% CI, 1.49-3.11], p-value =7.12 × 10-5) and OS (HR = 1.94 [95% CI, 1.21-3.11], p-value = 7.13 x 10-3) in the TCGA-ACC cohort, and to OS in the ENSAT cohort (HR = 3.86 [95% CI, 1.61-9.25], p-value = 7.45 × 10-3) (Supplementary Table S3). The Kaplan-Meier analysis also presented a significant value relating to worse outcomes in OS both in the TCGA-ACC and the ENSAT cohorts (Figure 6A and 6B, respectively).
A
TCGA-ACC Hazard Ratio (95% CI)
B
ENSAT Hazard Ratio (95% CI)
0.1
1.0
15
0.1
1.0
15
Regulons and other covariates
Regulons and other covariates
Age
Age
Tumor.Stage
Tumor.Stage
NR5A1
NR5A1
CENPA
CENPA
associated, HR<1 not associated
associated, HR>1
associated, HR<1 not associated
associated, HR>1
C
other covariates
TCGA-ACC Samples (n=78)
Logrank P: 2.69 x 10-3
D
other covariates
78
NR5A1
1.0-
44
NR5A1
1.0
Logrank P: 1.83 x 10-4
Survival probability
0.8
ENSAT Samples (n=44)
61
0.8
31
Survival probability
0.6
41
0.6
0.4-
21
0.4
21
0.2
11
Positive dES: 38 (20)
0.2
undetermined: 4 (2)
Positive dES: 20 (14)
1
0.0-
Negative dES: 36 (5)
1
undetermined: 8 (4)
0.0
Negative dES: 16 (1)
-2
-1
0
1
₪
Regulon activity (dES)
HSP
LSP
Proliferation-
C1A-
C1B-
Stage1-
Stage2
Stage3-
Stage4
0
50
100
150
200
-2
-1
0
1
C1A
C1B-
Stage1-
Stage2-
Stage3-
Stage4-
0
50
100
150
200
Months
E
Regulon activity (dES)
Months
F
TCGA-ACC Samples (n=78)
78
CENPA
1.0
Logrank P: 1.44 x 10-10
44
CENPA
1.0
Logrank P: 3.33 x 10-5
61
Survival probability
0.8
ENSAT Samples (n=44)
0.8-
31
Survival probability
0.6
41
0.6
0.4-
21
0.4
21
0.2
11
Positive dES: 27 (20)
0.2
undetermined: 18 (3)
Positive dES:
16 (12)
1
0.0-
Negative dES: 33 (4)
undetermined. 9 (5)
1
0.0-
Negative dES: 19 (2)
-2
-1
0
1
2
HSP
LSP
Proliferation
C1A
C1B
Stage1
Stage2-
Stage3
Stage4
0
50
100
150
200
Regulon activity (dES)
-2
-1
0
1
N
C1A,
C1B.
Stage1
Stage2.
Stage3.
Stage4-
0
50
100
150
200
Months
Regulon activity
(dES)
Months
This regulon presented a significant relation both with HSP (dES = 1.49, p-value = 1.65 x 10-3) and Proliferation IS phenotypes (dES = 1.33, p-value = 1.36 x 10-3) and clustered in the RC3. Concerning the immune features analyzed, its activity showed a significant negative correla- tion with leukocyte fraction (p = - 0.52, p-value = 4.84 x 10-6), TGF-ß response (p = - 0.41, p-value = 1.71 × 10-2, the lowest correlation among the 369 regulons), Macrophage Regulation (p = - 0.52, p-value = 6.84 x 10-6), Lymphocyte Infiltration Signature Score (p =- 0.56, p-value = 4.90 x 10-7), TCR Shannon (p = - 0.61, p-value = 2.14 × 10-3), and TCR Richness (p = - 0.47, p-value = 2.87 x 10-4). Moreover, this regulons’ activity showed a significant positive correlation with the Proliferation and Wound Healing sig- natures (p = 0.32, p-value = 5.84 x 10-3, and p = 0.35, p-value = 2.07 x 10-3, respectively) (Supplementary Table S4).
Within the 369 regulons, we found 150 pairs of regulons that significantly shared targets (Supplementary Table S8) and were therefore called “Duals” [50]. Despite being the second regulon with the most targets (n = 248, Supplementary Figure S6), the NR5A1 regulon did not show significant target sharing with other prognostic regulons.
Of the inferred Duals, Centromeric Protein A (CENPA), Small Nuclear RNA Activating Complex Polypeptide 4 (SNAPC4), and LIN54 [a component of the LINC (Linker of Nucle- oskeleton and Cytoskeleton) complex, an essential regulator of cell cycle genes], were the regulons with the most associations (9, 8, and 7, respectively) (Supplementary Table S8). Interestingly, CENPA had the highest score in the correlation with Proliferation (p = 0.90, p-value = 4.01 × 10-27) and Wound Healing signatures (p =0.79, p-value = 2.55 × 10-15), and a significant negative correlation with leukocyte fraction (p = - 0.55, p-value = 1.44 × 10-6), Macrophage Regulation (p = - 0.42, p-value = 4.29 x 10-4), and Lymphocyte Infiltration Signature Score (p = - 0.48, p-value = 2.67 × 10-5) (Supplementary Table S4).
We found that CENPA was the main regulator of the MKI-67 gene (MI = 0.77), the most common histopathologic marker of proliferation, followed by FOXM1 (MI = 0.59), MXD3 (MI = 0.48), and DNMT1 (MI = 0.34) (Supplementary Table S1). MXD3 and DNMT1 are CENPA Duals (Supplementary Table S8), and FOXM1, MXD3, and DNMT1 are CENPA- positive targets.
CENPA was significantly associated with both the Steroid IP and the Proliferation IS (p = 1.45, p-value = 1.65 x 10-3, and p = 1.32, p-value = 1.36 x 10-3, respectively) and clustered together with RC3 (Supplementary Table S4). In the Cox’s analysis, CENPA related to worse outcomes in the TCGA-ACC cohort in PFI (HR = 2.23 [95% CI, 1.54-3.21], p-value = 4.55 × 10-5) and in OS (HR = 3.24 [95% CI, 1.88-5.57], p-value = 1.82 × 10-4), in addition to the ENSAT cohort for OS (HR = 2.52 [95% CI, 1.39-4.59], p-value = 7.45 x 10-3) (Supplementary Table S3). In the Kaplan-Meier analysis, CENPA presented a significant relation in the TCGA-ACC and ENSAT cohorts (Figure 6A and B, respectively). Figure 7A presents CENPA’s targets, and Figure 7B shows their distribution in a karyogram.
In the TCGA-ACC cohort, the three regulons with higher negative Spearman cor- relation with CENPA’s activity were THRB (Thyroid Hormone Receptor Beta), STAT5B (Signal Transducer and Activator of Transcription 5B), and ZBTB4 (Zinc Finger And BTB Do- main Containing 4) (p =- 0.87, -0.86, and -0.85, respectively, all p-values < 2.2 × 10-16). All of them from RC1 related to good prognosis in the TCGA-ACC and ENSAT cohorts (Supplementary Tables S3 and S4) and presented a negative correlation with the Prolifera- tion signature (p = - 0.72, -0.75 and -0.68, respectively, all p-values < 10-9).
A
E2F2
ERCC6L
RRM1
GMNN
B CENPA targets Karyogram
COQ8A
SKA2
PBK
SGO1
KIF23
ECHDC3
NEK2
LIG1
chr1
TP73
EIF4EBP3
PTTG1
SEPT14P12
RRM2
MFSD6
VEPH1
TICRR
IQGAP3
MAOB
AC091057.1
chr2
NCAPG
BUB1
chr3
PURB
NOAPH
MYRIP
KPNĄ2
SCN11ASUV39H2
UBERS
MKI67
DLGAP5
BCL2L12
chr4
KAT14
SGO2
NDRG4
STK33
chr5
LMNB1
KIF11
AURKB
PSRC1
TMEM144
AC026401.STIL
FOXM
RHEBL1
CDCA5
chr6
SLC30A2
PIF1
chr7
JBE2T ARHGEF28 LRR1
FAM72C
CAB39L
BIRC5
HISTIHIE
ZBTB4
chr8
ZWINT
KNL1
NTF3
STATSA
KIF2C
SHISA4
chr9
DEPDC1 MCM7
KIF18B
chr10
CDC45
FEN1
ARSD
CDS2
C2orf40 CDC20
TEDC1
CENPA
TALT
CST3
C1QTNF2
chr11
CDCA2
RAB37
TACC3
chr12
CDC7 CCNB2 RAD51
CEP55
MCM10
CDK1
chr13
DEPDC1B
CCNA2
NKAIN2
chr14
PCMTD1 CDC25C
SPAG5
TOP2A
RM12
UBE2C
SPC25
chr15
PCLAF
CDCA7
LINC00634
TUBB6
chr16
TYMS
RFC2
CKAP2L
SYT14
SLC24A1
NUF2
PHGDH
CHMP4CMCM3
CENPL
chr17
EXO1
CDCA3FBX02
PLKY
CENPM
CFAP70
STAT5B
PLAC1
chr18
EZH2
DNMT1
TAF1B
chr19
PRC1
RACGAP1
NDC80
MELK
CDKN3
CDCA8
CENPF
KIFCY
chr20
GPR3
HAUS8
FBXO25
FPGT-
MXD3
ASPM
TNNI3K
KIF4A
APOBEC3B
0
chr22
BUB1B
EZH1
KNSTRN
ANK1
DNA2
chrX
LMOD1AURKA
VRK1
TEX30
KIF14
KIF20A
TROAP
0 Mb
50 Mb
100 Mb
150 Mb
200 Mb
250 Mb
Legend
NIPAL2
DMTN
MYLK2
CCNB1
CPQ
PPM1L
UHRF1
STMN1
Regulator
PEX12
SKA3
Positive Targets
HJURP
GTSE1
ANP32B
Association Direction
Negative Targets
GSAP
Negative (n=44)
Positive (n=120)
4. Discussion
The TCGA-ACC cohort is the largest ACC cohort, thus being the most suitable for regulatory network inference. The inferred network had a good balance between positive and negative targets (Supplementary Figure S2), with 611 regulons comprising more than 15 positive and 15 negative targets. We described regulons of prognostic value and provided an overview of the main regulons that control the expression of target genes and, consequently, the observed phenotypes associated with OS. To achieve this, we made functional annotations associating the regulon activity with molecular phenotypes and immunological characteristics.
Here we present 369 regulons associated with OS, where 188 are related to a good prognosis and 181 to a worse ACC prognosis. Within these regulons, we identified sub- groups with specific relationships for each set of molecular characteristics and functional annotations, resulting in four RCs. Regulon cluster 1 and RC3 are more related to the Steroid Phenotype, while RC2 and RC4 are more related to the molecular Proliferation Score as defined by Zheng et al. (2016) [7]. As the high Score Proliferation classification was present in 22 of the 47 HSP cases and only 4 of the 31 LSP cases, we distinguished these two features, the Steroid and Proliferation classifications, by analyzing them independently, removing the interference from each other.
RC1 and RC3 are related to Steroid IP, with RC1 associated with LSP and RC3 with HSP. These two clusters presented opposite profiles, both in activity and functional annotations. These clusters were at high activity in the immune and proliferation pathways, as previously shown for the Steroid Phenotype [6]. The strong correlation between RC1 activity and the pro-immune features (Figure 4, Supplementary Figures S3 and S4), in addition to increased activity in early tumor stages (Figure 5B), suggests that regulons in the RC1 may play a role
in activating the immune system against different targets at early phases of ACC, which is progressively lost by opposing forces of the RC3 regulons in more advanced stages.
We also observed that RC2 and RC4 follow opposite patterns in terms of activity and characteristics. RC2, with a good prognosis, is related to low scores in the Prolif- eration IS, while RC4, with a poor prognosis, is related to high scores in this pheno- type. We identified from the Hallmarks enrichment (Figure 4, Supplementary Figure S3, Supplementary Table S7) that these prognostic clusters are particularly associated with the PI3K Akt mTOR and Hedgehog signaling pathways, which may help uncover risk factors other than immune response and Steroid Phenotype [5,6,60].
Our functional analysis using the MSigDb Hallmarks presented an overall panorama, providing insights into the biology and function of each group of inferred regulons. How- ever, the MSigDb Hallmarks represent themes of gene sets rather than pathways [30], which may result in supposed contradictory findings. For example, the classical P53 pathway activates DNA repair in response to DNA damage [61,62]. However, in our analysis, the P53 pathway and the DNA repair Hallmarks presented opposite activities (i.e., low for P53 and high for DNA repair in both RC3 and RC4). Increased activity of DNA repair is observed after increased proliferation [62,63], which may explain the observed profiles in RC3 and RC4. In addition, some mutations, such as in the TP53 gene, may alter the P53 pathway and its DNA repair activity [61,64]. The overall pattern of immune response and proliferation pathways was coherent with the good and poor prognosis clusters. Specific functional correlations need further clarification using in vitro and in vivo studies.
The TCGA pan-cancer analysis, including the ACC cohort by Thorsson et al. (2018) [14], shows that TFs regulating immune modulators tend to be shared between different tissue- of-origin malignancies, in contrast to somatic mutations. They also highlight some immune- related TFs shared among the tumors, of which three appeared in our analysis with association to good prognosis in ACC: FLI1 (Friend Leukemia Virus Integration 1), STAT5A (Signal Transducer and Activator of Transcription 5A), and IKZF1 (Ikaros Family Zinc Finger Protein 1). For example, in our results, IKZF1 appeared related to good prognosis (HR = 0.61 [95% CI, 0.40-0.94], p-value = 2.65 × 10-2), and with the strongest association with the LSP (Steroid IP score of -1.76, p-value = 1.65 x 10-3) and the highest positive cor- relation with leukocyte fraction (p = 0.83, p-value = 5.47 x 10-18), Macrophage Regulation (p = 0.86, p-value = 1.12 × 10-21), and Lymphocyte Infiltration Signature Score (p = 0.83, p-value = 1.56 x 10-18) among the 369 regulons (Supplementary Table S4). Finding TFs shared between other tumors and ACC with importance to the immune activation may help to optimize clinical and research efforts in this rare carcinoma.
In ACC, the NR5A1 TF, an important player in adrenal development and ACC tu- morigenesis [65,66], generated one of the largest regulons in our analysis with 248 targets. Consistent with studies on the overexpression of this TF relating to a worse outcome and increased steroid metabolism [23,24,67], the NR5A1 regulon correlated with low OS in the multivariate Cox and Kaplan-Meier analyses in the TCGA-ACC cohort, which was also confirmed in the ENSAT cohort (Figure 6). Interestingly, of the 369 prognostic regulons, the NR5A1 regulon showed the lowest correlation with the TGF-ß response. Although steroid hormones are well-known immunosuppressors [68,69], the pathways by which they dampen the immune response in the ACC microenvironment are not completely understood. Further studies are needed to examine the immune variables directly altered by NR5A1 overexpression in the context of ACC.
Remarkably, we identified CENPA with interesting correlations with many prolifer- ative markers. CENPA is a histone H3-like protein involved in centromeric nucleosome formation [70]. In this study, the CENPA regulon showed the highest association with the Proliferation signature. In addition, CENPA also appears as the main regulator of MKI67 (Ki-67), a common prognostic and proliferation marker widely used in cancer histopathology [22].
CENPA was described as relating to proliferation and prognosis in ovarian cancer [71], being crucial in prostate cancer growth [72], regulating metabolic reprogramming in colon
cancer cells leading to its growth [73], and being associated with immune infiltration and prognosis in lung cancer [74]. Notably, CENPA overexpression was also associated with proliferation and metastasis in kidney carcinoma by activating the WNT/ß-Catenin signal- ing pathway [75], an important pathway described in the progression of ACC [3]. In our analysis, CENPA appeared positively enriched in this pathway (Supplementary Table S7). Further studies may elucidate whether the relationship in renal carcinoma between CENPA activity and the WNT/ß-Catenin signaling pathway may apply to ACC.
The three regulons with the highest negative correlation with CENPA’s activity were THRB, STAT5B, and ZBTB4, all from RC1 and of good prognosis. Interestingly, ZBTB4 is known to act as a tumor suppressor in various types of cancers (prostate cancer [76]; glioma [77]; colorectal cancer [78]; breast cancer [79]; Ewing sarcoma [80]), but more impor- tantly, it is being investigated as a promissory anti-cancer druggable target. For example, Kim et al. (2012) [76] showed that Methyl 2-cyano-3,11-dioxo-18ß-olean-1,12-dien-30-oate (CDODA-Me) could increase ZBTB4 expression in vivo and in vitro, having an antitumori- genic activity in prostate cancer. Whether the increased expression of ZBTB4 could inhibit CENPA and its proliferative activity may be an interesting area for future studies.
It should be noted that our study had some limitations. Our selection strategy may have removed some interesting regulons, which are not listed in the 369 prognostic regulons, thereby excluding their possible role in the ACC regulatory scenario. For example, FOXM1 is assigned as a prognostic marker in ACC by Yuan et al. (2018) [81]. This TF presented 111 targets, of which only 11 were downregulated, thus not passing through our filter of at least 15 positive and 15 negative targets. We also used protein-coding genes, excluding possible important regulators, such as miRNAs or lncRNAs. Despite these limitations, our results may offer a reference for future studies aiming to understand transcriptional alterations in ACC, prognostic markers, or therapeutic targets.
5. Conclusions
In conclusion, we have generated a regulatory network for ACC, evaluated the reg- ulons inferred in concern to OS, clustered them in four RCs, and investigated how they relate to characteristics associated with worse outcomes like the steroid phenotype and the immune and proliferation pathways. The list of prognostic regulons, and their charac- terization, may open new research avenues and relevant questions to be answered in this hard-to-treat and aggressive malignancy.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14215279/s1, reference [14,30,31]. Figure S1: Distribution of the positive and negative targets of the regulons; Figure S2: Consensus matrix and silhouettes of regulon clusters; Figure S3: Boxplots comparing the Hallmark enrichment scores in the regulon clusters; Figure S4: Boxplots comparing the Spearman’s correlation between the immune features described by Thorsson et al. (2018) [14] and the regulon activity in the clusters; Figure S5: A heatmap showing regulon activity in the ENSAT cohort (n = 44 ACC samples); Figure S6: Targets of the NR5A1 regulon and the related karyogram; Table S1: ACC regulatory network presenting the mutual information between the 1605 transcription factors and their targets; Table S2: Number of positive and negative targets for each regulon; Table S3: Results for the multivariate Cox analysis for overall survival in the TCGA-ACC and ENSAT cohorts, and PFI in the TCGA-ACC cohort; Table S4: A master table resuming the data generated for the 369 prognostic regulons; Table S5: Activity of the 369 prognostic regulons in each sample for the TCGA-ACC and ENSAT cohorts; Table S6: Master Regulator Analysis (MRA) and two-tails gene set enrichment analysis (GSEA2) for the enrichment between the 369 regulons and the Steroid and Proliferation Phenotypes; Table S7: Regulon activity in MsigDb Hallmarks; Table S8: Duals correlation.
Author Contributions: Conceptualization, J.C.D.M. and M.A.A.C .; methodology, J.C.D.M. and M.A.A.C .; software, M.A.A.C .; formal analysis, J.C.D.M., M.A.A.C., J.F.d.M. and B.C.F .; data curation, J.C.D.M .; writing-original draft preparation, J.C.D.M. and B.C.F .; writing-review and editing, J.C.D.M., J.M.M., J.S.S., L.M.A., J.F.d.M., B.C.F. and M.A.A.C .; visualization, J.C.D.M., J.M.M., J.S.S. and M.A.A.C .; supervision, L.M.A., J.F.d.M., B.C.F. and M.A.A.C .; project administration, M.A.A.C., L.M.A., J.F.d.M. and B.C.F .; funding acquisition, B.C.F. All authors have read and agreed to the published version of the manuscript.
Funding: J.C.D.M., J.M.M. and J.S.S.R. received scholarships from the BIG DATA innovation program from the Associação Hospitalar de Proteção à Infância Raul Carneiro-AHPIRAC (2020). M.A.A.C., L.M.A. and J.F.d.M. are funded by the Conselho Nacional de Desenvolvimento Científico e Tec- nológico (CNPq). M.A.A.C. is also funded by NAPI Bioinformática-Fundação Araucária (FA). J.C.D.M., J.M.M. and J.S.S. were funded by Associação Hospitalar de Proteção à Infância Dr. Raul Carneiro (AHPIRAC).
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data and software used in this study are publicly available in the sources described in the Section 2. All scripts and R data (rtni object) to generate results, figures and tables for this study are available on the GitHub repository (https://github.com/sysbiolab/Sup_ Material_Muzzi2022) (accessed on 30 August 2022).
Acknowledgments: Data used in this manuscript were obtained from The Cancer Genome Atlas (TCGA) ACC database and of the European Network for the Study of Adrenal Tumors (ENSAT).
Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
References
1. Wasserman, J.D .; Zambetti, G.P .; Malkin, D. Towards an understanding of the role of p53 in adrenocortical carcinogenesis. Mol. Cell Endocrinol. 2012, 351, 101-110. [CrossRef] [PubMed]
2. Else, T .; Kim, A.C .; Sabolch, A .; Raymond, V.M .; Kandathil, A .; Caoili, E.M .; Jolly, S .; Miller, B.S .; Giordano, T.J .; Hammer, G.D. Adrenocortical carcinoma. Endocr. Rev. 2013, 35, 282-326. [CrossRef] [PubMed]
3. Cosentini, D .; Grisanti, S .; Volta, A.D .; Laganà, M .; Fiorentini, C .; Perotti, P .; Sigala, S .; Berruti, A. Immunotherapy failure in adrenocortical cancer: Where next? Endocr. Connect. 2018, 7, E5-E8. [CrossRef] [PubMed]
4. Baechle, J.J .; Hanna, D.N .; Sekhar, K.R .; Rathmell, J.C .; Rathmell, W.K .; Baregamian, N. Integrative computational immunogenomic profiling of cortisol-secreting adrenocortical carcinoma. J. Cell Mol. Med. 2021, 25, 10061-10072. [CrossRef]
5. Landwehr, L .- S .; Altieri, B .; Schreiner, J .; Sbiera, I .; Weigand, I .; Kroiss, M .; Fassnacht, M .; Sbiera, S. Interplay between glucocorti- coids and tumor-infiltrating lymphocytes on the prognosis of adrenocortical carcinoma. J. Immunother. Cancer 2020, 8, e000469. [CrossRef]
6. Muzzi, J.C.D .; Magno, J.M .; Cardoso, M.A .; de Moura, J .; Castro, M.A.A .; Figueiredo, B.C. Adrenocortical carcinoma steroid profiles: In silico pan-cancer analysis of TCGA data uncovers immunotherapy targets for potential improved outcomes. Front. Endocrinol. 2021, 12, 672319. [CrossRef]
7. Zheng, S .; Cherniack, A.D .; Dewal, N .; Moffitt, R.A .; Danilova, L .; Murray, B.A .; Lerario, A.M .; Else, T .; Knijnenburg, T.A .; Ciriello, G .; et al. Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer Cell 2016, 29, 723-736. [CrossRef]
8. Wirapati, P .; Sotiriou, C .; Kunkel, S .; Farmer, P .; Pradervand, S .; Haibe-Kains, B .; Desmedt, C .; Ignatiadis, M .; Sengstag, T .; Schütz, F .; et al. Meta-analysis of gene expression profiles in breast cancer: Toward a unified understanding of breast cancer subtyping and prognosis signatures. Breast Cancer Res. 2008, 10, R65. [CrossRef]
9. de Reyniès, A .; Assié, G .; Rickman, D.S .; Tissier, F .; Groussin, L .; René-Corail, F .; Dousset, B .; Bertagna, X .; Clauser, E .; Bertherat, J. Gene expression profiling reveals a new classification of adrenocortical tumors and identifies molecular predictors of malignancy and survival. J. Clin. Oncol. 2009, 27, 1108-1115. [CrossRef]
10. Galon, J .; Costes, A .; Sanchez-Cabo, F .; Kirilovsky, A .; Mlecnik, B .; Lagorce-Pagès, C .; Tosolini, M .; Camus, M .; Berger, A .; Wind, P .; et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science 2006, 313, 1960-1964. [CrossRef]
11. Gajewski, T.F .; Corrales, L .; Williams, J .; Horton, B .; Sivan, A .; Spranger, S. Cancer immunotherapy targets based on understanding the T cell-inflamed versus non-T cell-inflamed tumor microenvironment. Adv. Exp. Med. Biol. 2017, 1036, 19-31. [CrossRef] [PubMed]
12. Galon, J .; Bruni, D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat. Rev. Drug Discov. 2019, 18, 197-218. [CrossRef] [PubMed]
13. Hegde, P.S .; Karanikas, V .; Evers, S. The where, the when, and the how of immune monitoring for cancer immunotherapies in the era of checkpoint inhibition. Clin. Cancer Res. 2016, 22, 1865-1874. [CrossRef] [PubMed]
14. Thorsson, V .; Gibbs, D.L .; Brown, S.D .; Wolf, D .; Bortone, D.S .; Ou Yang, T .- H .; Porta-Pardo, E .; Gao, G.F .; Plaisier, C.L .; Eddy, J.A .; et al. The immune landscape of cancer. Immunity 2018, 48, 812-830.e14. [CrossRef]
15. Le Tourneau, C .; Hoimes, C .; Zarwan, C .; Wong, D.J .; Bauer, S .; Claus, R .; Wermke, M .; Hariharan, S .; Von Heydebreck, A .; Kasturi, V .; et al. Avelumab in patients with previously treated metastatic adrenocortical carcinoma: Phase 1b results from the JAVELIN solid tumor trial. J. Immunother. Cancer 2018, 6, 111. [CrossRef]
16. Carneiro, B.A .; Konda, B .; Costa, R.B .; Costa, R.L.B .; Sagar, V .; Gursel, D.B .; Kirschner, L.S .; Chae, Y.K .; Abdulkadir, S.A .; Rademaker, A .; et al. Nivolumab in metastatic adrenocortical carcinoma: Results of a phase 2 trial. J. Clin. Endocrinol. Metab. 2019, 104, 6193-6200. [CrossRef]
17. Habra, M.A .; Stephen, B .; Campbell, M .; Hess, K .; Tapia, C .; Xu, M .; Rodon Ahnert, J .; Jimenez, C .; Lee, J.E .; Perrier, N.D .; et al. Phase II clinical trial of pembrolizumab efficacy and safety in advanced adrenocortical carcinoma. J. Immunother. Cancer 2019, 7, 253. [CrossRef]
18. Raj, N .; Zheng, Y .; Kelly, V .; Katz, S.S .; Chou, J .; Do, R.K.G .; Capanu, M .; Zamarin, D .; Saltz, L.B .; Ariyan, C.E .; et al. PD-1 blockade in advanced adrenocortical carcinoma. J. Clin. Oncol. 2020, 38, 71-80. [CrossRef]
19. McGregor, B.A .; Sonpavde, G.P. Rare Genitourinary Malignancies: Current Status and Future Directions of Immunotherapy. Eur. Urol. Focus 2020, 6, 14-16. [CrossRef]
20. Mete, O .; Erickson, L.A .; Juhlin, C.C .; de Krijger, R.R .; Sasano, H .; Volante, M .; Papotti, M.G. Overview of the 2022 WHO classification of adrenal cortical tumors. Endocr. Pathol. 2022, 33, 155-196. [CrossRef]
21. Hodgson, A .; Pakbaz, S .; Mete, O. A diagnostic approach to adrenocortical tumors. Surg. Pathol. Clin. 2019, 12, 967-995. [CrossRef] [PubMed]
22. Fassnacht, M .; Dekkers, O.M .; Else, T .; Baudin, E .; Berruti, A .; de Krijger, R.R .; Haak, H.R .; Mihai, R .; Assie, G .; Terzolo, M. European society of endocrinology clinical practice guidelines on the management of adrenocortical carcinoma in adults, in collaboration with the European network for the study of adrenal tumors. Eur. J. Endocrinol. 2018, 179, G1-G46. [CrossRef] [PubMed]
23. Mete, O .; Asa, S.L .; Giordano, T.J .; Papotti, M .; Sasano, H .; Volante, M. Immunohistochemical biomarkers of adrenal cortical neoplasms. Endocr. Pathol. 2018, 29, 137-149. [CrossRef] [PubMed]
24. Lalli, E. Adrenocortical development and cancer: Focus on SF-1. J. Mol. Endocrinol. 2010, 44, 301-307. [CrossRef]
25. Crona, J .; Beuschlein, F. Adrenocortical carcinoma-Towards genomics guided clinical care. Nat. Rev. Endocrinol. 2019, 15, 548-560. [CrossRef]
26. Castro, M.A.A.A .; de Santiago, I .; Campbell, T.M .; Vaughn, C .; Hickey, T.E .; Ross, E .; Tilley, W.D .; Markowetz, F .; Ponder, B.A.J.J .; Meyer, K.B. Regulators of genetic risk of breast cancer identified by integrative network analysis. Nat. Genet. 2016, 48, 12-21. [CrossRef]
27. Robertson, A.G .; Kim, J .; Al-Ahmadie, H .; Bellmunt, J .; Guo, G .; Cherniack, A.D .; Hinoue, T .; Laird, P.W .; Hoadley, K.A .; Akbani, R .; et al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell 2017, 171, 540-556.e25. [CrossRef]
28. Kamoun, A .; de Reyniès, A .; Allory, Y .; Sjödahl, G .; Robertson, A.G .; Seiler, R .; Hoadley, K.A .; Groeneveld, C.S .; Al-Ahmadie, H .; Choi, W .; et al. A Consensus molecular classification of muscle-invasive bladder cancer. Eur. Urol. 2020, 77, 420-433. [CrossRef]
29. Cangiano, M .; Grudniewska, M .; Salji, M.J .; Nykter, M .; Jenster, G .; Urbanucci, A .; Granchi, Z .; Janssen, B .; Hamilton, G .; Leung, H.Y .; et al. Gene regulation network analysis on human prostate orthografts highlights a potential role for the JMJD6 regulon in clinical prostate cancer. Cancers 2021, 13, 2094. [CrossRef]
30. Liberzon, A .; Birger, C .; Thorvaldsdóttir, H .; Ghandi, M .; Mesirov, J.P .; Tamayo, P. The molecular signatures database hallmark gene set collection. Cell Syst. 2015, 1, 417-425. [CrossRef]
1. Assié, G .; Letouzé, E .; Fassnacht, M .; Jouinot, A .; Luscap, W .; Barreau, O .; Omeiri, H .; Rodriguez, S .; Perlemoine, K .; René-Corail, F .; et al. Integrated genomic characterization of adrenocortical carcinoma. Nat. Genet. 2014, 46, 607-612. [CrossRef]
32. Colaprico, A .; Silva, T.C .; Olsen, C .; Garofano, L .; Cava, C .; Garolini, D .; Sabedot, T.S .; Malta, T.M .; Pagnotta, S.M .; Castiglioni, I .; et al. TCGAbiolinks: An R/bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016, 44, e71. [CrossRef] [PubMed]
33. Silva, T.C .; Colaprico, A .; Olsen, C .; D’Angelo, F .; Bontempi, G .; Ceccarelli, M .; Noushmehr, H. TCGA workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages. F1000Research 2016, 5, 1542. [CrossRef] [PubMed]
34. Mounir, M .; Lucchetta, M .; Silva, T.C .; Olsen, C .; Bontempi, G .; Chen, X .; Noushmehr, H .; Colaprico, A .; Papaleo, E. New functionalities in the TCGA biolinks package for the study and integration of cancer data from GDC and GTEx. PLOS Comput. Biol. 2019, 15, e1006701. [CrossRef] [PubMed]
35. Liu, J .; Lichtenberg, T .; Hoadley, K.A .; Poisson, L.M .; Lazar, A.J .; Cherniack, A.D .; Kovatich, A.J .; Benz, C.C .; Levine, D.A .; Lee, A.V .; et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell 2018, 173, 400-416.e11. [CrossRef] [PubMed]
36. Morgan, M .; Shepherd, L. Annotation Hub: Client to Access AnnotationHub Resources, version 2.22.0; R Package Bioconductor Package Maintainer, 2020. [CrossRef]
37. Love, M.I .; Huber, W .; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [CrossRef] [PubMed]
38. Fletcher, M.N.C.C .; Castro, M.A .; Wang, X .; de Santiago, I .; O’Reilly, M .; Chin, S .- F .; Rueda, O.M .; Caldas, C .; Ponder, B.A.J.J .; Markowetz, F .; et al. Master regulators of FGFR2 signalling and breast cancer risk. Nat. Commun. 2013, 4, 2464. [CrossRef]
39. Lambert, S.A .; Jolma, A .; Campitelli, L.F .; Das, P.K .; Yin, Y .; Albu, M .; Chen, X .; Taipale, J .; Hughes, T.R .; Weirauch, M.T. The human transcription factors. Cell 2018, 172, 650-665. [CrossRef]
40. Margolin, A.A .; Wang, K .; Lim, W.K .; Kustagi, M .; Nemenman, I .; Califano, A. Reverse engineering cellular networks. Nat. Protoc. 2006, 1, 662-671. [CrossRef]
41. Cohen, J. Statistical Power Analysis for the Behavioral Sciences; Taylor and Francis: Hoboken, NJ, USA, 1988; ISBN 9781134742707.
42. Benjamini, Y .; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 1995, 57, 289-300. [CrossRef]
43. Cox, D.R. Regression models and life-tables. J. R. Stat. Soc. Ser. B Stat. Methodol. 1972, 34, 187-202. [CrossRef]
44. Groeneveld, C.S .; Chagas, V.S .; Jones, S.J.M .; Robertson, A.G .; Ponder, B.A.J .; Meyer, K.B .; Castro, M.A.A. RTN survival: An R/bioconductor package for regulatory network survival analysis. Bioinformatics 2019, 35, 4488-4489. [CrossRef] [PubMed]
45. Kaplan, E.L .; Meier, P. Nonparametric estimation from incomplete observations. J. Am. Stat. Assoc. 1958, 53, 457-481. [CrossRef]
46. Mantel, N. Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemother Rep. 1966, 50, 163-170. [PubMed]
47. Peto, R .; Peto, J. Asymptotically efficient rank invariant test procedures. J. R. Stat. Soc. Ser. A Gen. 1972, 135, 185. [CrossRef]
48. Wilkerson, M.D .; Hayes, D.N. ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics 2010, 26, 1572-1573. [CrossRef]
49. Dolgalev, I. Msigdbr: MSigDB Gene Sets for Multiple Organisms in a Tidy Data Format. Available online: https://igordot.github. io/msigdbr/ (accessed on 30 June 2022).
50. Chagas, V.S .; Groeneveld, C.S .; Oliveira, K.G .; Trefflich, S .; De Almeida, R.C .; Ponder, B.A.J .; Meyer, K.B .; Jones, S.J.M .; Robertson, A.G .; Castro, M.A.A. RTNduals: An R/bioconductor package for analysis of co-regulation and inference of dual regulons. Bioinformatics 2019, 35, 5357-5358. [CrossRef]
51. Davis, S .; Meltzer, P.S. GEOquery: A bridge between the gene expression omnibus (GEO) and bioconductor. Bioinformatics 2007, 23, 1846-1847. [CrossRef]
52. R Core Team. R: A Language and Environment for Statistical Computing; R Core Team: Vienna, Austria.
53. Huber, W .; Carey, V.J .; Gentleman, R .; Anders, S .; Carlson, M .; Carvalho, B.S .; Bravo, H.C .; Davis, S .; Gatto, L .; Girke, T .; et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods 2015, 12, 115-121. [CrossRef]
54. Wilcoxon, F. Individual comparisons by ranking methods. Biom. Bull. 1945, 1, 80. [CrossRef]
55. Mann, H.B .; Whitney, D.R. On a test of whether one of two random variables is stochastically larger than the other. Ann. Math. Stat. 1947, 18, 50-60. [CrossRef]
56. Kruskal, W.H .; Wallis, W.A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 1952, 47, 583-621. [CrossRef]
57. Dunn, O.J. Multiple comparisons using rank sums. Technometrics 1964, 6, 241-252. [CrossRef]
58. Gu, Z .; Eils, R .; Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016, 32, 2847-2849. [CrossRef]
59. Castro, M.A .; Wang, X .; Fletcher, M.N .; Meyer, K.B .; Markowetz, F. RedeR: R/bioconductor package for representing modular structures, nested networks and multiple levels of hierarchical associations. Genome Biol. 2012, 13, R29. [CrossRef]
60. Fiorentini, C .; Grisanti, S .; Cosentini, D .; Abate, A .; Rossini, E .; Berruti, A .; Sigala, S. Molecular drivers of potential immunotherapy failure in adrenocortical carcinoma. J. Oncol. 2019, 2019, 6072863. [CrossRef]
61. Sengupta, S .; Harris, C.C. P53: Traffic cop at the crossroads of DNA repair and recombination. Nat. Rev. Mol. Cell Biol. 2005, 6, 44-55. [CrossRef]
62. Jaber, S .; Toufektchan, E .; Lejour, V .; Bardot, B .; Toledo, F. P53 downregulates the Fanconi anaemia DNA repair pathway. Nat. Commun. 2016, 7, 11091. [CrossRef]
63. Kotoula, V .; Krikelis, D .; Karavasilis, V .; Koletsa, T .; Eleftheraki, A.G .; Televantou, D .; Christodoulou, C .; Dimoudis, S .; Korantzis, I .; Pectasides, D .; et al. Expression of DNA repair and replication genes in non-small cell lung cancer (NSCLC): A role for thymidylate synthetase (TYMS). BMC Cancer 2012, 12, 342. [CrossRef]
64. Ge, Z .; Leighton, J.S .; Wang, Y .; Peng, X .; Chen, Z .; Chen, H .; Sun, Y.Y .; Yao, F .; Li, J .; Zhang, H .; et al. Integrated genomic analysis of the ubiquitin pathway across cancer types. Cell Rep. 2018, 23, 213-226.e3. [CrossRef]
65. Pianovski, M.A.D .; Cavalli, L.R .; Figueiredo, B.C.L .; Santos, S.C .; Doghman, M .; Ribeiro, R.C .; Oliveira, A.G .; Michalkiewicz, E .; Rodrigues, G.A .; Zambetti, G .; et al. SF-1 overexpression in childhood adrenocortical tumours. Eur. J. Cancer 2006, 42, 1040-1043. [CrossRef]
66. Doghman, M .; Karpova, T .; Rodrigues, G.A .; Arhatte, M .; de Moura, J .; Cavalli, L.R .; Virolle, V .; Barbry, P .; Zambetti, G.P .; Figueiredo, B.C .; et al. Increased steroidogenic factor-1 dosage triggers adrenocortical cell proliferation and cancer. Mol. Endocrinol. 2007, 21, 2968-2987. [CrossRef] [PubMed]
67. Sbiera, S .; Schmull, S .; Assié, G .; Voelker, H .- U .; Kraus, L .; Beyer, M .; Ragazzon, B .; Beuschlein, F .; Willenberg, H.S .; Hahner, S .; et al. High diagnostic and prognostic value of steroidogenic factor-1 expression in adrenal tumors. J. Clin. Endocrinol. Metab. 2010, 95, E161-E171. [CrossRef] [PubMed]
68. Kronfol, Z .; Starkman, M .; Schteingart, D.E .; Singh, V .; Zhang, Q .; Hill, E. Immune regulation in cushing’s syndrome: Relationship to hypothalamic-pituitary-adrenal axis hormones. Psychoneuroendocrinology 1996, 21, 599-608. [CrossRef]
69. Jefferies, W.M. Cortisol and immunity. Med. Hypotheses 1991, 34, 198-208. [CrossRef]
70. Howman, E.V .; Fowler, K.J .; Newson, A.J .; Redward, S .; MacDonald, A.C .; Kalitsis, P .; Choo, K.H.A. Early disruption of centromeric chromatin organization in centromere protein A (Cenpa) null mice. Proc. Natl. Acad. Sci. USA 2000, 97, 1148-1153. [CrossRef]
71. Han, J .; Xie, R .; Yang, Y .; Chen, D .; Liu, L .; Wu, J .; Li, S. CENPA is one of the potential key genes associated with the proliferation and prognosis of ovarian cancer based on integrated bioinformatics analysis and regulated by MYBL2. Transl. Cancer Res. 2021, 10, 4076-4086. [CrossRef]
72. Saha, A.K .; Contreras-Galindo, R .; Niknafs, Y.S .; Iyer, M .; Qin, T .; Padmanabhan, K .; Siddiqui, J .; Palande, M .; Wang, C .; Qian, B .; et al. The role of the histone H3 variant CENPA in prostate cancer. J. Biol. Chem. 2020, 295, 8537-8549. [CrossRef] [PubMed]
73. Liang, Y .- C .; Su, Q .; Liu, Y .- J .; Xiao, H .; Yin, H .- Z. Centromere protein A (CENPA) regulates metabolic reprogramming in the colon cancer cells by transcriptionally activating karyopherin subunit alpha 2 (KPNA2). Am. J. Pathol. 2021, 191, 2117-2132. [CrossRef]
74. Zhou, H .; Bian, T .; Qian, L .; Zhao, C .; Zhang, W .; Zheng, M .; Zhou, H .; Liu, L .; Sun, H .; Li, X .; et al. Prognostic model of lung adenocarcinoma constructed by the CENPA complex genes is closely related to immune infiltration. Pathol. Res. Pract. 2021, 228, 153680. [CrossRef]
75. Wang, Q .; Xu, J .; Xiong, Z .; Xu, T .; Liu, J .; Liu, Y .; Chen, J .; Shi, J .; Shou, Y .; Yue, C .; et al. CENPA promotes clear cell renal cell carcinoma progression and metastasis via Wnt/ -catenin signaling pathway. J. Transl. Med. 2021, 19, 4-17. [CrossRef] [PubMed]
76. Kim, K .; Chadalapaka, G .; Pathi, S.S .; Jin, U .- H .; Lee, J .- S .; Park, Y .- Y .; Cho, S .- G .; Chintharlapalli, S .; Safe, S. Induction of the transcriptional repressor ZBTB4 in prostate cancer cells by drug-induced targeting of microRNA-17-92/106b-25 clusters. Mol. Cancer Ther. 2012, 11, 1852-1862. [CrossRef] [PubMed]
77. Dong, W .; Liu, X .; Yang, C .; Wang, D .; Xue, Y .; Ruan, X .; Zhang, M .; Song, J .; Cai, H .; Zheng, J .; et al. Glioma glycolipid metabolism: MSI2-SNORD12B-FIP1L1-ZBTB4 feedback loop as a potential treatment target. Clin. Transl. Med. 2021, 11, e411. [CrossRef]
78. Xiang, T .; He, K .; Wang, S .; Chen, W .; Li, H. Expression of zinc finger and BTB domain-containing 4 in colorectal cancer and its clinical significance. Cancer Manag. Res. 2020, 12, 9621-9626. [CrossRef] [PubMed]
79. Kim, K .; Chadalapaka, G .; Lee, S .- O .; Yamada, D .; Sastre-Garau, X .; Defossez, P .- A .; Park, Y .- Y .; Lee, J .- S .; Safe, S. Identification of oncogenic microRNA-17-92/ZBTB4/specificity protein axis in breast cancer. Oncogene 2012, 31, 1034-1044. [CrossRef] [PubMed]
80. Yu, Y .; Shang, R .; Chen, Y .; Li, J .; Liang, Z .; Hu, J .; Liu, K .; Chen, C. Tumor suppressive ZBTB4 inhibits cell growth by regulating cell cycle progression and apoptosis in Ewing sarcoma. Biomed. Pharmacother. 2018, 100, 108-115. [CrossRef]
81. Yuan, L .; Qian, G .; Chen, L .; Wu, C .- L .; Dan, H.C .; Xiao, Y .; Wang, X. Co-expression network analysis of biomarkers for adrenocortical carcinoma. Front. Genet. 2018, 9, 328. [CrossRef]