Integrated genomic characterization of adrenocortical carcinoma

Guillaume Assié1-4,22, Eric Letouzé5,22, Martin Fassnacht6-8, Anne Jouinot1-3, Windy Luscap1-3, Olivia Barreau1-4, Hanin Omeiri1-3, Stéphanie Rodriguez1-3, Karine Perlemoine1-3, Fernande René-Corail1-3, Nabila Elarouci5, Silviu Sbiera6,7, Matthias Kroiss8, Bruno Allolio7, Jens Waldmann9, Marcus Quinkler10, Massimo Mannelli11, Franco Mantero12, Thomas Papathomas13, Ronald De Krijger13, Antoine Tabarin14,15, Véronique Kerlan15,16, Eric Baudin15,17, Frédérique Tissier1-3,18, Bertrand Dousset1-4,19, Lionel Groussin1-4, Laurence Amar20, Eric Clauser21, Xavier Bertagna1-4,15, Bruno Ragazzon1-3, Felix Beuschlein6, Rossella Libé1-4,15, Aurélien de Reyniès5,23 & Jérôme Bertherat1-4,15,23

npg

Adrenocortical carcinomas (ACCs) are aggressive cancers originating in the cortex of the adrenal gland1. Despite overall poor prognosis, ACC outcome is heterogeneous2,3. We performed exome sequencing and SNP array analysis of 45 ACCs and identified recurrent alterations in known driver genes4,5 (CTNNB1, TP53, CDKN2A, RB1 and MEN1) and in genes not previously reported in ACC (ZNRF3, DAXX, TERT and MED12), which we validated in an independent cohort of 77 ACCs. ZNRF3, encoding a cell surface E3 ubiquitin ligase6, was the most frequently altered gene (21%) and is a potential new tumor suppressor gene related to the ß-catenin pathway. Our integrated genomic analyses further identified two distinct molecular subgroups with opposite outcome. The C1A group of ACCs with poor outcome displayed numerous mutations and DNA methylation alterations, whereas the C1B group of ACCs with good prognosis displayed specific deregulation of two microRNA clusters. Thus, aggressive and indolent ACCs correspond to two distinct molecular entities driven by different oncogenic alterations.

Tumors of the adrenal cortex are common (with a prevalence of 1-10%), but most are benign adenomas discovered incidentally7. In contrast, ACCs are rare (with an annual incidence of 0.7-2.0 cases per million8,9) and aggressive, with a 5-year survival rate lower than

35% in most series2,10-12. Currently, surgery is the only curative therapy available. Medical treatments, including the adrenolytic drug mitotane and cytotoxic chemotherapy, show limited therapeutic potential13. Unraveling the genomic alterations underlying ACC is crucial for the development of more effective therapies and to help predict individual outcomes. Study of familial forms of ACC identified TP53-inactivating mutations and IGF2 overexpression as key drivers14. CTNNB1- activating mutations are also common, occurring in about one-fifth of cases15. Transcriptome studies16-18 have shown the existence of two molecular subgroups strongly associated with survival, which we named C1A (poor prognosis) and C1B (good prognosis)16. Alterations in copy number19-22, DNA methylation23-25 and microRNA (miRNA) expression26-31 associated with malignancy and prognosis have also been described. However, integrative and extensive molecular charac- terization of the disease is needed to grasp the full spectrum of driver genes and altered pathways. We have now analyzed a discovery cohort of 45 ACCs through a combination of genomic approaches, including exome sequencing, SNP arrays, DNA methylation analysis, mRNA expression arrays and miRNA sequencing (Supplementary Table 1). We subsequently validated candidate driver genes by targeted sequenc- ing and SNP arrays in an independent cohort of 77 ACCs, recruited from 16 clinical centers from the European Network for the Study of Adrenal Tumors (ENSAT; Supplementary Table 1).

1INSERM U1016, Institut Cochin, Paris, France. 2CNRS UMR 8104, Paris, France. 3Université Paris Descartes, Sorbonne Paris Cité, Paris, France. 4Center for Rare Adrenal Diseases, Department of Endocrinology, Assistance Publique-Hôpitaux de Paris, Hôpital Cochin, Paris, France. 5Programme Cartes d’Identité des Tumeurs (CIT), Ligue Nationale Contre Le Cancer, Paris, France. 6Medizinische Klinik und Poliklinik IV, Klinikum der Universität München, University of Munich, Munich, Germany. 7Endocrine and Diabetes Unit, Department of Internal Medicine I, University Hospital of Würzburg, Würzburg, Germany. 8Comprehensive Cancer Center Mainfranken, University of Würzburg, Würzburg, Germany. 9Visceral, Thoracic and Vascular Surgery, University Hospital Giessen and Marburg, Marburg, Germany. 10Department of Clinical Endocrinology, Charité Campus Mitte, Charité University Medicine, Berlin, Germany. 11Department of Experimental and Clinical Biomedical Sciences, University of Florence, Florence, Italy. 12Endocrinology Unit, Department of Medicine, University of Padova, Padova, Italy. 13Department of Pathology,

Josephine Nefkens Institute, Erasmus MC University Medical Center, Rotterdam, The Netherlands. 14Department of Endocrinology, Diabetes and Metabolic Diseases, University Hospital of Bordeaux, Bordeaux, France. 15Rare Adrenal Cancer Network COMETE, Paris, France. 16Department of Endocrinology, Diabetes and Metabolic Diseases, University Hospital of Brest, Brest, France. 17Department of Nuclear Medicine and Endocrine Oncology, Institut Gustave Roussy, Université Paris-Sud, Villejuif, France. 18Department of Pathology, Assistance Publique-Hôpitaux de Paris, Hôpital Pitié-Salpétrière, Pierre et Marie Curie Université, Paris, France.

19Department of Digestive and Endocrine Surgery, Assistance Publique-Hôpitaux de Paris, Hôpital Cochin, Paris, France. 20Hypertension Unit, Assistance Publique-Hôpitaux de Paris, Hôpital Européen Georges Pompidou, Paris, France. 21Oncogenetic Laboratory, Assistance Publique-Hôpitaux de Paris, Hôpital Cochin, Paris, France. 22These authors contributed equally to this work. 23These authors jointly directed this work. Correspondence should be addressed to J.B. (jerome.bertherat@cch.aphp.fr).

Received 8 November 2013; accepted 17 March 2014; published online 20 April 2014; doi:10.1038/ng.2953

Figure 1 Pattern of somatic copy number alterations in ACC. (a) Frequency of gains, losses, LOH, high-level amplifications and homozygous deletions along the genome. The left y axis indicates the frequency of low-amplitude changes (gains, losses and LOH); the right y axis indicates the frequency of high-amplitude changes (amplifications and homozygous deletions). Labels indicate regions with amplifications or homozygous deletions in >5% of cases. (b) Heat maps of copy number profiles in three regions with recurrent homozygous deletions but no established tumor suppressor gene. The absolute copy number profile of each sample harboring a homozygous deletion in each region is represented by color.

a

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18 19| 20 |21 22|

60%

10%

Gains

30%

5%

High-level amplifications

0%

0%

Losses

30%

5%

Homozygous deletions

LOH

60%

10%

3q13.31

4q34.3

TERT

CDKN2A

RB1

22q12.1

b

3q13.1

4q34.3

22q12.1

Absolute copy number

>2

2

1

0

GAP43

LOC285194

LINC00290

MGC45800

CHEK2

ZNRF3

KREMEN1

HH

1

-

I

TNEM3

I

I

LSAMP

MIR1305

XBP1

C22orf31

npg

Somatic copy number alterations were profiled by SNP arrays in 121 ACCs. Chromosomal instability, measured as the fraction of aberrant arms (FAA score)32, was very high in most samples (median FAA = 0.85, range of 0.2-1). Chromosomes 4, 5, 7, 8, 12, 16, 19 and 20 showed signif- icant gains, whereas chromosome arm 1p and chromosome 22 showed significant deletion (q value < 0.05; Fig. la and Supplementary Table 2). Loss of heterozygosity (LOH) was very frequent, with 16 of 22 auto- somes showing LOH in ≥30% of cases. As previously reported33, LOH encompassing the IGF2 locus at 11p15 was very frequent (82%) and was mostly copy neutral. Most tumors were either hypodiploid (33%) with a large number of chromosome losses or polyploid (43%) with a large number of copy-neutral LOH events, a pattern that may be explained by the accumulation of deletions followed by cell fusion or endorepli- cation (Supplementary Fig. 1). We identified significantly recurrent high-level amplifications of TERT (5p15.33) and CDK4 (12q14) and homozygous deletions of CDKN2A (9p21.3), RB1 (13q14) and ZNRF3 (22q12.1) (qvalue <0.05; Fig. 1 and Supplementary Table 3). Recurrent homozygous deletions were also found at 3q13.1 and 4q34.3 (Fig. 1b). The 3q13.1 region is frequently deleted in osteosarcoma34, and the long noncoding RNA LOC285194 (encoded by TUSC7) was shown to be a p53-regulated tumor suppressor35. Deletions in the 4q34.3 region had no single region of overlap, but the deletion peak was centered on the long noncoding RNA locus LINC00290, as also described in a recent pan-cancer study36. These findings suggest that the deregulation of long noncoding RNAs, whose implication in cancer is increasingly being appreciated37,38, may be important in ACC.

We performed exome capture followed by paired-end massively parallel sequencing on genomic DNA from 45 ACCs (discovery cohort) and matched normal samples. The baits for exon cap- ture targeted 335,765 exons from 20,975 genes totaling 70 Mb of

non-redundant sequence. The mean coverage depth was 90x per sam- ple, with 84% of targets covered to a depth of ≥25x (Supplementary Fig. 2). We identified 3,153 somatic mutations in the 45 tumors (Supplementary Table 4), 1,881 of which occurred in 2 tumors with a hypermutation phenotype. Excluding these highly mutated sam- ples, the remaining tumors harbored a median of 6 silent and 24 non-silent mutations (range of 3 to 71 total mutations), correspond- ing to a mean somatic mutation rate in coding sequences of 0.60 mutations per megabase (Supplementary Fig. 3a). Having higher numbers of mutations was associated with worse 5-year survival rate (P = 0.003, Wilcoxon rank-sum test), higher Weiss score (P = 0.02, ANOVA) and higher ENSAT stage (P = 0.06, ANOVA; Fig. 2a). The mutation spectrum of ACC was characterized by a predominance of C>T transitions, as found in most solid tumors39,40 (Fig. 2b); however, the two hypermutated samples displayed a remarkably high frequency of C>A transversions (Supplementary Fig. 3b). As reported in hepatocellular carcinomas32,41, the prevalence of T>C mutations was significantly higher on the transcribed strand than on the untranscribed strand (P = 6.2 x 10-7; Fig. 2b). This bias may be introduced by DNA adducts on adenine being less efficiently repaired by transcription-coupled nucleotide excision repair (NER) on the untranscribed strand.

We used MutSigCV40 to identify genes harboring more nonsynony- mous mutations than expected by chance given gene size, sequence context and the mutation rate of each tumor. Four genes were signifi- cantly enriched for mutations (q value < 0.05): CTNNB1, TP53, DAXX and MEN1 (Supplementary Table 5). Five other genes displayed homozygous deletions and/or damaging mutations (indels or non- sense or missense mutations predicted to be damaging by PolyPhen-2 software) in three or more tumors (Fig. 2c). These genes were further

Figure 2 Mutational landscape of ACC. (a) Number of mutations in exons per megabase of DNA in the 45 ACC cases. Clinical features associated with the number of mutations are indicated by color (5-year survival: black, deceased; gray, alive). Missing values are indicated by diagonal lines. (b) Relative proportions of the 6 possible base-pair substitutions on the transcribed and non- transcribed strands among the 1,118 point mutations identified in the exome sequences of 43 non-hypermutated ACCs (means with s.d .; binomial test). (c) Genes altered by damaging mutations and/or homozygous deletions in at least 3 of 45 ACCs. The number of damaging mutations (nonsense, frameshift or missense mutations predicted to be damaging by PolyPhen-2) is plotted on the x axis, and the number of homozygous deletions is plotted on the y axis. The size of each circle indicates the total number of damaging events, and the significance of mutation recurrence, calculated by MutSig software, is represented by color.

a

25

Synonymous

Number of mutations

20

Nonsynonymous

per Mb

15

10

1

0.

5

0

5-year survival

Weiss score

4

ENSAT stage

6

7 5

7

8

7

7

3

9

6

6

7

5

7

5

2

4

8

5 6

8

4 6

8

8

6

3

9

5

8

4 6

E

IL 6

3 3

2

3

9

3

7

3

2

8

2

4 3

4

2

2

3

4

3

1

4

2

4

2

2

2

4

2

4

4 . 4 + 3

2

4

2 2

4

2

4

2 2 2

2 2

2

2

2

4

2

2

1

1

4

b

Transcribed

C

8

ZNRF3

log10 (MutSig P value)

Percentage of mutations

0.3

Non-transcribed

Number of

homozygous deletions

-6

6

C

-1

0.2

P = 6.2 × 10

7

4

CDKN2A

2

RB1

TP53

Number of events

10

0.1

MCAM MEN1

0

MED12 DAXX

3

CTNNB1

0

C>T

C>A

C>G

T>C

T>A

T>G

0

2

4

6

8

10

Number of damaging mutations

npg

screened by targeted sequencing in the validation cohort of 77 tumors (Supplementary Table 6).

In the entire series of 122 ACCs (discovery cohort and valida- tion cohort), 9 genes (ZNRF3, CTNNB1, TP53, CDKN2A, RB1, MEN1, DAXX, MED12 and TERT) displayed damaging muta- tions, homozygous deletions or high-level amplifications in ≥5% of ACCs (Fig. 3a and Supplementary Table 7). These genes were altered at similar frequencies in the discovery and validation cohorts (Supplementary Table 8). ZNRF3 was the most frequently altered gene (21%). This gene encodes a cell surface transmembrane E3 ubiq- uitin ligase that acts as a negative feedback regulator of Wnt/B-catenin signaling by promoting the degradation of the LRP6 and Frizzled receptors6. We found homozygous deletions of ZNRF3 in 19 tumors (Fig. 1b) and mutations in 7 tumors: 4 nonsense alterations affecting sequence upstream of the transmembrane domain, presumably lead- ing to a truncated protein unable to anchor properly to the cell mem- brane, and 3 alterations located within or close to the RING domain (Supplementary Fig. 4). Five of the seven mutations were homozygous following loss of the second allele by LOH (Supplementary Table 7). CTNNB1 and ZNRF3 alterations were mutually exclusive in our series (P = 0.007, Fisher’s exact test). Tumors with altered ZNRF3 showed transcriptional activation of ß-catenin targets, but this activation was weaker than in CTNNB1-mutated tumors (Supplementary Fig. 5).

Mutations of APC were also identified in 2% of cases. Taken together, 47 tumors (39%) harbored an alteration affecting the Wnt/ß-catenin pathway (Fig. 3b). The p53-Rb pathway was the second most fre- quently altered pathway. Inactivating mutations or homozygous dele- tions of TP53 (16%), CDKN2A (11%) and RB1 (7%) and high-level amplifications of CDK4 (2%) and MDM2 (2%) altered this pathway in 40 cases (33%). TP53 mutations preferentially affected the DNA- binding domain or the oligomerization domain (Supplementary Fig. 4). Pathway enrichment analysis also identified chromatin remodeling as a recurrently altered function in ACC. In addition to mutations in MEN1 (7%) and DAXX (6%), 2 tumors had ATRX alter- ations, and 13 genes associated with chromatin remodeling42 were mutated in a single case (Supplementary Table 9). ATRX and DAXX mutations have been associated with the alternative lengthening of telomeres (ALT) phenotype43. Consistent with this observation, we found the ALT phenotype in 4 of 5 ACCs with ATRX or DAXX muta- tions and in only 1 of the 15 other ACCs tested (P = 0.0049, Fisher’s exact test; Supplementary Table 10). The MED12 gene (encoding mediator complex subunit 12) harbored damaging mutations in six tumors. Recurrent MED12 mutations were previously reported in leiomyomas44 and prostate cancer45. In ACC, MED12 mutations affected the C-terminal region of the protein (Supplementary Fig. 4), which was found to physically interact with ß-catenin, leading to the recruitment of mediator to Wnt-responsive genes46. MED12 mutations may there- Percentage of tumors with an event 0 5 10 15 20 L fore perturb the transcription of ß-catenin target genes.

Figure 3 Overview of somatic alterations in genes frequently altered in ACC. (a) Mutations and focal copy number alterations are depicted with color for the 9 genes with an alteration frequency above 5% (red line) in the complete series of 122 ACCs (discovery and validation samples). (b) Major pathways commonly altered in ACC by mutation, high-level amplification or homozygous deletion. Frequencies of alterations are expressed as percentages of all cases, with backgrounds in red for activated genes and in blue for inactivated genes. Darker colors reflect more frequent alterations.

a

Missense Splice site

Nonsense

Homozygous deletion

Indel

Focal amplification

ß-catenin pathway

ZNRF3

21%

CTNNB1

16%

TP53

16%

p53/Rb signaling

CDKN2A

11%

RB1

7%

Chromatin remodeling

MEN1

7%

DAXX

6%

MED12

5%

TERT

6%

b

ß-catenin pathway

p53/Rb signaling

Chromatin remodeling

ZNRF3 21%

APC 2%

CDK4 2%

RB1 7%

MEN1 7%

CTNNB1 16%

CDKN2A 11%

MDM2 2%

TP53

DAXX 6%

ATRX

16%

4%

Figure 4 miRNA expression patterns in ACC. (a) Unsupervised classification of 45 ACCs identifies 3 stable miRNA expression clusters. The consensus matrix represents the similarity between samples. Consensus index values range from O (highly dissimilar profiles; white) to 1 (highly similar profiles; dark blue). Samples are ordered on the x and y axes according to consensus clustering, which is depicted above the heat map. (b) Heat map of miRNA expression profiles. The expression levels of the miRNAs with the most variant expression are shown by color. Samples are ordered by miRNA expression subgroup. The expression levels of the same miRNAs in normal adrenal samples are represented on the right. (c) Box-and-whisker plots show the distribution of miR-506-514 expression levels relative to each sample group (NT, non-tumoral adrenal; ACA, adrenocortical adenoma; ACCs are classified by miRNA expression subgroup). Middle bar, median; box, interquartile range; bars extend to 1.5 times the interquartile range. (d) Deregulation of the DLK1-MEG3 cluster in Mil ACCs. The mean expression levels of miRNAs belonging to the DLK1-MEG3 cluster are indicated, together with the presence of an LOH at the locus (gray, LOH; white, no LOH) and the methylation levels at the promoter CpG island (yellow, fully methylated; green, hemimethylated; dark green, unmethylated). Samples are ordered as in b.

To obtain a comprehensive view of the genomic landscape of ACC, we analyzed DNA methylation, mRNA expression and miRNA expres- sion in tumors from the discovery set (Supplementary Table 11). Using the recursively partitioned mixture model (RPMM)47, we identified four DNA methylation-based tumor clusters. Two correspond to the previously described ‘CIMP-high’ and ‘CIMP-low’ subgroups23 (CIMP, CpG island methylator phenotype), associated with poor prognosis. The ‘non-CIMP’ tumors fell into two clusters, one of which showed widespread hypomethylation of CpG sites located outside CpG islands (Supplementary Fig. 6). Consensus clustering of mRNA expression profiles confirmed the existence of two main transcriptional clusters strongly correlated with survival-the aggressive C1A and indolent C1B clusters16. C1A tumors could be classified into three subgroups strongly associated with DNA methylation clusters (Supplementary Fig. 7).

npg

We assayed miRNA expression levels in 45 ACCs of the discovery cohort and 3 normal adrenal samples by Illumina sequencing. As previously reported26,28,31, MIR483 (encoding hsa-miR-483), which is located within intron 2 of the IGF2 locus, was overexpressed in ACC (Supplementary Table 12). Consensus clustering identified three stable tumor clusters (Fig. 4a). Cluster Mil showed the largest miRNA expression differences relative to normal adrenal samples (Supplementary Fig. 8 and Supplementary Table 12). Mil tumors

a

b

Mi1

Mi2

Mi3

Mi1

Mi2

Mi3

MIR-506-514

DLK1-MEG3

miRNA expression

Normal tissue

C

Normalized expression of MIR-506-514

0.7

Low

High

0.6

d

0.5

miRNA cluster

0.4

expression

0.3

LOH 14q32

0.2

DNA

0.1

methylation

0

·

NT

ACA

Mi3

Mi1

Mi2

Mi3

Mi1

Mi2

ACC

were characterized by the upregulation of 11 miRNAs belonging to the miRNA-506-514 cluster (Xq27.3) and by the downregulation of 38 miRNAs belonging to the imprinted DLK1-MEG3 cluster (14q32.2) (Fig. 4b). The miRNA-506-514 cluster, which was also weakly overexpressed in the Mi2 subgroup (Fig. 4c), has an oncogenic role in melanoma48. The DLK1-MEG3 cluster includes 2 long noncoding RNAs (MEG3 and MEG8) and 54 miRNAs expressed from the maternally inherited homolog49. SNP array analysis identi- fied LOH of chromosome arm 14q in all Mil tumors (compared to 1 of 23 non-Mil samples; P = 1.1 x 10-6, Fisher’s exact test), associated with a shift in MEG3 promoter methylation from hemimethylation to full methylation (P = 5.1 x 10-7, Wilcoxon rank-sum test; Fig. 4d). Thus, loss of the maternal unmethylated allele results in silencing of the DLK1-MEG3 miRNA cluster in Mil tumors. The DLK1- MEG3 cluster is also silenced in pituitary adenomas and ovarian cancers50,51 and appears to have tumor-suppressive activity52,53. Interestingly, MEG3 is highly expressed in normal adrenal cortex, pituitary gland and ovary compared to other human tissues (Supplementary Fig. 9), suggesting that this cluster has a key role in endocrine tissues.

Our integrative study identified substantial overlap between the different omics classifications (Fig. 5 and Supplementary Table 11).

Figure 5 Coordinated analysis of genomic features in ACC subtypes. (a) mRNA expression, DNA methylation and miRNA expression clusters delineate two main molecular groups: C1A and C1B. C1A ACCs can be classified into three subgroups according to methylation profile: CIMP-high, CIMP-low and non-CIMP. C1B ACCs can be classified into two subgroups according to miRNA expression profile (Mil or Mi2). miRNA deregulation, mutation rate and alterations in key genes and pathways are shown for each tumor of the discovery set analyzed on all platforms. (b) Overall survival rate in the various molecular subgroups of ACC defined by omics integration.

a

mRNA subgroup

b

C1A

C1B

C1A

C1B

Omics subgroups

mRNA

DNA methylation

DNA methylation subgroup

1.0

CIMP-high

C1B Mi1

miRNA

CIMP-low

Deregulated miRNAs

miR-506-514

Non-CIMP

Overall survival rate

0.8

C1B Mi2

DLK1-MEG3

miRNA cluster

Mi1

0.6

ZNRF3

Mi2

Driver genes

CTNNB1

Mi3

0.4

C1A non-CIMP

TP53

Drivers and pathways

Altered

Wnt/ß-catenin

Not altered

0.2

C1A CIMP-low

Pathways

p53/Rb

Deregulated miRNAs

Log-rank

Chromatin

0

P = 3.67 x 10

C1A CIMP-high

remodeling

Down

Up

Mutation rate

Mutation rate (per Mb)

0

24

48

72

96

120

144

I

Time (months)

0.05

1.4

npg

Transcriptome clusters were strongly correlated with subgroups based on DNA methylation (P = 6.3 x 10-5, x2 test) and miRNA expression (P = 3.4 x 10-5, x2 test). Indeed, the C1A subgroup with poor prog- nosis included almost all CIMP and Mi3 tumors. In contrast, C1B tumors with good prognosis were generally non-CIMP and belonged to the Mi1 or Mi2 miRNA cluster. The mutation rate was significantly higher in the C1A group (mean = 0.75 mutations per megabase) than in the C1B group (mean = 0.32 mutations per megabase; P = 7.5 x 10-4, Wilcoxon rank-sum test), and the key genes and pathways iden- tified by exome and SNP analysis were altered primarily in the C1A group. In contrast, miRNA deregulation was mainly observed in the C1B group.

In conclusion, we have identified new genes and pathways that are recurrently altered in ACC and have refined the molecular clas- sification of the disease. Our data strongly suggest that ZNRF3 is a new tumor suppressor gene related to the ß-catenin pathway. Further functional studies are needed to investigate the role of long noncoding RNAs and miRNAs in ACC.

URLs. European Network for the Study of Adrenal Tumors (ENSAT), http://www.ensat.org/; dbSNP, http://www.ncbi.nlm.nih.gov/SNP/; 1000 Genomes Project, http://www.1000genomes.org/; Illumina CASAVA software, http://support.illumina.com/sequencing/ sequencing_software/casava/downloads.ilmn; Geneious software, http://www.geneious.com/; Integrative Genomics Viewer, http://www. broadinstitute.org/igv/; Bioconductor, http://www.bioconductor. org/; Genome Alteration Print (GAP), http://bioinfo-out.curie.fr/ projects/snp_gap/; Gene Expression Omnibus (GEO), http://www. ncbi.nlm.nih.gov/geo/; European Genome-phenome Archive (EGA), http://www.ebi.ac.uk/ega/.

METHODS

Methods and any associated references are available in the online version of the paper.

Accession codes. SNP array, DNA methylation, mRNA and miRNA expression data have been deposited in the Gene Expression Omnibus (GEO) under accession GSE49280. Exome sequencing data have been deposited in the European Genome-phenome Archive (EGA), which is hosted by the European Bioinformatics Institute (EBI), under accession EGAS00001000665.

Note: Any Supplementary Information and Source Data files are available in the online version of the paper.

ACKNOWLEDGMENTS

We thank F. Letourneur, R. Pelletier and S. Jacques from the Genomic platform of the Cochin Institute and P. Nietscke from the Paris Descartes University Bioinformatics platform for their technical support, J. Métral and J. Godet from the Ligue Nationale Contre le Cancer for organization of the Cartes d’Identité des Tumeurs (CIT) program, P.F. Plouin for organization of the COMETE network, the tumor bank of Cochin Hospital (B. Terris and M. Sibony) for help in sample collection, the Oncogenetic Department of Cochin Hospital (V. Duchossoy), A. Steel for management of the ENSAT database, the members of our laboratories and the COMETE and ENSAT networks for support and discussions, and all the staffs of the clinical and pathology departments who were involved in patient care. This study is part of the CIT Program from La Ligue Nationale Contre le Cancer. It was supported by funding from the Programme Hospitalier de Recherche Clinique to the COMETE network (grant AOM95201), the Seventh Framework Programme (FP7/2007-2013) under grant agreement 259735, Institut National du Cancer Recherche Translationelle 2009-RT-02, the Institut National du Cancer (to the Rare Adrenal Cancer Network COMETE), INSERM (G.A. receives a Contrat d’Interface) and the Conny-Maeva Charitable Foundation (to the laboratory of J.B.).

AUTHOR CONTRIBUTIONS

G.A., E.L., A.d.R. and J.B. conceived the study, designed the experiments, analyzed the data and wrote the manuscript. E.L., A.d.R., G.A. and W.L. performed the

bioinformatics analyses. G.A., E.L., A.d.R., B.R. and J.B. provided analytical advice. A.d.R. and J.B. supervised the study. R.L., A.J., G.A., E.L., W.L. and N.E. managed the data. F.T. reviewed the histopathology. A.J. performed the Sanger sequencing experiments. K.P., F.R .- C., S.R., G.A., A.J., O.B., S.S., H.O. and B.R. performed the molecular analyses. E.C. performed the microsatellite instability experiments. M.F. coordinated the ENSAT centers. M.F., J.B., G.A., M.K., B.A., J.W., M.Q., M.M., F.M., T.P., R.D.K., A.T., V.K., E.B., B.D., L.G., L.A., X.B., F.B. and R.L. recruited the subjects. All authors discussed the results and implications and commented on the manuscript.

COMPETING FINANCIAL INTERESTS

The authors declare no competing financial interests.

Reprints and permissions information is available online at http://www.nature.com/ reprints/index.html.

1. Schteingart, D.E. et al. Management of patients with adrenal cancer: recommendations of an international consensus conference. Endocr. Relat. Cancer 12, 667-680 (2005).

2. Abiven, G. et al. Clinical and biological features in the prognosis of adrenocortical cancer: poor outcome of cortisol-secreting tumors in a series of 202 consecutive patients. J. Clin. Endocrinol. Metab. 91, 2650-2655 (2006).

3. Fassnacht, M. et al. Limited prognostic value of the 2004 International Union Against Cancer staging classification for adrenocortical carcinoma: proposal for a Revised TNM Classification. Cancer 115, 243-250 (2009).

4. Wandoloski, M., Bussey, K.J. & Demeure, M.J. Adrenocortical cancer. Surg. Clin. North Am. 89, 1255-1267 (2009).

5. Soon, P.S.H. & Sidhu, S.B. Molecular basis of adrenocortical carcinomas. Minerva Endocrinol. 34, 137-147 (2009).

6. Hao, H .- X. et al. ZNRF3 promotes Wnt receptor turnover in an R-spondin-sensitive manner. Nature 485, 195-200 (2012).

7. Grumbach, M.M. et al. Management of the clinically inapparent adrenal mass (‘incidentaloma’). Ann. Intern. Med. 138, 424-429 (2003).

8. Kebebew, E., Reiff, E., Duh, Q .- Y., Clark, O.H. & McMillan, A. Extent of disease at presentation and outcome for adrenocortical carcinoma: have we made progress? World J. Surg. 30, 872-878 (2006).

9. Kerkhofs, T.M.A. et al. Adrenocortical carcinoma: a population-based study on incidence and survival in the Netherlands since 1993. Eur. J. Cancer 49, 2579-2586 (2013).

10. Luton, J .- P. et al. Clinical features of adrenocortical carcinoma, prognostic factors, and the effect of mitotane therapy. N. Engl. J. Med. 322, 1195-1201 (1990).

11. Icard, P. et al. Adrenocortical carcinomas: surgical trends and results of a 253-patient series from the French Association of Endocrine Surgeons study group. World J. Surg. 25, 891-897 (2001).

12. Allolio, B. & Fassnacht, M. Clinical review: adrenocortical carcinoma: clinical update. J. Clin. Endocrinol. Metab. 91, 2027-2037 (2006).

13. Fassnacht, M. et al. Combination chemotherapy in advanced adrenocortical carcinoma. N. Engl. J. Med. 366, 2189-2197 (2012).

14. Libè, R., Fratticci, A. & Bertherat, J. Adrenocortical cancer: pathophysiology and clinical management. Endocr. Relat. Cancer 14, 13-28 (2007).

15. Tissier, F. et al. Mutations of ß-catenin in adrenocortical tumors: activation of the Wnt signaling pathway is a frequent event in both benign and malignant adrenocortical tumors. Cancer Res. 65, 7622-7627 (2005).

16. de Reyniès, A. et al. Gene expression profiling reveals a new classification of adrenocortical tumors and identifies molecular predictors of malignancy and survival. J. Clin. Oncol. 27, 1108-1115 (2009).

17. Giordano, T.J. et al. Molecular classification and prognostication of adrenocortical tumors by transcriptome profiling. Clin. Cancer Res. 15, 668-676 (2009).

18. Laurell, C. et al. Transcriptional profiling enables molecular classification of adrenocortical tumours. Eur. J. Endocrinol. 161, 141-152 (2009).

19. Barreau, O. et al. Clinical and pathophysiological implications of chromosomal alterations in adrenocortical tumors: an integrated genomic approach. J. Clin. Endocrinol. Metab. 97, E301-E311 (2012).

20. Stephan, E.A. et al. Adrenocortical carcinoma survival rates correlated to genomic copy number variants. Mol. Cancer Ther. 7, 425-431 (2008).

21. De Martino, M.C. et al. Molecular screening for a personalized treatment approach in advanced adrenocortical cancer. J. Clin. Endocrinol. Metab. 98, 4080-4088 (2013).

22. Ronchi, C.L. et al. Single nucleotide polymorphism array profiling of adrenocortical tumors-evidence for an adenoma carcinoma sequence? PLOS ONE 8, e73959 (2013).

23. Barreau, O. et al. Identification of a CpG island methylator phenotype in adrenocortical carcinomas. J. Clin. Endocrinol. Metab. 98, E174-E184 (2013).

24. Rechache, N.S. et al. DNA methylation profiling identifies global methylation differences and markers of adrenocortical tumors. J. Clin. Endocrinol. Metab. 97, E1004-E1013 (2012).

25. Fonseca, A.L. et al. Comprehensive DNA methylation analysis of benign and malignant adrenocortical tumors. Genes Chromosom. Cancer 51, 949-960 (2012).

26. Soon, P.S.H. et al. miR-195 and miR-483-5p identified as predictors of poor prognosis in adrenocortical cancer. Clin. Cancer Res. 15, 7684-7692 (2009).

27. Tömböl, Z. et al. Integrative molecular bioinformatics study of human adrenocortical tumors: microRNA, tissue-specific target prediction, and pathway analysis. Endocr. Relat. Cancer 16, 895-906 (2009).

28. Patterson, E.E., Holloway, A.K., Weng, J., Fojo, T. & Kebebew, E. MicroRNA profiling of adrenocortical tumors reveals miR-483 as a marker of malignancy. Cancer 117, 1630-1639 (2011).

29. Schmitz, K.J. et al. Differential expression of microRNA-675, microRNA-139-3p and microRNA-335 in benign and malignant adrenocortical tumours. J. Clin. Pathol. 64, 529-535 (2011).

30. Özata, D.M. et al. The role of microRNA deregulation in the pathogenesis of adrenocortical carcinoma. Endocr. Relat. Cancer 18, 643-655 (2011).

31. Chabre, O. et al. Serum miR-483-5p and miR-195 are predictive of recurrence risk in adrenocortical cancer patients. Endocr. Relat. Cancer 20, 579-594 (2013).

32. Guichard, C. et al. Integrated analysis of somatic mutations and focal copy-number changes identifies key genes and pathways in hepatocellular carcinoma. Nat. Genet. 44, 694-698 (2012).

33. Gicquel, C. et al. Rearrangements at the 11p15 locus and overexpression of insulin- like growth factor-II gene in sporadic adrenocortical tumors. J. Clin. Endocrinol. Metab. 78, 1444-1453 (1994).

34. Pasic, I. et al. Recurrent focal copy-number changes and loss of heterozygosity implicate two noncoding RNAs and one tumor suppressor gene at chromosome 3q13.31 in osteosarcoma. Cancer Res. 70, 160-171 (2010).

35. Liu, Q. et al. LncRNA loc285194 is a p53-regulated tumor suppressor. Nucleic Acids Res. 41, 4976-4987 (2013).

36. Zack, T.I. et al. Pan-cancer patterns of somatic copy number alteration. Nat. Genet. 45, 1134-1140 (2013).

37. Du, Z. et al. Integrative genomic analyses reveal clinically relevant long noncoding RNAs in human cancer. Nat. Struct. Mol. Biol. 20, 908-913 (2013).

38. Cheetham, S.W., Gruhl, F., Mattick, J.S. & Dinger, M.E. Long noncoding RNAs and the genetics of cancer. Br. J. Cancer 108, 2419-2425 (2013).

39. Greenman, C. et al. Patterns of somatic mutation in human cancer genomes. Nature 446, 153-158 (2007).

40. Lawrence, M.S. et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature 499, 214-218 (2013).

41. Alexandrov, L.B. et al. Signatures of mutational processes in human cancer. Nature 500, 415-421 (2013).

42. Gonzalez-Perez, A., Jene-Sanz, A. & Lopez-Bigas, N. The mutational landscape of chromatin regulatory factors across 4623 tumor samples. Genome Biol. 14, r106 (2013).

43. Heaphy, C.M. et al. Altered telomeres in tumors with ATRX and DAXX mutations. Science 333, 425 (2011).

44. Mäkinen, N. et al. MED12, the mediator complex subunit 12 gene, is mutated at high frequency in uterine leiomyomas. Science 334, 252-255 (2011).

45. Barbieri, C.E. et al. Exome sequencing identifies recurrent SPOP, FOXA1 and MED12 mutations in prostate cancer. Nat. Genet. 44, 685-689 (2012).

46. Kim, S., Xu, X., Hecht, A. & Boyer, T.G. Mediator is a transducer of Wnt/B-catenin signaling. J. Biol. Chem. 281, 14066-14075 (2006).

47. Houseman, E.A. et al. Model-based clustering of DNA methylation array data: a recursive-partitioning algorithm for high-dimensional data arising as a mixture of beta distributions. BMC Bioinformatics 9, 365 (2008).

48. Streicher, K.L. et al. A novel oncogenic role for the miRNA-506-514 cluster in initiating melanocyte transformation and promoting melanoma growth. Oncogene 31, 1558-1570 (2012).

49. Benetatos, L. et al. The microRNAs within the DLK1-DIO3 genomic region: involvement in disease pathogenesis. Cell. Mol. Life Sci. 70, 795-814 (2013).

50. Cheunsuchon, P. et al. Silencing of the imprinted DLK1-MEG3 locus in human clinically nonfunctioning pituitary adenomas. Am. J. Pathol. 179, 2120-2130 (2011).

51. Zhang, L. et al. Genomic and epigenetic alterations deregulate microRNA expression in human epithelial ovarian cancer. Proc. Natl. Acad. Sci. USA 105, 7004-7009 (2008).

52. Haga, C.L. & Phinney, D.G. MicroRNAs in the imprinted DLK1-DIO3 region repress the epithelial-to-mesenchymal transition by targeting the TWIST1 protein signaling network. J. Biol. Chem. 287, 42695-42707 (2012).

53. Zhou, Y. et al. Activation of p53 by MEG3 non-coding RNA. J. Biol. Chem. 282, 24731-24742 (2007).

odu

ONLINE METHODS

Tumor samples. We included 130 ACCs belonging to 2 different cohorts in this study. Tumors from a discovery cohort (core set of 45 cases, 53 cases ana- lyzed on at least 1 omics platform) were used for integrated genomic analyses: exome sequencing (45 tumor-normal pairs), SNP array (45 tumors), DNA methylation analysis (51 tumors), mRNA expression array (47 tumors) and miRNA sequencing (45 tumors). Tumors from a validation cohort of 77 cases were analyzed by targeted sequencing and SNP array.

The 53 ACCs of the discovery cohort were collected by the Cochin Hospital team from patients prospectively included in the COMETE network54 (Supplementary Table 1). Each tumor was snap frozen immediately after surgery, and the tumor was dissected by the pathologist and then kept in liquid nitrogen until use. Diagnosis of ACC, including analysis of sections on an adjacent block, was performed by an expert pathologist for each tumor. Non-tumoral adrenal tissues and adrenocortical adenomas were collected, diagnosed by the pathologist and stored similarly in the COMETE tumor bank. Minimal endocrine workup and hormone assays were performed as previously reported16. Metastasis was diagnosed by systematic imaging inves- tigations, principally abdominal and chest computed tomography (CT) scans and, when appropriate, bone scintigraphy and magnetic resonance imaging (MRI). Tumor staging was performed using ENSAT staging criteria3. For each case, malignancy diagnosis and tumor weight, size and classification were determined by pathological examination. Malignancy was assessed according to Weiss criteria. After surgery, patients were examined at least twice a year for 2 years and annually thereafter. Patients were followed until the date of their death, their last examination or the end of the follow-up period. Signed informed consent for molecular analysis of tumor and leukocyte DNA and for access to data collected was obtained from all patients, and the study was approved by the Institutional Review Board of Cochin Hospital.

The 77 ACCs of the validation cohort were collected by 16 European centers in the ENSAT network (Supplementary Table 1). For all cases, malignancy was confirmed on the basis of pathological examination by expert pathologists. Minimal clinical information from the ENSAT database55 included age, sex, tumor side and size, existence of hormonal secretion, ENSAT tumor stage, disease-free survival and overall survival. Signed informed consent for analysis of tumor and leukocyte DNA and for access to data collected was obtained from all patients, and the study was approved by the local institutional review board of each clinical center.

DNA and RNA extraction. Tumor samples (10-50 mg) were powdered under liquid nitrogen. DNA was extracted and purified by cesium chloride gradi- ent ultracentrifugation or proteinase K digestion and ethanol extraction, and samples were cleaned using columns (Qiagen). DNA concentrations were determined using a Nanodrop ND-1000 instrument (Nyxor Biotech).

RNA was extracted using RNAble (Eurobio) and cleaned by RNeasy column (Qiagen). Aliquots of RNA were analyzed by electrophoresis on a Bioanalyzer 2100 (version A.02 S1292, Agilent Technologies) and quantified using a Nanodrop ND-1000. Stringent criteria for RNA quality were applied to rule out degradation, especially a 28S/18S ratio above 1.8. A similar protocol was applied for miRNA extraction, using miRNeasy columns (Qiagen).

SNP array analysis. Samples from 121 ACCs were analyzed with Illumina HumanCNV610-Quad v1.0 (n = 7), HumanOmniExpress v12 (n = 38) or HumanCore-12v1 (n = 76) BeadChips. Hybridization was performed by Integragen, according to the instructions provided by the array manufacturer. Raw fluorescent signals were imported into BeadStudio software (Illumina) and normalized as previously described56 to obtain log R ratio (LRR) and B allele frequency (BAF) values. The tQN normalization procedure was then applied to correct for asymmetry in BAF signals due to bias between the two dyes used in Illumina assays57.

Genomic profiles were divided into homogeneous segments by applying the circular binary segmentation algorithm58 to both LRR and BAF values. We then used the Genome Alteration Print (GAP) method59 to determine the ploidy of each sample, the level of contamination with normal cells and the allele-specific copy number of each segment. Chromosome aberrations were defined using empirically determined thresholds32 as follows: gain, copy number ≥ ploidy + 1; loss, copy number ≤ ploidy - 1; high-level amplification,

copy number > ploidy + 2; homozygous deletion, copy number = 0. Finally, we considered a segment to have undergone LOH when the copy number of the minor allele was equal to 0. Lists of homozygous deletions and focal amplifica- tions, defined by at least five consecutive probes, were generated and verified manually to remove doubtful events. Significantly recurrent copy number changes were identified using the GISTIC2.0 algorithm60. The overall genomic instability of each sample was quantified as the FAA score, i.e., the proportion of chromosome arms altered on more than 40% of their length.

Whole-exome sequencing. Exome sequencing was performed in 45 pairs of ACCs and matched normal samples (from blood or adipose tissue). In all tumor samples, the ratio of tumor cells to non-tumor cells was determined to be >50%, and this value was confirmed by SNP array analysis. Sequence capture and exome sequencing were performed by Integragen as previously described32. In brief, Agilent in-solution enrichment (SureSelect Human All Exon Kit v4+UTR) with the provided biotinylated oligonucleotide probe library (Human All Exon v4+UTR-70 Mb) was used for DNA sequence cap- ture, enrichment and elution from the tumor and matched blood DNA sam- ples. The eluted fraction was amplified by PCR and sequenced on an Illumina HiSeq 2000 sequencer as paired-end 75-bp reads61. Image analysis and base calling were performed using Illumina Real-Time Analysis (RTA) Pipeline version 1.12.4.2 with default parameters.

Analysis pipeline and mutation annotation. Illumina CASAVA 1.8.2 soft- ware was used to align reads against the hg19 genome build (GRCh37) and to detect single-nucleotide variants and indels. The alignment algorithm used was ELANDv2. Targeted regions were sequenced to an average depth of 90.3x, with 99% of the regions covered by ≥1x, 93.0% covered by ≥10x and 84% covered by ≥25x (Supplementary Fig. 2). We applied a previously described pipeline32, with some modifications, to generate a list of somatic variants located in coding regions plus two intronic bases corresponding to consensus splice sites. Quality control filtering removed variants sequenced in <10 reads, with <3 variant calls or with QPHRED of <10. We used relaxed thresholds to detect somatic variants in order to minimize the risk of filter- ing out true somatic mutations. Variants were considered to be of somatic origin when the frequency of variant reads was ≥10% in the tumor and <5% in the normal counterpart, with a significant enrichment of variant calls in the tumor as assessed through a Fisher’s exact test. Common polymorphisms (with reported frequency of >1%) were removed by comparison with dbSNP135, the 1000 Genomes Project database and a proprietary database of exomes from normal tissues. Missense variants were classified as ‘benign’, ‘possibly damaging’ or ‘probably damaging’ using PolyPhen-2 software. All nonsynony- mous somatic mutations, except those encountered in the two hypermutated samples, were validated visually using the Integrative Genomics Viewer62. The mean number of somatic mutations per sample (6 silent and 24 non- silent) comprised exonic single-nucleotide variants and indels. The muta- tion rate in each tumor was evaluated by dividing the number of somatic mutations by the number of exonic bases covered by ≥10x in both the tumor and normal samples. To identify mutations with transcriptional strand bias, we used a binomial test to assess whether the proportions of each type of somatic substitution occurring on the transcribed strand differed significantly from 0.5. For all variants, validation status determined by visual examination using the Integrative Genomics Viewer and Sanger sequencing is provided (Supplementary Table 4).

Identification of significantly mutated genes and pathways. We used MutSigCV to identify genes harboring more nonsynonymous mutations than expected by chance given gene size, sequence context and the mutation rate of each tumor40. We used as genomic covariates the mean expression level of each gene in our ACC expression data set, the DNA replication time63 and the HiC statistic of chromatin state64 available in MutSig reference files. We also searched for gene sets harboring more damaging mutations than expected by chance. Given set G of all the genes sequenced with sufficient coverage, the set S of tumor samples (of size n) and any gene set P, we calculated the prob- ability of observing a number of mutations equal to or greater than the number observed in P across n samples according to a binomial law B(k, p), with k = n × L(P) and the mutation rate p = A(G, S)/(n × L(G)), where L(X) is the

npg

npg

sum of the lengths (in base pairs) of all genes or exons from a gene set X and A(G, S) is the total number of mutations observed in all targeted sequences across all samples from S.

Targeted sequencing of candidate drivers. Nine genes (CTNNB1, ZNRF3, TP53, CDKN2A, RB1, DAXX, MEN1, MED12 and MCAM) displayed ≥3 damaging events (indels or nonsense or missense mutations predicted to be damaging by PolyPhen-2 software) in the 45 ACCs of the discovery cohort analyzed by SNP array and exome sequencing. These genes, together with the known tumor suppressor APC, were sequenced in the validation cohort of 77 ACCs and matched blood samples. Multiplex PCR probes were designed using AmpliSeq v2.0 software (Life Technologies). Starting from 20 ng of genomic DNA, the selected genes were captured by PCR and sequenced using a Personal Genome Machine (Life Technologies) at the Genomics Platform of the Cochin Institute, using Ion 314 chips for normal DNA and Ion 316 chips for tumor DNA (Life Technologies). Read alignment against the hg19 refer- ence genome and variant calling were performed using Torrent Suite Software v3.6.63335 (Life Technologies). Sequence variants were annotated using ANNOVAR65 v521 and filtered using similar criteria as those used for exome sequencing data.

Validation of selected mutations by Sanger sequencing. Genomic DNA from samples in both the exome and mutation validation cohorts was subjected to Sanger dye-terminator sequencing. Reference sequences were downloaded from the Ensembl Genome Browser. Primer design was performed using Primer-BLAST66. PCR amplifications were performed with AmpliTaq Gold (Applied Biosystems, Life Technologies). PCR products were purified and sequenced at the Genomics Platform of the Cochin Institute, using an ABI 3730 XL (Applied Biosystems, Life Technologies). Sequence analysis was carried out using Geneious software. The somatic nature of mutations was confirmed by sequencing tumors and their paired normal tissue. Overall, we validated 125 mutations identified by exome sequencing (out of 132 vari- ants tested, validation rate = 94.7%; Supplementary Table 4) and 42 variants identified by targeted sequencing (Supplementary Table 6). All mutations in candidate driver genes were validated.

Alternative lengthening of telomeres. The ALT phenotype was evaluated using the ALT-specific C-circle assay, as previously described67. In brief, ALT-specific C-circles (partially single-stranded DNA molecules specifically associated with ALT activity) were detected by isothermic amplification of C-circle complementary strand and hybridization with the sequence-specific 32P-labeled probe 32P-(CCCTAA)3. Experiments were performed by Capital Biosciences, starting from 200 to 500 ng of tumor DNA.

DNA methylation profiling. Whole-genome DNA methylation was analyzed in 51 ACCs of the discovery cohort and 30 adrenocortical adenomas using the Illumina Infinium HumanMethylation27 assay56, as previously described23. The ß-value DNA methylation scores for each locus were extracted together with detection P values from Illumina GenomeStudio software. As described elsewhere68, we replaced data points with detection P value > 0.05 with “NA” values, and we masked data points as “NA” for 4,484 probes that contained SNPs or overlapped with a repetitive element that was not uniquely aligned to the human genome or regions of insertions and deletions in the human genome. We used RPMM, developed for beta-distributed DNA methylation measurements47, to identify homogeneous tumor subgroups with similar methylation profiles. RPMM-based clustering was performed as implemented in the Bioconductor RPMM package, with default settings, using the 10% most variant CpG sites.

mRNA expression profiling. Gene expression was analyzed for 47 ACCs of the discovery cohort using Affymetrix Human Gene 2.0 ST arrays (n = 44) or HG-U133_Plus_2.0 arrays (n = 33)16, with 31 ACCs profiled with both chip types. All samples were normalized in batches independ- ently for each chip type using the RMA (Robust Multi-array Average) algorithm (Bioconductor affy package), and probe set intensities were then averaged for each gene symbol. Each of the two resulting matrices

was row centered using the 31 ACCs in common to calculate the center. The two matrices were then concatenated on the common gene symbols, and one ACC per duplicate was kept for further analyses. We used consensus clus- tering69 to identify homogeneous gene expression clusters. We selected the 20% most variant probes and established consensus partitions of the data set in K clusters (for K = 2, 3, … , 8), on the basis of 1,000 resampling iterations of hierarchical clustering, with Pearson’s dissimilarity as the distance metric and Ward’s method for linkage analysis. We used the cumulative distribution functions (CDFs) of the consensus matrices to determine the optimal number of clusters (K = 2 and K = 4), considering both the shape of the functions and the area under the CDF curves, as previously described69. The Bioconductor ConsensusClusterPlus package was used for consensus clustering analysis.

miRNA expression profiling. miRNA expression levels were assayed in 45 ACCs from the discovery cohort and 3 normal samples via Illumina sequenc- ing, as previously described70. In brief, small RNAs (<100 bases in length) were purified from total RNA samples using the miRNeasy kit (Qiagen), and multi- plexed miRNA libraries were prepared using a previously described protocol71 and sequenced on an Illumina HiSeq 2000 sequencer. Image analysis, base calling, demultiplexing and conversion of BCL to FASTQ format were carried out using Illumina CASAVA 1.8.2 software. Adaptor sequences were removed using mirExpress software72, and FASTA files for each sample were processed by miRanalyzer0.3 software73 to quantify read counts for each miRNA refer- enced in mirBase74 v18. In total, 765 miRNAs were expressed (>10 reads) in at least 2 samples and were used for unsupervised classification. miRNA count data were log2 transformed, divided by the total number of reads in each sample and centered on the mean expression level of each gene. Consensus clustering was then performed as described above, and the optimal number of clusters (K = 3) was determined from CDF curves. We used the Bioconductor DESeq package75 to identify differentially expressed miRNAs from count data.

Concordance between different omics platforms. We performed three tests to verify the proper pairing of samples analyzed by the different omics platforms. First, we verified that the genotypes obtained from exome sequencing and SNP array data were consistent. Then, we generated virtual copy number profiles from mRNA expression and DNA methylation array data, and we verified their consist- ency with the copy number profiles obtained from SNP arrays. Finally, we used the copy number and expression of genes located on sex chromosomes to predict the sex of each sample, and we verified agreement with clinical annotations.

54. Plouin, P .- F., Gimenez-Roqueplo, A .- P. & Bertagna, X. COMETE, a network for the study and management of hypersecreting adrenal tumors. Bull. Acad. Natl. Med. 192, 73-82 (2008).

55. Stell, A. & Sinnott, R. The ENSAT registry: a digital repository supporting adrenal cancer research. Stud. Health Technol. Inform. 178, 207-212 (2012).

56. Bibikova, M. et al. Genome-wide DNA methylation profiling using Infinium® assay. Epigenomics 1, 177-200 (2009).

57. Staaf, J. et al. Normalization of Illumina Infinium whole-genome SNP data improves copy number estimates and allelic intensity ratios. BMC Bioinformatics 9, 409 (2008).

58. Venkatraman, E.S. & Olshen, A.B. A faster circular binary segmentation algorithm for the analysis of array CGH data. Bioinformatics 23, 657-663 (2007).

59. Popova, T. et al. Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays. Genome Biol. 10, R128 (2009).

60. Mermel, C.H. et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 12, R41 (2011).

61. Gnirke, A. et al. Solution hybrid selection with ultra-long oligonucleotides for massively parallel targeted sequencing. Nat. Biotechnol. 27, 182-189 (2009).

62. Robinson, J.T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24-26 (2011).

63. Chen, C .- L. et al. Impact of replication timing on non-CpG and CpG substitution rates in mammalian genomes. Genome Res. 20, 447-457 (2010).

64. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 326, 289-293 (2009).

65. Wang, K., Li, M. & Hakonarson, H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 38, e164 (2010).

66. Ye, J. et al. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics 13, 134 (2012).

67. Henson, J.D. et al. DNA C-circles are specific and quantifiable markers of alternative-lengthening-of-telomeres activity. Nat. Biotechnol. 27, 1181-1185 (2009).

68. Hinoue, T. et al. Genome-scale analysis of aberrant DNA methylation in colorectal cancer. Genome Res. 22, 271-282 (2012).

69. Monti, S., Tamayo, P., Mesirov, J. & Golub, T. Consensus clustering: a resampling- based method for class discovery and visualization of gene expression microarray data. Mach. Learn. 52, 91-118 (2003).

70. Jung, A.C. et al. A poor prognosis subtype of HNSCC is consistently observed across methylome, transcriptome and miRNome analysis. Clin. Cancer Res. 19, 4174-4184 (2013).

71. Vigneault, F. et al. High-throughput multiplex sequencing of miRNA. Curr. Protoc. Hum. Genet. Chapter 11, Unit 11.12.1-10 (2012).

72. Wang, W .- C. et al. miRExpress: analyzing high-throughput sequencing data for profiling microRNA expression. BMC Bioinformatics 10, 328 (2009).

73. Hackenberg, M., Rodríguez-Ezpeleta, N. & Aransay, A.M. miRanalyzer: an update on the detection and analysis of microRNAs in high-throughput sequencing experiments. Nucleic Acids Res. 39, W132-W138 (2011).

74. Kozomara, A. & Griffiths-Jones, S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 39, D152-D157 (2011).

75. Anders, S. & Huber, W. Differential expression analysis for sequence count data. Genome Biol. 11, R106 (2010).