EJE

Steroidogenic factor-1 regulates a core set of target genes to promote malignancy in adrenocortical carcinoma

João C.D. Muzzi,1,2 Carmen Ruggiero,3,4 Mabrouka Doghman-Bouguerra,3,4 Maísa E. Colodel,1 Jessica M. Magno, 1,2 Jean S.S. Resende, 1,2 Nelly Durand,3,4 Juliana F. de Moura,5

Larissa M. Alvarenga,5 Luciane R. Cavalli, 1,6 Bonald C. Figueiredo, 1,7 Enzo Lalli,3,4,8,* [D and Mauro A.A. Castro2,*( 2,* İD

1Pelé Pequeno Príncipe Research Institute, Curitiba 80250-060, Brazil

2Bioinformatics and Systems Biology Laboratory, Federal University of Paraná (UFPR), Curitiba 81520-260, Brazil

3Institut de Pharmacologie Moleculaire et Cellulaire CNRS UMR 7275, Valbonne 06560, France

4Universite Côte d’Azur, Valbonne 06560, France

5Department of Basic Pathology, Federal University of Paraná (UFPR), Curitiba 81530-990, Brazil

6Department of Oncology, Lombardi Comprehensive Cancer Center, Georgetown University, Washington, DC 20057, United States

7Faculdades Pequeno Príncipe, Curitiba, PR 80230-020, Brazil

8Inserm, Valbonne 06560, France

*Corresponding author: Institut de Pharmacologie Moléculaire et Cellulaire CNRS-Inserm-Université Côte d’Azur 660 route des Lucioles 06560 Valbonne, France. Email: ninino@ipmc.cnrs.fr (E.L.); Bioinformatics and Systems Biology Laboratory, Federal University of Paraná (UFPR) Centro Politécnico UFPR (SEPT - BLOCO B) 1225 Alcides Vieira Arcoverde Street Jardim das Américas 81520-260 Curitiba, Brazil. Email: mauro.castro@ufpr.br (M.A.A.C.)

Abstract

Objective: Gene dosage is at the core of the biological activity of the steroidogenic factor-1 (SF-1/NR5A1) transcription factor. Its overexpression in adrenocortical carcinoma (ACC) is associated with enhanced proliferation and invasive capacities, steroid modulation, immune suppression, and poor prognosis. Surprisingly, 3 independent studies showed less than 10% agreement in identifying SF-1-regulated genes in the same ACC cell line, raising concerns about technical reproducibility and methodological consistency. This study aimed to reconcile discrepancies in SF-1- regulated gene identification across independent studies using a systematic approach.

Design and Methods: We reanalysed datasets from those studies using an in silico SF-1 regulon obtained from ACC TCGA data as an external reference to evaluate transcriptional patterns. Additionally, we assessed how threshold selection impacts the overlap between experiments and optimized this process. Furthermore, we performed functional experiments to evaluate how variations in SF-1 dosage impact target gene expression.

Results: Our analysis revealed comparable transcriptional patterns across all studies, reconciling transcriptional signatures and phenotypes. Threshold optimization identified consensus sets of genes responsive to SF-1 perturbations. Functional experiments confirmed that variations in SF-1 dosage significantly impact gene expression, explaining discrepancies in previous studies, and evidenced negative autoregulation of the SF-1 transcript by its encoded protein both in ACC cells and in a mouse model of Sf-1 overexpression in the adrenal cortex.

Conclusions: Our findings deepen our understanding of SF-1 regulatory activity in ACC and demonstrate that dosage is critical for observed gene expression patterns. Our integrative approach improves reproducibility and biological interpretation, offering a framework to reconcile cross-study findings.

Keywords: SF-1/NR5A1, adrenocortical carcinoma, gene regulation, bioinformatics, threshold optimization, reproducibility, regulon

Significance

Three independent previously published studies surprisingly showed only limited overlap in SF-1 target genes in the same adrenocortical carcinoma (ACC) cell line. Here, we present a unified analytical approach aimed to reconcile divergences be- tween those studies. Consensus between cross-study findings was obtained by threshold optimization and regulon analysis. SF-1 dosage contributed to the observed heterogeneity in SF-1 responsive genes. Remarkably, functional experiments showed that SF-1 is able to regulate expression of its own transcript in ACC cells. These findings have important implications in our understanding of the transcriptional regulatory mechanisms implicated in the pathogenesis of adrenocortical tumours.

Introduction

For several transcription factors (TFs), gene dosage is a major determinant of their biological function, and nuclear receptors

are overrepresented among the TFs that are more sensitive to gene dosage effects.1 Steroidogenic factor-1 (SF-1; NR5A1) is a nuclear receptor TF that has a pivotal role in the development

of adrenal glands and gonads and in controlling steroid hor- mone production, being also implicated in the pathogenesis of adrenocortical tumours. The first indication of the import- ance of SF-1 dosage for its biological function was provided by genetic inactivation experiments in mice.2 Sensitivity to SF-1 dosage also exists in humans, where heterozygote NR5A1 mutations are most commonly associated with go- nadal, but not adrenal, defects.3 SF-1 overexpression enhances adrenocortical carcinoma (ACC) cell proliferation and inva- sion and is linked to the development of adrenocortical tu- mours in mice.4-6 In humans, SF-1 overexpression is widely observed in paediatric adrenocortical tumours,7 while in adults it is associated with increased proliferation, enhanced steroido- genesis, immune suppression, and poor prognosis in ACC.8-11

Given its critical role in ACC, 3 independent studies aimed to identify the regulatory targets of SF-1: Ferraz-de-Souza et al.,12 Ehrlund et al.,13 and Doghman et al.14 All 3 studies utilized H295R cells, the “gold standard” cell line in ACC research, to evaluate differentially expressed genes (DEGs) following SF-1 dosage perturbation via knockdown or overexpression ex- periments. However, subsequently, Ruggiero et al.15 found a poor overlap when comparing results among those 3 stud- ies, with less than 10% of DEGs being shared. That diver- gence was supposed to be due to differences in experimental and analytical methods, as well as to variations in the precise levels of SF-1 dosage perturbation produced in the different experiments. 15

While a lack of agreement among studies raises questions about the reliability and consistency of results in identifying TF targets in ACC, it may also offer an opportunity to explore reconciliation strategies. Variations in experimental techni- ques, such as different methods of SF-1 knockdown or overex- pression, can influence the results and contribute to the lack of consistency between studies. It may also reflect the complexity of gene regulation in ACC, perhaps influenced by multiple TFs cooperating in a network. Another important aspect to con- sider is the biological heterogeneity of ACC cells. This hetero- geneity can result in variable gene expression profiles among samples, impacting the statistical power and thresholds used in each study, making it difficult to identify a consensus SF-1 target set.

Here, we aimed to address the discrepancies among the pre- vious findings and propose a unified approach to reconcile them. We reanalysed gene expression raw data and applied standardized analytical methods across the 3 datasets,12-14 in- cluding normalization techniques, data preprocessing, gene annotations, and statistical tests. Using a previously defined in silico SF-1 regulon obtained from ACC TCGA data16 as an external reference for comparisons, we demonstrated that, despite the apparent divergences, those experiments ex- hibit consistent phenotypes. Additionally, we evaluated how threshold selection impacts the overlap between experiments, and we propose a new threshold optimization approach to identify common gene sets between expression phenotypes. Finally, using functional experimental approaches, we demon- strate that SF-1 dosage contributes to biological heterogeneity, producing different expression levels of its target genes. We also show that SF-1 autoregulates its own transcript in ACC cells and in a mouse model of adrenocortical neoplasia caused by SF-1 overexpression. Altogether, these findings contribute to explain the previously observed discrepancies and offer an integrative approach to exploring SF-1 gene regulation in ACC.

Materials and methods

Gene expression data

We used the R package GEOquery v 2.70.0 to retrieve pheno- typic and expression data published by Ferraz-de-Souza et al.12 (GSE26529) and Doghman et al.14 (GSE43035) for SF-1 knockdown and overexpression experiments. Both studies made the gene expression matrix available in log2 GC-RMA ex- pression values. For experiments that utilized the Affymetrix Human Gene 1.0 ST Array, we reannotated the probes using the R package hugene10sttranscriptcluster.db v 8.8.0.17 For ex- periments that employed the Affymetrix Human Genome U133 Plus 2.0 Array, probe reannotation was carried out using the R package hgu133plus2.db v 3.13.0.18 We filtered probes that tar- geted the same genes by the coefficient of variation among the samples. Additionally, we filtered gene biotypes, retaining protein-coding genes. Differential expression phenotypes were obtained by comparing treatment and control groups as defined by the authors, using the R package limma v3.58.1.19

The raw data from Ehrlund et al.13 were obtained from the ArrayExpress data repository (accession number E-MEXP-3259). We used the standard protocol for a separate channel analysis of 2-color data described in the R package limma user’s guide.19 Background correction was applied using the normal-exponential method with an offset of 50. Normalization was done in 2 steps: first within arrays using the print-tip loess method (as described in ref. 13) and then be- tween arrays using the quantile method as described in limma’s vignette.19 We filtered gene biotypes, retaining protein-coding genes. Differential expression phenotypes were obtained by comparing treatment and control groups as defined by the au- thors, using the R package limma v3.58.1.19

NR5A1 regulon annotation

The NR5A1 (SF-1) regulon was retrieved from Muzzi et al.16 Briefly, those authors inferred a transcriptional network for 1600 transcription factors (TFs) using TCGA-ACC RNA-seq data20 from protein-coding genes and the RTN package to call TF-target interactions. The RTN package estimates mutual information between TFs and candidate target genes, filtering non-significant associations through permutation and bootstrap-based procedures. In a subsequent step, the RTN package uses the Algorithm for the Reconstruction of Accurate Cellular Networks (ARACNe) and the Data Processing Inequality (DPI) algorithms to refine the network by identifying and removing indirect interactions.25 We re- trieved the TF-target genes annotated for the NR5A1 regu- lon. The research complied with the Declaration of Helsinki.

NR5A1 regulon activity in published SF-1 knockdown and overexpression experiments

Regulon activities and associated P-values were estimated for each of the 5 SF-1 experiments (Ferraz-de-Souza et al.,12 Ehrlund et al.,13 and Doghman et al.14 SF-1 knockdown datasets; Ferraz- de-Souza et al.12 and Doghman et al.14 SF-1 overexpression data- sets) using a 2-tailed gene set enrichment analysis (GSEA2).21-23 P-values were calculated through 10 000 permutations and cor- rected for multiple tests by the Benjamini-Hochberg method.26

Threshold optimization in consensus gene set analysis We generated a pairwise ranked list of genes between 2 differential expression phenotypes, either between the SF-1

knockdown or overexpression phenotypes, ordered by abso- lute log2 fold change (log2FC) values for each experiment in- dividually. We then employed an incremental approach to identify the gene set size that maximizes the overlap between phenotypes. This involved sequentially selecting gene sets with top “n” genes from each ranked list of genes and incre- menting “n” until the entire genes in the ranked lists were evaluated. For each iteration, the size of the intersection be- tween the 2 gene sets was recorded. We then applied a lower-tail hypergeometric test to assess the overlap between the gene sets.27,28 The P-values were regularized and then adjusted using the Bonferroni method.29 The gene set size that minimized the regularized P-value was considered optimal.

For the overexpression experiments, as there were only 2 conditions, a single optimal threshold was obtained during the optimization process, which was used as the cut-off to de- fine the consensus gene set between SF-1 overexpression ex- periments. For the knockdown experiments, which involved 3 conditions, 3 thresholds were derived from pairwise com- parisons. Among these, the lowest threshold was selected as the cut-off to define the consensus gene set between SF-1 knockdown experiments. The direction of the perturbation ef- fect of each experiment was determined based on the log2FC calculated by comparing the treatment and control groups as defined by the original authors12-14 (see “Gene expression data”). A schematic workflow of the procedures used in our analysis is shown in Figure 1.

Cell culture

Parental H295R/TR cells4 were cultured in Dulbecco’s modi- fied Eagle’s medium/F12 supplemented with penicillin- streptomycin, 2% NuSerum (BD Biosciences), 1% ITS+ (BD Biosciences), and blasticidin (5 µg/mL; Cayla InvivoGen). Wild-type and AF-2 mutant H295R/TR SF-1 cells4 were cul- tured in the same medium with additional zeocin (100 µg/ mL; Cayla InvivoGen).

Immunoblots

Cells were treated for 24, 48, or 72 h with increasing concentra- tions of doxycycline (Dox; Sigma-Aldrich): 100, 300, 500, 700, 1000 ng/ml or with vehicle (ethanol) as a control. Protein ex- tracts were prepared by homogenization of cells and tissues in Laemmli buffer (50 mM Tris-HCl [pH 6.8], 10% glycerol, 2% sodium dodecyl sulphate [SDS], and 0.02% bromophenol

blue) containing 5% B-mercaptoethanol, as described.5 Immunoblots were performed as previously described.5 Primary antibodies used were rabbit polyclonal anti-SF-1 (07-618, Millipore; 1:1000) and mouse monoclonal anti-GAPDH (CB1001, Sigma-Aldrich; 1:2000). Secondary antibodies were HRP-conjugated anti-mouse (NA-931, GE HealthCare; 1:5000) and anti-rabbit (NA-934, GE HealthCare; 1:5000).

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

H295R/TR, H295R/TR SF-1, and H295R/TR SF-1 AF-2 mu- tant cells4 were treated for various lengths of time with Dox at a concentration of 1 µg/mL, and total RNA was extracted by the RNeasy Mini kit (Qiagen #74104). Reverse transcription- quantitative polymerase chain reaction was performed as de- scribed5 using TATA-binding protein (TBP) as a reference transcript. Human primer sequences are listed in Table S1. Results were calculated using the 4ACT threshold cycle method.30

SF-1 transgenic (TR) mice in the C57BL/6 background were previously described.4,6 All experimental protocols were ap- proved by the French Ministry of Research (number 00621.02) and were conducted according to French and European laws for animal experimentation. Total RNA from mouse adrenal glands was isolated according to the Chomczynski-Sacchi method with minor modifications.5 Reverse transcription-quantitative polymerase chain reaction was performed as described3 using mouse TATA-binding pro- tein (Tbp) as a reference transcript. Mouse primer sequences are listed in Table S2. Data were analysed using the Relative Expression Software Tool (REST; https://www.gene- quantification.de/rest.html).31

Northern blot

Total RNA was extracted according to the Chomczynski-Sacchi method.5 Fifteen micrograms of total RNA were loaded onto a 1% agarose/MOPS 1X/1% formaldehyde gel, submitted to electrophoresis for 3 h at 120 V, and transferred to a nylon membrane (Hybond N+, Amersham), as described.32 RNA was fixed to the membrane by baking for 30 min at 80 ℃. Probes located in the SF-1 coding region (SF-1 CDS: obtained by KpnI-EcoRI digestion of the pSG.SF1 plasmid33) and 3’ UTR (SF-1 3’UTR: obtained by PCR of pSG.SF1 with

Figure 1. Schematic workflow for the reconciliation and analysis of SF-1 perturbation studies. The workflow outlines the approach to investigate the low (<10%) overlap in DEG lists from 3 independent studies of SF-1 (Ferraz-de-Souza et al., 12 Ehrlund et al., 13 and Doghman et al.14). The strategy for this cross-study reconciliation was based on using the activity of an independently derived in silico NR5A1 regulon (from Muzzi et al.16) as a functional reference. This multi-step process enabled threshold optimization, confirmation of phenotype convergence, and the generation of robust consensus target sets.

Study 1 (Ferraz et al.)

Cross-study Reconciliation

SF1 DEGs list

Three studies targeted SF-1 regulatory genes in ACC using H295R cells

Study 2 (Doghman et al.)

in silico NR5A1 Regulon Activity (Muzzi et al.)

<10%

I

· Threshold optimization

Study 3 (Ehrlund et al.)

· Phenotype convergence

· Consensus target sets

primers 5’-GGGGATCCGGCTCCGCATGGTTCATTTT-3’ and 5’-GGCTCGAGATTTATTCCAGGGCTGTGGGG-3’ and subsequent cloning in pGEM-T Easy [Promega]) were radioactively labelled with [a-32P]-dATP (Amersham) by ran- dom priming using the DecaLabel DNA labelling kit (Thermo Fisher Scientific). Hybridization was performed in 0.5 M so- dium phosphate (pH 7.2)/5% SDS/1 mM EDTA at 65 ℃. The membrane was washed in 2x SSC/0.1% SDS, 1x SSC/0.1% SDS, and 0.1x SSC/0.1% SDS at 65 ℃ before expos- ure in a PhosphorImager (GE HealthCare).

Results

The NR5A1 regulon is activated in SF-1 overexpression and repressed in SF-1 knockdown phenotypes

A regulatory unit (regulon) comprises a transcription factor (TF) and its target genes.16,21-23 Targets are referred to as posi- tive or negative based on their expression relative to the TF: positive if their expression parallels that of the TF and negative if it diverges. Once defined, the activity of a regulon can be as- sessed in a differential expression phenotype using a 2-tailed Gene Set Enrichment Analysis (GSEA2), yielding a differential enrichment score (dES).21-23 A schematic representation of

regulon activity readouts is presented in Figure 2A. A regulon is considered activated if its positive targets are among the most highly expressed genes, and its negative targets are among the most repressed genes in the differential expression phenotype, resulting in a positive dES (Figure 2A, top). Conversely, it is considered repressed for the opposite distri- bution pattern, resulting in a negative dES (Figure 2A, middle). If no clear pattern emerges in the target distribution, the regu- lon activity is assigned as undefined (Figure 2A, bottom).

We use the NR5A1 (SF-1) regulon reconstructed from the TCGA-ACC cohort16 as an external reference to assess the ex- perimental agreement between the SF-1 perturbation results from Ferraz-de-Souza et al.,12 Ehrlund et al.,13 and Doghman et al.14 We calculated the activity of the NR5A1 regulon in the differential expression phenotypes generated from the 3 datasets (Figure 2B-C). The NR5A1 regulon pre- sented a significant activation in both SF-1 overexpression ex- periments. For the Doghman et al. dataset, the calculated dES was 0.94 (Figure 2B), and for the Ferraz-de-Souza et al. data- set, the dES was 0.87 (Figure 1C), both with P-values below .001. Conversely, all SF-1 knockdown experiments showed negative enrichments, indicating repression of the NR5A1 regulon: - 0.73 dES for the Doghman et al. dataset (Figure 2D), -0.81 dES for Ferraz-de-Souza et al. dataset

Figure 2. Activity of the NR5A1 regulon in SF-1 modulated phenotypes. (A) Schematic "toy" regulon activity readouts illustrating potential outcomes of a 2-tailed Gene Set Enrichment Analysis (GSEA2), showing expected patterns for activated, repressed, and undefined regulon activities, respectively. (B-F) NR5A1 regulon activity readouts in SF-1 overexpression and knockdown experiments: (B) SF-1 overexpression from Doghman et al., 14 (C) SF-1 overexpression from Ferraz-de-Souza et al., 12 (D) SF-1 knockdown from Doghman et al., 14 (E) SF-1 knockdown from Ferraz-de-Souza et al., 12 and (F) SF-1 knockdown from Ehrlund et al. 13 Top: Differential expression phenotypes were used to generate the ranked list of genes; the differential expression was replicated according to the authors, comparing treatment and control groups. Middle: Distribution of regulon targets on the ranked list of genes. Bottom: GSEA2 enrichment scores. P-values were calculated by permutation analysis and adjusted by the Benjamini-Hochberg method.

A

dES>0, activated

B

Phenotype

4

SF1 Overexpression (Doghman et al.)

C

Phenotype

4

2

2

SF1 Overexpression (Ferraz-de-Souza et al.)

TF

Enrichment score

0

0

-2

-2

targets

Regulon (Muzzi et al.)

negative (167) - positive (69)

Regulon (Muzzi et al.)

negative (167) - positive (66)

dES<0, repressed

0.4-

NR5A1

0.4-

NR5A1

TF

Enrichment score

Enrichment score

Enrichment score

0.2

0.2

targets

0.0

0.0

dES=0, undefined

-0.2

-0.2

Enrichment score

TF

-0.4

dES = 0.94

Adj. p-value < 0.001

-0.4

dES = 0.87

Adj. p-value < 0.001

Position in the ranked list of genes

targets

0

5000

10000

15000

0

5000

10000

15000

Position in the ranked list of genes

Position in the ranked list of genes

D

E

F

Phenotype

2

SF1 Knockdown (Doghman et al.)

Phenotype

2

SF1 Knockdown (Ferraz-de-Souza et al.)

Phenotype

3

2

SF1 Knockdown (Ehrlund et al.)

0

1

1

-2

0

0

-1

Regulon (Muzzi et al.)

negative (167) - positive (66)

-1

Regulon (Muzzi et al.)

negative (167) - positive (66)

-2

Regulon (Muzzi et al.)

negative (164) - positive (66)

0.4-

NR5A1

0.4

NR5A1

0.4-

NR5A1

Enrichment score

Enrichment score

Enrichment score

0.2

0.2

0.2

0.0

0.0

0.0

-0.2

-0.2

-0.2

-0.4.

dES = - 0.73

Adj. p-value < 0.001

-0.4.

dES = - 0.81

Adj. p-value < 0.001

-0.4

dES =- 0.67

Adj. p-value < 0.001

0

5000

10000

15000

0

5000

10000

15000

0

5000

10000

15000

Position in the ranked list of genes

Position in the ranked list of genes

Position in the ranked list of genes

(Figure 2E), and -0.67 dES for Ehrlund et al. dataset (Figure 2F), with P-values below .001. This consistent pattern across the studies indicates convergence between all phenotypes.

Threshold optimization reveals consistent overlap between SF-1 phenotypes

Having established that the Ferraz-de-Souza et al.,12 Ehrlund et al.,13 and Doghman et al.14 experiments converged to the same differential expression phenotypes, we aimed to identify the SF-1 modulated genes that are shared between the pheno- types. Given that the overlap between different experiments can be influenced by arbitrary threshold selection criteria, we applied a standardized approach to determine the optimal threshold for pairwise comparison. The results from our threshold optimization approach are shown in Figure 3. The optimal threshold is represented by the gene set size in which 2 experiments exhibit the most significant overlap. Results are

summarized in Table S3. For the SF-1 overexpression experi- ments, the optimal threshold was obtained by selecting the top 3221 genes from the ordered lists of Doghman et al. and Ferraz-de-Souza et al. datasets, with an observed overlap of 1096 genes (adjusted P-value < 10-4, hypergeometric test), of which 82.57% (905/1096) displayed the same log2FC direction in both experiments, referred as “consistency” in phenotype comparison. For the SF-1 knockdown experiments, the optimal threshold for the Doghman et al. and the Ferraz-de-Souza et al. datasets was 1854 genes, with an overlap of 421 (adjusted P-value < 10-4) and 72.68% consistency (306/421). The com- parison between the Doghman et al. and the Ehrlund et al. data- sets showed an optimal threshold of 1945 genes, with 375 overlapping (adjusted P-value < 10-4) and 86.4% consistency (324/375). Lastly, the Ferraz-de-Souza et al. and the Ehrlund et al. datasets had an optimal threshold of 1563 genes, with 220 overlapping (adjusted P-value < 10-4) and 79.55% consist- ency (175/220). The lists of shared genes are provided in Table S4, representing consensus gene sets.

Figure 3. Threshold optimization for assessing the overlap between SF-1 phenotypes. (A) Ferraz-de-Souza et al.12 and Doghman et al. 14 SF-1 overexpression experiments. (B) Ferraz-de-Souza et al.12 and Doghman et al.14 SF-1 knockdown experiments. (C) Ehrlund et al.13 and Doghman et al. 14 SF-1 knockdown experiments. (D) Ferraz-de-Souza et al.12 and Ehrlund et al.13 SF-1 knockdown experiments. Top: Pairwise differential expression (DE) phenotypes were used to generate a ranked list of genes, ordered by absolute log2 fold change (log2FC), for each experiment individually. Bottom: The red line (right y-axis) indicates the intersection between the ranked list of genes for different gene set sizes. For each gene set size, the overlap between gene lists was assessed by a hypergeometric test (left y-axis). The optimal threshold is indicated by a green circle (ie, at the lowest P-value), representing the gene set size in which the 2 experiments exhibit the most significant overlap.

A

1.0

B

1.0

DE phenotypes (absolute log2FC)

Ferraz’s SF1 Overexpression

DE phenotypes (absolute log2FC)

Ferraz’s SF1 Knockdown

0.0

2.5

0.0

2.5

Doghman’s SF1 Overexpression

Doghman’s SSF1 knockdown

0.0

0.0

Position in the ranked list of genes

Position in the ranked list of genes

Log10 P-value (Hypergeo test)

0

Genes in the intersection (x1000)

Log10 P-value (Hypergeo test)

0

Genes in the intersection (x1000)

15

15

-20

-5

10

10

-10

-40

5

5

-15

-60

0

0

0

5000

10000

15000

0

5000

10000

15000

Gene set size

Gene set size

C

D

2.5

1.0

DE phenotypes (absolute log2FC)

Doghman’s SF1 knockdown

DE phenotypes (absolute log2FC)

Ferraz’s SF1 knockdown

0.0

NO

4

A0

Ehrlund’s SF1 knockdown

Ehrlund’s SF1 knockdown

0.0

0.0

Position in the ranked list of genes

Position in the ranked list of genes

Log10 P-value (Hypergeo test)

0.0

Genes in the intersection (x1000)

Log10 P-value (Hypergeo test)

0

Genes in the intersection (x1000)

15

15

2.5

-1

10

10

-5.0

5

-2

5

7.5

0

-3

0

0

5000

10000

15000

0

5000

10000

15000

Gene set size

Gene set size

Different categories of genes are regulated by SF-1 in conditions of basal and increased dosage

Gene Ontology analysis of the overlapping genes regulated by SF-1 knockdown and overexpression in the Ferraz-de-Souza et al., Ehrlund et al., and Doghman et al. studies revealed that remarkably different categories of target genes are regu- lated by SF-1 in conditions of basal and increased dosage in ACC cells (Table S5). Only 46 SF-1 target genes were common between the knockdown and the overexpression conditions (Figure 4A). Genes involved in cholesterol/steroid metabolism were significantly enriched among the genes regulated by SF-1 knockdown, while genes involved in lipid metabolism, fatty acid metabolism, and cell adhesion were significantly enriched consequently to SF-1 overexpression (Figure 4B).

SF-1 dosage modulates the expression patterns of SF-1 target genes in ACC cells

Next, we evaluated the expression of a battery of SF-1 dosage- dependent target genes as a function of SF-1 dosage in H295R/ TR SF-1 cells, a subclone of the H295R cell line in which SF-1 overexpression can be induced in a doxycycline-dependent manner.4 We selected twelve genes involved in critical processes linked to tumour biology (eg, composition of the extracellu- lar matrix, cytoskeleton remodelling, cell proliferation and invasion, intracellular signalling, resistance to chemother- apy) whose expression was previously described to be regu- lated either positively or negatively by SF-1 overexpression in that cell line12,14: ADAMTS9, ARHGAP26, AXL, CADM3, COL18A1, DUSP6, EPAS1, FATE1, FGD4, MAPK4, NOV/ CCN3, and VAV2 (Table S6). Four of them (ADAMTS9, MAPK4, NOV/CCN3, and VAV2) were also regulated by SF-1 knockdown,12,14 and 8 of them (ADAMTS9, ARHGAP26,

CADM3, DUSP6, EPAS1, FATE1, NOV/CCN3, and VAV2) were retained in our consensus overexpression gene sets (Table S4). In Figure S1, we show that Dox treatment did not increase SF-1 expression in the H295R/TR parental cell line, which is expected as this cell line does not harbour a Dox-inducible SF-1 transgene. However, in H295R/TR SF-1 cells, in which SF-1 expression can be increased by Dox treatment, SF-1 protein expression was increased by 24, 48, and 72 h of Dox treatment at different concentrations, as also shown previously.4,34 To gain deeper insight into the regulatory kinetics of our battery of SF-1 dosage-dependent target genes, we monitored their expression by RT-qPCR at 24, 48, and 72 h after the start of Dox treatment. The results presented in Figure 5 show different expression patterns in re- sponse to treatment: (1) For genes upregulated by SF-1 over- expression, the expression of some of them (AXL, COL18A1, FATE1) was increasingly upregulated in proportion with the Dox treatment time, while others (EPAS1, CADM3, ARHGAP26) were significantly upregulated only at the lon- gest (72 h) time. Interestingly, the expression of other genes (FGD4, MAPK4) peaked at a shorter (24 or 48 h) time and decreased thereafter, while the expression of VAV2 was upre- gulated already at 24 h and remained constant over the Dox treatment time. (2) Genes downregulated by SF-1 overexpres- sion also displayed different repression kinetics, already start- ing at 24 h and being prolonged until 72 h of Dox treatment (DUSP6, ADAMTS9), or being increasingly repressed after SF-1 overexpression (NOV/CCN3). Taken together, these re- sults show that SF-1 dosage impacts the expression kinetics of its target genes, generating heterogeneous responses that may contribute to the discrepancies observed by Ruggiero et al.15 Results in Figure S2 show the expression patterns of the 12 SF-1 dosage-dependent genes in the Ferraz-de-Souza et al.,12 Ehrlund

Figure 4. Different sets of genes regulated by SF-1 knockdown and overexpression in ACC cells. (A) Gene overlap between the knockdown (green) and overexpression conditions in H295R cells. (B) Gene Ontology analysis of gene categories significantly enriched after SF-1 knockdown (circles) and overexpression (triangles). The log10 of the Benjamini-Hochberg FDR for the enrichment P-value is indicated with a colour scale.

A

KD

Over

B

Cholesterol metabolism

-log 10 (pvalue)

57

46

1050

5

4

Steroid metabolism

3

2

Fatty acid biosynthesis

class

· Knockdown

Overexpression

Size of each list

Lipid metabolism-

count

20

1096

1096

40

548

60

103

Cell adhesion -

0

KD

Over

Number of elements: specific (1) or shared by 2, 3, … lists

4

8

12

16

1107

2 (46)

1

Figure 5. Kinetics of mRNA expression of twelve SF-1 dosage-dependent target genes in ACC H295R/TR SF-1 cells. The time course of AXL, COL18A1, FATE1, EPAS1, CADM3, ARHGAP26, FGD4, MAPK4, VAV2, DUSP6, ADAMTS9, and NOV/CCN3 mRNA abundance after Dox treatment (1 ug/mL) of H295R/TR SF-1 cells. Data are shown as mean ± SEM of at least 3 independent experiments relative to basal expression. * P <. 05; ** P <. 01; *** P <. 001, analysis of variance (ANOVA) with Dunnett's post-test for multiple comparisons.

AXL

COL18A1

FATE1

relative expression (fold)

20

8-

8

*

*

15

*

6.

6

*

10

4

4

5.

2

2.

0

0

0

Dox

24 h

48 h

72 h

24 h

48 h

72 h

24 h

48 h

72 h

EPAS1

CADM3

ARHGAP26

relative expression (fold)

6-

*

10-

5-

*

*

8

4

4

6

3

2

4

2-

2

1.

0

0

0

Dox

24 h

48 h

72 h

24 h

48 h

72 h

24 h

48 h

72 h

FGD4

MAPK4

VAV2

relative expression (fold)

2.5

2.5

3.

*

2.0

2.0

*

1.5

1.5

2

1.0

1.0

1

0.5

0.5

0

0

0

Dox

24 h

48 h

72 h

24 h

48 h

72 h

24 h

48 h

72 h

DUSP6

ADAMTS9

NOV

relative expression (fold)

0.8

0.8

1.0-

0.6

0.6

0.8.

0.6

0.4

0.4

0.4

0.2

0.2

0.2.

0

0

0

Dox

24 h

48 h

72 h

24 h

48 h

72 h

24 h

48 h

72 h

et al.,13 and Doghman et al.14 datasets, indicating that all ex- pression patterns agree with our in vitro results for comparable experiments, that is, except for DUSP6, ADAMTS9, and NOV/ CCN3, all the other genes exhibit positive log2 fold changes in the SF-1 overexpression phenotypes.

Negative autoregulation of SF-1 expression in human adrenocortical cells and mouse adrenal gland

Autoregulation represents another source of variation that could introduce heterogeneity in SF-1 downstream effects. Positive SF-1 autoregulation has been demonstrated in mouse foetal adrenal development, in which this

transcription factor upregulates its own expression after ini- tiation by a Hox-Pbx1-Prep1 complex.35 To gain further in- sights into SF-1 transcriptional mechanisms, we evaluated its autoregulation in H295R/TR SF-1 cells. The overexpres- sion of a Dox-dependent SF-1 transgene negatively regulated the expression of the endogenous SF-1 mRNA (Figure 6A). Downregulation of the endogenous SF-1 transcript is de- pendent on SF-1 transcriptional activity, since it did not take place when the expression of a transcriptionally silent AF-2 SF-1 mutant4 was activated by Dox treatment in H295R/TR cells (Figure 6A). No variation of SF-1 mRNA ex- pression was caused by Dox treatment in parental H295R/TR cells (Figure 6A). Downregulation of the endogenous SF-1 transcript by SF-1 overexpression in H295R/TR SF-1 cells

Figure 6. SF-1 transcript and protein regulation in human ACC cells and in mouse adrenal gland. (A) Transgenic (red histograms), endogenous (blue histograms), and total (yellow histograms) SF-1 mRNA abundance after Dox treatment (1 µg/mL for 72 h) of H295R/TR, H295R/TR SF-1, and H295R/TR SF-1 AF2 mutant cells, respectively. Results are calculated as the ratio of expression relative to basal and shown as mean ± SEM of at least 3 independent experiments. * P <. 05; ** P <. 01; *** P < . 001, 1-way analysis of variance (ANOVA) with Dunnett's post-test for multiple comparisons. (B) Total RNA was extracted from H295R/TR SF-1 cells under basal or Dox (1 µg/mL, 72 h treatment) conditions. The Northern blots (15 µg total RNA per lane) were probed with labelled cDNA probes for the SF-1 CDS region or 3'UTR. Gel migration of 28S and 18S rRNA is indicated. (C) Expression of total (Sf-1 pan#1, Sf-1 pan#2) and endogenous (Sf-1 mouse) Sf-1 mRNA in Sf-1 transgenic (TR) mouse adrenal glands (n=26) calculated by the Relative Expression Software Tool (REST). The control group is represented by wild-type (WT) mice. * P <. 05; *** P <. 001, Pairwise Fixed Reallocation Randomization Test. (D) Immunoblot showing SF-1 protein levels in adrenal glands of 6-week-old WT or Sf-1 TR mice. GAPDH was used as a loading control. SF-1 protein abundance was quantified after normalization to GAPDH about the expression level in WT males and shown as mean ± SEM. n=4. **** P <. 0001, 1-way analysis of variance (ANOVA) with Šídák's post-test for multiple comparisons.

A

B

CDS

3’ UTR

SF-1 mRNA expression (fold)

64


transgene

32

**

endogenous

28S

16

**

total

8

4

2

*

1

*

18S

0.5

TR

WT

AF2 mut

Dox

-

+

-

+

C

WT

Sf-1 TR

C

*

SF-1

Sf-1 mRNA expression

8

GAPDH

male

female

male

female

4

*


2

2.0-



M WT

1

SF-1/GAPDH

1.5

M TR

0.5

F WT

1.0

F TR

0.25

0.5

0.125

0

Sf-1 pan#1

Sf-1 pan#2

Sf-1 mouse

male

female

was confirmed by Northern blot, using probes for the SF-1 CDS and for the SF-1 3’UTR (Figure 6B).

We also studied the impact of Sf-1 overexpression in a YAC transgenic mouse model bearing multiple copies of the rat Sf-1 gene (Sf-1 TR).4,6,36 We have previously shown that Sf-1 TR mice develop adrenocortical dysplasia and tumours.4,6 By using 2 different couples of primers amplifying different fragments in both the mouse and the rat Sf-1 coding sequence (Sf-1 pan#1 and pan#2) to assay for total Sf-1 mRNA expression, we showed that its transcript was moderately, but significantly, upregulated in the adrenal glands of Sf-1 TR mice compared to wild-type (1.307 fold on average by Sf-1 pan#1 RT-qPCR; 1.476 fold on average by Sf-1 pan#2 RT-qPCR). In contrast, endogenous Sf-1 mRNA expression, as measured by a couple of primers which amplify specifically the mouse transcript (Sf-1 mouse), was significantly downregulated in TR adrenals compared to wild-type (0.660-fold on average) (Figure 6C). Conversely, RT-qPCR performed by primers amplifying specifically the transgenic rat Sf-1 transcript showed its expression at high levels only in the adrenal glands of Sf-1 TR mice (Ct = 22-24). Consistent with the mRNA results, the SF-1 protein was moderately but significantly overexpressed in the ad- renal glands of Sf-1 TR animals of both sexes (1.6-fold on average for both males and females) compared to wild type (Figure 6D). These results are consistent with the data we previously published on the quantification of SF-1 protein overexpression in the adrenal glands of Sf-1 TR mice.4

Discussion

One of the pillars of scientific methodology is reproducibility. Although biological variability and statistical noise are ex- pected sources of divergence between experiments, studies with similar experimental designs should exhibit some level of consistency. Concerns regarding the reproducibility of find- ings have gained attention in various fields, including molecu- lar biology.37 Ruggiero et al.15 brought attention to the limited overlap between SF-1 targets identified in 3 independent in vi- tro studies using the same ACC cell line model.12-14 Our study aimed to investigate whether it is possible to reconcile those apparently discordant previous studies by taking a new ap- proach to analyse existing gene expression data. The differen- ces in gene expression platforms, as well as in the methods used in these studies to produce SF-1 knockdown and overex- pression, may partially explain the observed discrepancies. Ferraz-de-Souza et al.12 employed GeneChip Human Gene 1.0 ST Arrays, while Doghman et al.14 used the same array for their knockdown experiments and GeneChip HG-U133 Plus 2.0 for their overexpression studies. Ehrlund et al.13 uti- lized the Operon Human 35k microarray. Those platforms vary in the methods in how they quantify gene expression, po- tentially contributing to the differences in intensity measurements.38-40 In addition, each study applied different normalization methods and thresholds to identify DEGs. For example, Ehrlund et al.13 focused on the top 35 up- and

downregulated genes by log2FC, whereas Ferraz-de-Souza et al.12 used an adjusted P-value <. 05, and Doghman et al.14 selected genes with absolute log2FC > 2 and unadjust- ed P-value <. 05. Each study applied different knockdown and overexpression methods, affecting SF-1 interference lev- els and, consequently, target gene expression. Despite these differences and inherent variability even within the same cell line across laboratories, 41-43 some level of concordance is expected given that all studies aimed to assess SF-1 function in the same ACC cell line model.

Here, we used the NR5A1 regulon, inferred by Muzzi et al.16 from TCGA-ACC RNA-seq data,2º as an external ref- erence to compare the 3 aforementioned studies. A regulon is defined as a regulatory module composed of a TF and its sig- nificantly associated target genes. As each target can be linked to multiple TFs, regulation can occur as a result of both direct and indirect interactions.25 The NR5A1 regulon was con- structed by applying the DPI algorithm, which removes the weakest interaction in any triplet of 2 TFs and a target gene, thus preserving the dominant TF-target pairs.25 Therefore, the composition of the filtered regulon is expected to be less complex and enriched with the most significant and directed interactions.22 Our analysis assessed the activity of the NR5A1 regulon within differential expression phenotypes and took into account the complete list of genes with their log- fold changes, not just those identified as DEGs, avoiding arbi- trary threshold choices in gene selection. The results indicated a consistent agreement between all studies, showing positive regulon activity in the SF-1 overexpression phenotypes and negative regulon activity in the SF-1 knockdown phenotypes. These findings suggest a convergence across the 3 studies, des- pite the limited overlap between the DEGs defined in each study. Notably, the NR5A1 regulon was independently de- rived from bulk RNA-seq data from ACC tumour samples, which may reflect the heterogeneity of tumour, stromal, and immune cells, in contrast with the controlled environment of in vitro cell cultures. Yet, all studies were consistently aligned by the NR5A1 regulon activity readouts. Given the conver- gence that we achieved in our phenotypic analysis, in the sub- sequent approach, we aimed to identify consensus gene lists across the studies, using a unified threshold criteria for select- ing genes. To standardize the comparisons, we implemented a threshold optimization approach, ranking genes and system- atically evaluating gene set sizes and intersections to find the threshold that maximized the significance of the overlap be- tween phenotypes. We found substantially larger gene sets than those reported by Ruggiero et al.,15 providing shared gene lists that represent SF-1 consensus target gene sets.

Once we addressed these major reproducibility issues, we in- vestigated possible sources of variation that may contribute to the biological heterogeneity of SF-1 regulation. We assessed the effect of changes in SF-1 dosage and time of exposure at differ- ent dosages. We focused on 12 genes (AXL, COL18A1, FATE1, EPAS1, CADM3, ARHGAP26, FGD4, MAPK4, VAV2, DUSP6, ADAMTS9, and NOV/CCN3), involved in critical processes linked to tumour biology, each of which was previously described12-14 to be regulated by SF-1 in a dosage-dependent manner (Table S6). We showed that SF-1 dosage impacts the kinetics of expression of its target genes, generating different regulatory patterns and responses. For ex- ample, the VAV2 gene was consistently upregulated in all SF-1 overexpression time points, in agreement with Ferraz-de-Souza et al.,12 Ehrlund et al., 13 and Doghman et al.14 where VAV2

was upregulated after SF-1 overexpression and downregulated in SF-1 knockdown experiments, while other genes, such as MAPK4, displayed a more heterogeneous response to varia- tions in SF-1 dosage over time.

An important limitation of our study is the investigation of the transcriptional response to SF-1 only in the H295R cell line. While this represents the most highly differentiated and most widely used cell model for ACC, recently a few other hu- man cell lines have been described which display varying de- grees of adrenocortical differentiation and expression of SF-1.44 We acknowledge that it would be important to also characterize dosage-dependent SF-1 target genes in those new cell lines and compare them to H295R in future studies, but in our present work, we have focused our functional stud- ies in the H295R cell line with the purpose to validate and ex- tend the results obtained in our bioinformatic analysis of 3 distinct datasets all derived from that same cell line. Also, we only tested gene expression variations in conditions of SF-1 overexpression but not knockdown. While this repre- sents another limitation of our work, SF-1 overexpression has a direct relevance to the mechanisms of ACC pathogenesis.4,6-8 However, despite these limitations, our re- sults indicate that dosage is central to the SF-1 biological activ- ity, determining different levels of up- and downregulation of target genes. The temporal aspect is also relevant, as demon- strated by varying levels of gene activation or repression at specific time points. This highlights the existence of low- and high-responder target genes according to their modulation de- gree in response to the same SF-1 dosage. Notably, our cell cul- ture and transgenic mouse experiments show that SF-1 negatively autoregulates its own transcript. This autoregula- tion underscores the fine-tuning responses to SF-1 activity. Genetic inactivation experiments in mice provided the first in- dications of SF-1 dosage sensitivity,2,45 a characteristic that also exists in humans. Collectively, our findings contribute to the understanding that SF-1 functions as a “Goldilocks” transcription factor,11 which modulates different sets of target genes according to its dosage.

In conclusion, our results reconcile the apparent discrepan- cies between the 3 in vitro studies12-14 and provide insights into the complexity of SF-1 transcriptional regulatory mecha- nisms in ACC cells. These findings highlight the need for care- ful control of SF-1 levels in future studies, as dosage variations can significantly impact experimental results. Our approach generated unified lists of shared SF-1 responsive genes between experiments, which can be further explored in future research. The use of an in silico regulon to support in vitro findings of- fers a robust method for reconciling transcriptional signatures and phenotypes. This approach could be applied to other TFs and datasets in gene regulatory studies.

Acknowledgments

We thank John Achermann for help in retrieving data from the Ferraz-de-Souza et al. study and Leslie Heckert for pro- viding Sf-1 transgenic mice. This article is based on work from COST Action CA20122 Harmonis@tion, supported by COST.

Supplementary material

Supplementary material is available at European Journal of Endocrinology online.

Funding

This work was funded by Agence Nationale de la Recherche (ANR), grant number ANR20-CE14-0007 (Goldilocks) to E.L., CNRS IIPACT International Research Network to E.L. and B.C.F., and the BIG DATA innovation programme funded by Behring Foundation. The authors would like to ac- knowledge the support of the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), funded through CNPq grants (316622/2021-4 and 440412/2022-6 to M.A.A.C; 310263/2022-0 to J.F.d.M). M.A.A.C. is also funded by CAPES (88882.632783/2021-2001) and Fundação Araucária (NAPI Bioinformática).

Authors’ contributions

João C.D. Muzzi (Conceptualization [equal], Data curation [equal], Formal analysis [equal], Methodology [equal], Visualization [equal], Writing-original draft [equal]), Carmen Ruggiero (Methodology [equal], Writing-review & editing [equal]), Mabrouka Doghman-Bouguerra (Methodology [equal], Writing-review & editing [equal]), Maísa E. Colodel (Methodology [equal]), Jessica M. Magno (Visualization [equal], Writing-review & editing [equal]), Jean S.S. Resende (Visualization [equal], Writing-review & editing [equal]), Nelly Durand (Methodology [equal], Writing-review & editing [equal]), Juliana F. de Moura (Visualization [equal], Writing-review & editing [equal]), Larissa M. Alvarenga (Visualization [equal], Writing- review & editing [equal]), Luciane R. Cavalli (Visualization [equal], Writing-review & editing [equal]), Bonald C. Figueiredo (Conceptualization [equal], Funding acquisi- tion [equal], Methodology [equal], Supervision [equal]), Enzo Lalli (Conceptualization [equal], Formal analysis [equal], Funding acquisition [equal], Methodology [equal], Supervision [equal], Writing-original draft [equal], Writing-review & editing [equal]), and Mauro A.A. Castro (Conceptualization [equal], Funding acquisition [equal], Methodology [equal], Supervision [equal], Visualization [equal], Writing-review & editing [equal])

Data availability

All datasets and software used in this study are publicly access- ible through the sources detailed in the Methods section.

References

1. Ni Z, Zhou XY, Aslam S, Niu DK. Characterization of human dosage-sensitive transcription factor genes. Front Genet. 2019;10: 1208. https://doi.org/10.3389/fgene.2019.01208

2. Luo X, Ikeda Y, Parker KL. A cell-specific nuclear receptor is es- sential for adrenal and gonadal development and sexual differen- tiation. Cell. 1994;77(4):481-490. https://doi.org/10.1016/0092- 8674(94)90211-9

3. Domenice S, Machado AZ, Ferreira FM, et al. Wide spectrum of NR5A1-related phenotypes in 46,XY and 46,XX individuals. Birth Defects Res C Embryo Today. 2016;108(4):309-320. https://doi.org/10.1002/bdrc.21140

4. Doghman M, Karpova T, Rodrigues GA, et al. Increased steroido- genic factor-1 dosage triggers adrenocortical cell proliferation and cancer. Mol Endocrinol. 2007;21(12):2968-2987. https://doi.org/ 10.1210/me.2007-0260

5. Ruggiero C, Doghman-Bouguerra M, Sbiera S, et al. Dosage- dependent regulation of VAV2 expression by steroidogenic factor-1

drives adrenocortical carcinoma cell invasion. Sci Signal. 2017; 10(469):eaal2464. https://doi.org/10.1126/scisignal.aal2464

6. de Late PL, El Wakil A, Jarjat M, et al. Vanin-1 inactivation antag- onizes the development of adrenocortical neoplasia in Sf-1 trans- genic mice. Endocrinology. 2014;155(7):2349-2354. https://doi. org/10.1210/en.2014-1034

7. Pianovski MA, Cavalli LR, Figueiredo BC, et al. SF-1 overexpres- sion in childhood adrenocortical tumours. Eur J Cancer. 2006; 42(8):1040-1043. https://doi.org/10.1016/j.ejca.2005.12.032

8. Sbiera S, Schmull S, Assie G, et al. High diagnostic and prognostic value of steroidogenic factor-1 expression in adrenal tumors. J Clin Endocrinol Metab. 2010;95(10):E161-E171. https://doi.org/ 10.1210/jc.2010-0584

9. Lalli E. Adrenocortical development and cancer: focus on SF-1. J Mol Endocrinol. 2010;44(6):301-307. https://doi.org/10.1677/ JME-10-0013

10. Muzzi JCD, Magno JM, Cardoso MA, et al. Adrenocortical carcin- oma steroid profiles: in silico pan-cancer analysis of TCGA data un- covers immunotherapy targets for potential improved outcomes. Front Endocrinol (Lausanne). 2021;12:672319. https://doi.org/ 10.3389/fendo.2021.672319

11. Relav L, Doghman-Bouguerra M, Ruggiero C, et al. Steroidogenic Factor 1, a Goldilocks transcription factor from adrenocortical or- ganogenesis to malignancy. Int J Mol Sci. 2023;24(4):3585. https:/ doi.org/10.3390/ijms24043585

12. Ferraz-de-Souza B, Hudson-Davies RE, Lin L, et al. Sterol O-acyltransferase 1 (SOAT1, ACAT) is a novel target of steroido- genic factor-1 (SF-1, NR5A1, Ad4BP) in the human adrenal. J Clin Endocrinol Metab. 2011;96(4):E663-E668. https://doi.org/ 10.1210/jc.2010-1853

13. Ehrlund A, Jonsson P, Vedin LL, et al. Knockdown of SF-1 and RNF31 affects components of steroidogenesis, TGFß, and Wnt/ B-catenin signaling in adrenocortical carcinoma cells. PLOS One. 2012;7(3):e32080. https://doi.org/10.1371/journal.pone.0032080

14. Doghman M, Figueiredo BC, Volante M, et al. Integrative analysis of SF-1 transcription factor dosage impact on genome-wide binding and gene expression regulation. Nucleic Acids Res. 2013;41(19): 8896-8907. https://doi.org/10.1093/nar/gkt660

15. Ruggiero C, Doghman M, Lalli E. How genomic studies have im- proved our understanding of the mechanisms of transcriptional regulation by NR5A nuclear receptors. Mol Cell Endocrinol. 2015;405:138-144. https://doi.org/10.1016/j.mce.2015.02.008

16. Muzzi JCD, Magno JM, Souza JS, et al. Comprehensive charac- terization of the regulatory landscape of adrenocortical carcin- oma: novel transcription factors and targets associated with prognosis. Cancers (Basel). 2022;14(21):5279. https://doi.org/ 10.3390/cancers14215279

17. MacDonald JW. hugene10sttranscriptcluster.db: Affymetrix hu- gene10 annotation data (chip hugene10sttranscriptcluster). R package version 8.7.0; 2017. https://doi.org/10.18129/B9.bioc. hugene 10sttranscriptcluster.db

18. Carlson M. hgu133plus2.db: Affymetrix Human Genome U133 Plus 2.0 Array annotation data (chip hgu133plus2). R package version 3.2.3; 2016. https://doi.org/10.18129/B9.bioc.hgu133plus2.db

19. Ritchie ME, Phipson B, Wu D, et al. Limma powers differential expres- sion analyses for RNA-Sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. https://doi.org/10.1093/nar/gkv007

20. Zheng S, Cherniack AD, Dewal N, et al. Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer Cell. 2016; 29(5):723-736. https://doi.org/10.1016/j.ccell.2016.04.002

21. Castro M, de Santiago I, Campbell T, et al. Regulators of genetic risk of breast cancer identified by integrative network analysis. Nat Genet. 2016;48(1):12-121. https://doi.org/10.1038/ng.3458

22. Fletcher M, Castro M, Wang X, et al. Master regulators of FGFR2 signalling and breast cancer risk. Nat Commun. 2013;4(1):2464. https://doi.org/10.1038/ncomms3464

23. Campbell TM, Castro MAA, de Santiago I, et al. FGFR2 risk SNPs confer breast cancer risk by augmenting oestrogen

responsiveness. Carcinogenesis. 2016;37(8):741-750. https:// doi.org/10.1093/carcin/bgw050

24. Groeneveld CS, Chagas VS, Jones SJM, et al. RTNsurvival: an R/ Bioconductor package for regulatory network survival analysis. Bioinformatics. 2019;35(21):4488-4489. https://doi.org/10.1093/ bioinformatics/btz269

25. Margolin AA, Wang K, Lim WK, et al. Reverse engineering cellular networks. Nat Protoc. 2006;1(2):662-671. https://doi.org/10.1038/ nprot.2006.106

26. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). 1995;57(1):289-300. https://doi.org/10. 1111/j.2517-6161.1995.tb02031.x

27. Fisher RA. On the interpretation of x2 from contingency tables, and the calculation of P. J R Stat Soc. 1922;85(1):87-94. https://doi.org/ 10.2307/2340521

28. Rivals I, Personnaz L, Taing L, Potier MC. Enrichment or deple- tion of a GO category within a class of genes: which test? Bioinformatics. 2007;23(4):401-407. https://doi.org/10.1093/ bioinformatics/btl633

29. Bonferroni CE. Teoria statistica delle classi e calcolo delle probabi- lità. Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze. 1936;8:3-62.

30. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-AACT method. Methods. 2001;25(4):402-408. https://doi.org/10.1006/meth.2001.1262

31. Pfaffl MW, Horgan GH, Dempfle L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30(9):e36. https://doi.org/10.1093/nar/30.9.e36

32. Sambrook J, Russell DW. Molecular Cloning: a Laboratory Manual. IIId ed. Cold Spring Harbor Laboratory Press; 2001.

33. Lehmann SG, Lalli E, Sassone-Corsi P. X-linked adrenal hypoplasia congenita is caused by abnormal nuclear localization of the DAX-1 protein. Proc Natl Acad Sci U S A. 2002;99(12):8225-8230. https:/ doi.org/10.1073/pnas.122044099

34. Ye P, Nakamura Y, Lalli E, Rainey WE. Differential effects of high and low Steroidogenic Factor-1 expression on CYP11B2 ex- pression and aldosterone production in adrenocortical cells. Endocrinology. 2009;150(3):1303-1309. https://doi.org/10.1210/ en.2008-1110

35. Zubair M, Ishihara S, Oka S, et al. Two-step regulation of Ad4BP/ SF-1 gene transcription during fetal adrenal development: initiation by a Hox-Pbx1-Prep1 complex and maintenance via autoregulation by Ad4BP/SF-1. Mol Cell Biol. 2006;26(10):4111-4121. https://doi. org/10.1128/MCB.02279-07

36. Karpova T, Presley J, Manimaran RR, et al. A Ftz-F1-containing yeast artificial chromosome recapitulates expression of Steroidogenic Factor 1 in vivo. Mol Endocrinol. 2005;19(10):2549-2563. https:// doi.org/10.1210/me.2005-0177

37. Baker M. 1,500 scientists lift the lid on reproducibility. Nature. 2016;533(7604):452-454. https://doi.org/10.1038/533452a

38. Maouche S, Poirier O, Godefroy T, et al. Performance comparison of two microarray platforms to assess differential gene expression in human monocyte and macrophage cells. BMC Genomics. 2008; 9(1):302. https://doi.org/10.1186/1471-2164-9-302

39. Petersen D, Chandramouli GV, Geoghegan J, et al. Three micro- array platforms: an analysis of their concordance in profiling gene expression. BMC Genomics. 2005;6(1):63. https://doi.org/10. 1186/1471-2164-6-63

40. Pedotti P, ‘t Hoen PA, Vreugdenhil E, et al. Can subtle changes in gene expression be consistently detected with different microarray platforms? BMC Genomics. 2008;9(1):124. https://doi.org/10. 1186/1471-2164-9-124

41. Ben-David U, Siranosian B, Ha G, et al. Genetic and transcriptional evolution alters cancer cell line drug response. Nature. 2018; 560(7718):325-330. https://doi.org/10.1038/s41586-018-0343-0

42. Kleensang A, Vantangoli M, Odwin-DaCosta S, et al. Genetic vari- ability in a frozen batch of MCF-7 cells invisible in routine authen- tication affecting cell function. Sci Rep. 2016;6(1):28994. https:// doi.org/10.1038/srep28994

43. Kinker GS, Greenwald AC, Tal R, et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity. Nat Genet. 2020;52(11):1208-1218. https://doi.org/10.1038/ s41588-020-00713-z

44. Sigala S, Rossini E, Abate A, et al. An update on adrenocortical cell lines of human origin. Endocrine. 2022;77(3):432-437. https://doi. org/10.1007/s12020-022-03112-w

45. Bland ML, Jamieson CA, Akana SF, et al. Haploinsufficiency of Steroidogenic Factor-1 in mice disrupts adrenal development lead- ing to an impaired stress response. Proc Natl Acad Sci U S A. 2000;97(26):4488-4493. https://doi.org/10.1073/pnas.080580497