Whole-genome and transcriptome analysis of advanced adrenocortical cancer highlights multiple alterations affecting epigenome and DNA repair pathways.

Adrenocortical cancer (ACC) is a rare cancer of the adrenal gland. Several driver mutations have been identified in both primary and metastatic ACCs, but the therapeutic options are still limited. We performed whole-genome and transcriptome sequencing on seven patients with metastatic ACC. Integrative analysis of mutations, RNA expression changes, mutation signature, and homologous recombination deficiency (HRD) analysis was performed. Mutations affecting CTNNB1 and TP53 and frequent loss of heterozygosity (LOH) events were observed in our cohort. Alterations affecting genes involved in cell cycle (RB1, CDKN2A, CDKN2B), DNA repair pathways (MUTYH, BRCA2, ATM, RAD52, MLH1, MSH6), and telomere maintenance (TERF2 and TERT) consisting of somatic and germline mutations, structural variants, and expression outliers were also observed. HRDetect, which aggregates six HRD-associated mutation signatures, identified a subset of cases as HRD. Genomic alterations affecting genes involved in epigenetic regulation were also identified, including structural variants (SWI/SNF genes and histone methyltransferases), and copy gains and concurrent high expression of KDM5A, which may contribute to epigenomic deregulation. Findings from this study highlight HRD and epigenomic pathways as potential therapeutic targets and suggest a subgroup of patients may benefit from a diverse array of molecularly targeted therapies in ACC, a rare disease in urgent need of therapeutic strategies.

INTRODUCTION

Adrenocortical cancer (ACC) is a rare genitourinary malignancy, with an incidence of 1–2 cases per million patient-years; of these, ∼40% present with metastatic disease (Ng and Libertino 2003). Although most cases are localized, metastatic disease has a poor prognosis, with a 5-yr survival rate of ∼15% (Fassnacht et al. 2009). In addition, in more than one-half of cases, the tumors secrete steroid hormones, which can lead to Cushing's syndrome and medical complications thereof (Ng and Libertino 2003).

The current standard treatment of ACC is limited to cytotoxic chemotherapy (etoposide, doxorubicin, and cisplatin [EDP]) and the adrenolytic agent mitotane. The FIRM-ACT phase 3 clinical trial has established the combination of mitotane and triplet chemotherapy as a standard treatment for advanced ACC, with a median overall survival of 14.8 mo (Fassnacht et al. 2012). There is no established second-line therapy, and small phase 2 trials evaluating different targeted agents in all comers have shown no clear clinical benefit. The vascular endothelial growth factor (VEGF) antibody bevacizumab in combination with capecitabine (Wortmann et al. 2010), VEGF receptor inhibitors axitinib (O'Sullivan et al. 2014) or sunitinib (Kroiss et al. 2012), epidermal growth factor receptor (EGFR) inhibitor erlotinib in combination with gemcitabine (Quinkler et al. 2008), or the insulin-like growth factor type 1 receptor (IGF1R) antibody cixutumumab (Lerario et al. 2014) or inhibitor linsitinib (Fassnacht et al. 2015) have shown limited activity. Immune checkpoint inhibitors have also been investigated. A phase Ib expansion cohort from the JAVELIN solid tumor trial showed a response rate of 6% to the PD-L1 antibody avelumab (Le Tourneau et al. 2018). A phase II study of the PD-1 antibody pembrolizumab had an overall response rate of 23%; two of six patients with mismatch repair–deficient tumors (identified by immunohistochemistry [IHC]) had a response to therapy (Raj et al. 2020). In this context, better treatments are sorely needed.

Although studies of targeted agents have shown poor overall outcomes, molecular characterization of ACC may inform genetic changes involved in cancer progression and reveal potential therapeutic targets. Studies leveraging whole-exome sequencing and single-nucleotide polymorphism (SNP) and expression arrays in ACC have revealed recurrent mutations in CTNNB1, TP53, MEN1, PRKAR1A, RPL22, TERF2, TERT, CCNE1, RB1, and NF1 genes, frequent whole-genome doubling events, large-scale loss of heterozygosity (LOH), and increased ERBB4, IGF2, WNT, retinoic acid receptor, GPCR (Gq, G, Gi), and PDGFR signaling (Giordano et al. 2003; Tissier et al. 2005; Assié et al. 2014; Zheng et al. 2016; Gara et al. 2018; Fojo et al. 2020).

Here we describe the genomic profiles from whole-genome and transcriptome analysis (WGTA) of seven advanced ACC tumors, with detailed clinical and therapeutic annotations. These patients underwent molecular analysis as part of the personalized medicine program at BC Cancer (NCT02155621; Pleasance et al. 2020). In addition to identifying previously described recurrent alterations affecting genes involved in the WNT pathway, DNA repair, and cell cycle regulation, our analysis revealed frequent alterations in genes involved in epigenetic regulation. Specifically, we identified multiple copy gains affecting KDM5A, a lysine demethylase that antagonizes the histone methyltransferase activity associated with the tumor suppressor MEN1 (Feng et al. 2017), a well-described tumor suppressor in ACC. Multiple potentially actionable therapeutic targets were revealed from profiling metastatic ACC tumors, and treatment and response to systemic therapies are reported.

RESULTS

Clinical Presentations and Treatment Outcomes

Baseline characteristics at the time of diagnosis of metastatic ACC are summarized in Table 1. Overall, six of seven patients were female. The median age was 40 (range, 19–54), and five of seven patients had symptoms of cortisol excess, with one patient also having symptoms of aldosterone excess and another also having symptoms of testosterone excess. The clinical course is outlined in Figure 1. Standard chemotherapy and mitotane were given to four of seven patients, with an additional patient receiving this treatment in the adjuvant setting (Case 3). There were no responses to treatment except for Case 5 (response to single-agent mitotane, which had also been used in the adjuvant setting) and Case 6 (response to EDP chemotherapy without mitotane). Case 5 also showed stable disease for 8.5 mo on olaparib. There were no adverse events leading to treatment discontinuation, except for Case 6 in whom an immune checkpoint inhibitor led to autoimmune hepatitis requiring treatment discontinuation and corticosteroids.

Table 1.

Presenting clinical and histopathological characteristics of seven patients with metastatic adrenocortical carcinoma

Case number	1	2	3	4	5	6	7
Age, gender	40, M	19, F	24, F	54, F	46, F	34, F	48, F
Elevated cortisol	No	Yes	Yes	No	Yes	Yes	Yes
Other hormone elevation	No	No	No	No	Aldosterone	Testosterone	No
Source of tissue for WGTA	Primary	Primary	Soft tissue metastasis	Lymph node	Lung	Primary	Primary
Histopathology score
>5 mitotic figures/50 HPF	1	1	1	1	1	1	1
Clear or vacuolated cells <25%	n/a	1	0	1	n/aa	1	n/a
Atypical mitotic figures	n/a	1	0	1	0	1	1
Necrosis	1	1	1	1	1	1	1
Capsular invasion	n/a	1	1	0	0	1	n/a
Nuclear grade III or IV	1	1	1	1	1	1	1
Diffuse architecture	1	1	1	1	1	1	n/a
Venous invasion	n/a	1	0	0	0	1	n/a
Sinusoidal invasion	n/a	1	0	1	0	1	n/a
Specimen type	Biopsy	Resection	Resection	Resection	Resection	Resection	Biopsy
Score: Weiss (modified Weiss)	4/9 (3/7)	9/9 (7/7)	5/9 (4/7)	7/9 (6/7)	*	9/9 (7/7)	4/9 (4/7)

(M) Male, (F) female, (WGTA) whole-genome and transcriptome analysis, (HPF) high-power field, (n/a) not applicable, (*) oncocytic variant.

Figure 1.

Treatment and outcomes. The gray bar represents the time from diagnosis of metastatic disease to the last follow-up or death, with colored sections representing the time on treatment for different systemic agents.

Genomic Analyses

Sequencing coverage, tumor content, and tumor ploidy for all seven cases is shown in Table 2. All patients had advanced disease, with four of the biospecimens for sequencing acquired from the primary site and three from metastatic sites. Somatic alterations including small mutations/indels, copy-number alterations (CNAs), structural variant data, and the gene expression (the number of reads per kilobase of transcript per million mapped reads [RPKM]) data for each case presented in the paper have been deposited at the Michael Smith Genome Sciences Centre at BC Cancer (BCGSC; https://bcgsc.ca/downloads/POG_ACC/) and are highlighted below.

Table 2.

Sequencing coverage, tumor content, and tumor ploidy

Case ID	Tumor genome sequencing coverage	RNA sequencing coverage (no. of paired reads)	Tumor content (%)	Ploidy	Ts/Tv
Case1	74×	N/A	40	Triploid	0.691
Case2	98×	331M	86	Triploid	0.749
Case3	84×	167M	54	Diploid	0.541
Case4	91×	192M	80	Triploid	0.79
Case5	59×	213M	84	Triploid	1.177
Case6	107×	286M	80	Tetraploid	0.964
Case7	96×	208M	80	Diploid	2.005

(Ts/Tv) Transition to transversion ratio, (N/A) not applicable, (M) million.

Somatic Copy-Number Alterations and Loss of Heterozygosity

Large-scale LOH affecting multiple whole chromosomes was observed in six of the seven samples (Fig. 2A). Chromosome 17 was affected in all cases (with lesser extent in Case 2) and Chromosomes 2, 4, and 6 were also largely affected in five cases. Frequent copy-number gains were observed in cases 4, 6, and 7 and these were accompanied by frequent copy-number losses in Cases 4 and 6 (Fig. 2B). Large-scale copy-number losses were also detected in Case 5. Gene amplification events were not seen in Case 7, and Cases 1 and 2 have very few amplification events. In Case 5, arm-level amplification was observed on Chromosome 12 (Ch12p) together with focal amplification events on the other arm (Ch12q).

Figure 2.

Loss of heterozygosity (LOH) events (A) and copy-number alterations (B) in seven ACC tumors. LOH events are indicated in green, copy gains are indicated in red, and copy losses are indicated in blue.

Evaluation of genes affected by copy gains and losses revealed recurrent gains affecting KDM5A. A gain of two or more copies of the KDM5A gene on Chromosome 12 was observed in four cases and a single copy gain was observed in one case (Case 4) (Fig. 3A). KDM5A encodes a demethylase that antagonizes the histone methyltransferase activity promoted by MEN1, a gene that is recurrently altered in ACC. KDM5A gene expression was higher in cases with multiple copies gains, supporting the hypothesis that KDM5A activity may be elevated and contribute to epigenetic dysregulation in cases with KDM5A copy gains (Fig. 3B). TERT promoter hotspot mutations were frequently identified in primary ACC tumors (Liu et al. 2014). We did not observe TERT hotspot mutations; however, all samples had copy gains affecting TERT on Chromosome 5 and one sample was found to have a somatic variant of unknown significance affecting TERT (R447C) (Fig. 3A). Copy gains on Chromosome 16 affecting TERF2, a member of the telomere nucleoprotein complex involved in telomerase regulation and telomere maintenance, was also observed in four cases.

Figure 3.

(A) OncoPrint describing alterations in previously identified driver genes and in genes involved in signaling processes frequently affected in our cohort. (B) KDM5A expression in cases with (≥2) and without (<2) multiple copy gains. (C) homologous recombination deficiency (HRD) score, single base substitution signatures (SBS3 and SBS8), structural variant signatures (SV signature 3 and 5), fraction of deletions with microhomology, and the aggregate HRDetect score in the cohort. (D) Hierarchical clustering of ACC tumors with transcriptome profiles (biopsy site is indicated for each case).

Copy losses affecting the genes involved in the cell cycle were observed including homozygous loss of CDKN2A and CDKN2B (Case 4) and RB1 (Cases 2 and 6). Copy gains, but not amplification, of CCNE1 were noted in Cases 6 and 7. NF1, another recurrently altered gene in ACC, was subjected to a homozygous deletion in one case (Case 3).

Single-Nucleotide Variants and Indels

The variants highlighted in this section are listed in Table 3. Several alterations predicted to affect diverse DNA repair pathways were observed. Somatic alterations in genes involved in the homologous recombination repair pathways include a homozygous in-frame deletion of ATM (Case 6), a heterozygous BRCA2 nonsense mutation (Case 4); other alterations affecting DNA repair pathways include a pathogenic germline mutation affecting a base excision repair gene, MUTYH (Case 1), and a homozygous frameshift mutation affecting MSH6 and a homozygous somatic mutation in MLH1 (Case 7; Fig. 3A).

Table 3.

Variants of interest

Gene	Chromosome	HGVS DNA reference	HGVS protein reference	Variant type	Predicted effect	dbSNP/dbVar ID	Genotype	ClinVar accession number
CTNNB1	3	c.133T > G	p.Ser45Ala	SNV	Missense	rs121913407	Het
CTNNB1	3	c.133T > G	p.Ser45Ala	SNV	Missense	rs121913407	Hom
TERT	5	c.1339C > T	p.Arg447Cys	SNV	Missense	rs1304418053	Het
ATM	11	c.1837_1839delGTG	p.Val613del	DEL	In-frame deletion		Hom	SCV002062047
BRCA2	13	c.2368G > T	p.Glu790*	SNV	Nonsense	rs398122746	Het
MUTYH	1	c.1145G > A	p.G382D	SNV	Missense	rs36053993	Het
MSH6	2	c.3644T > A	p.Leu1215*	SNV	Nonsense		Hom	SCV002062052
MLH1	3	c.884G > T	p.Ser295Ile	SNV	Splice region variant		Hom	SCV002062048
TP53	17	c.499C > T	p.Gln167*	SNV	Nonsense	rs1555526097	Het
TP53	17	c.742C > T	p.Arg248Trp	SNV	Missense	rs121912651	Hom
ARID1B	6	c.4084C > T	p.Gln1362*	SNV	Nonsense		Hom	SCV002062056
SMARCA4	19	c.276G > C	p.Gln92His	SNV	Missense		Het	SCV002062051
SMARCC2	12	c.3221delC	p.Pro1074fs	SNV	Frameshift variant	rs752788954	Het
KMT2C	7	c.8390delA	p.Lys2797fs	SNV	Frameshift variant	rs747256476	Het
KMT2A	11	c.2318dupC	p.Ser774fs	SNV	Frameshift variant		Hom	SCV002062055
ATRX	X	c.6235C > T	p.Arg2079*	SNV	Nonsense		Het	SCV002062053
ATRX	X	c.4377_4379delGGA	p.Glu1460del	DEL	Disruptive in-frame deletion	rs587780286	Het
ATRX	X	c.3904delA	p.Arg1302fs	SNV	Frameshift variant		Het	SCV002062054
PIK3CA	3	c.2086G > A	p.Gly696Arg	SNV	Missense		Het	SCV002062049
JAK2	9	c.2385G > T	p.Arg795Ser	SNV	Missense		Hom	SCV002062050

(HGVS) Human Genome Variation Society, (SNV) single-nucleotide variant, (Het) heterozygous, (Hom) homozygous.

Somatic mutations affecting CTNNB1 and TP53 are the most recurrent driver events in ACC (Zheng et al. 2016) and were observed in Cases 1, 5, 6, and 7, and co-occurred in Case 1 (Fig. 3A). The TP53 variants comprised a nonsense mutation, a missense mutation associated with loss of function and a homozygous deletion. The recurrent gain-of-function S45P mutation affecting CTNNB1 was detected in two samples, with a third sample affected by an in-frame deletion of exons 2–3 of CTNBB1 that deletes the first 80 amino acids of the protein, leading to the loss of inhibitory amino-terminal phosphorylation sites (Maharjan et al. 2018).

Mutations affecting MEN1, a gene involved in epigenetic regulation through its association with the histone methyltransferase MLL1/KMT2A, are recurrent in ACC but were not observed in this cohort. On the other hand, we observed several somatic mutations affecting SWI/SNF chromatin remodeling genes (ARID1B, SMARCA4, SMARCC2) and histone methyltransferase genes (KMT2C, KMT2A, and ATRX) (Fig. 3A).

Alterations affecting genes involved in oncogenic signaling pathways, including a somatic missense mutation of unknown significance in PIK3CA (Case 2), were detected (Fig. 3A). A novel homozygous missense mutation affecting the kinase domain of JAK2 was also observed in Case 5 and may affect JAK2-mediated cytokine signaling (Fig. 3A).

Structural Variant Analysis

Few structural variants affecting cancer genes were observed in this case series. Case 1 was found to harbor a translocation between CDKN2C and FAF1 (1:50978798|1:51436870), which may disrupt the function of CDKN2C involved in cell cycle inhibition. Case 5 harbored a structural variant in RAD52 (12:1029459|12:1095898), which is predicted to cause a duplication of one-half of the gene and the creation of an additional truncated RAD52 protein that lacks the carboxy-terminal 181–412 residues (Fig. 3A). The same case was also found to harbor an interchromosomal translocation affecting STK11 and a noncoding region that was predicted to disrupt the first 300 residues of STK11 and is an inferred loss-of-function mutation. A complete list of high-quality structural variants (see Methods) have been made publicly accessible for download (https://bcgsc.ca/downloads/POG_ACC/).

Transcriptome Analysis

EZH2, histone methyltransferase, and catalytic subunit of the PRC2 complex were highly expressed in two patients compared to the Cancer Genome Atlas (TCGA) ACC data set (96th percentile in Cases 2 and 91st percentile in Case 6). Both cases also had homozygous loss of RB1, in agreement with previous findings that loss of this tumor suppressor promotes EZH2 expression in ACC (median EZH2 expression is 7.4 RPKM in RB1 mutant vs. 2.6 RPKM in RB1 wild-type cases in our cohort) (Drelon et al. 2016).

IGF1R signaling has been shown to be frequently activated in ACCs (Weber et al. 1997), and high expression of IGF1R (Cases 3 and 7; 100th and 95th percentile, respectively), IGF2R (Cases 3, 6, and 7; 100th, 85th, and 89th percentile, respectively), IGF2 (Case 5; 82th percentile) and IGF1 (Cases 2, 3, and 5; 86th, 91st, and 85th percentile, respectively) compared to the TCGA ACC data set was observed, indicating potential increased activation of the IGF1R signaling axis.

Clustering based on genomic, methylation, and transcriptome profiles has revealed molecular subgroups of ACC as well as signatures associated with adrenocortical differentiation (Zheng et al. 2016). Transcriptome clustering using previously defined gene sets identified three different expression subtypes, with the rarest steroid-phenotype-low + proliferation expression subtype lacking in our cohort (Fig. 3A). Three patients (Cases 4, 5, and 7) clustered with the steroid-phenotype-high group, two patients (Cases 2 and 6) with the steroid-phenotype-high + proliferation group, and one patient with the steroid-phenotype-low group. Evaluation of the adrenocortical differentiation gene expression signature score (ADS) revealed four of six patients with positive ADS values supporting the differentiated state of these cases, one of which harbored a CTNNB1 somatic mutation and exhibited a very high ADS value.

Previously described expression signatures were not closely aligned with driver mutations in our cohort; thus, we sought to determine if hierarchical clustering in our cohort may reveal expression profiles associated with driver status. Hierarchical clustering of transcriptome profiles from six patients with available RNA sequencing (RNA-seq) distinguished three pairs (Fig. 3D). The three pairs were characterized by different driver gene mutations, TP53 mutation or loss was detected in pair 1 (Cases 7 and 3), CTNNB1 mutation or deletion was identified in pair 2 (Cases 5 and 6), whereas tumors in pair 3 (Cases 2 and 4) had no mutations in either TP53 or CTNNB1 but had other gene alterations. In addition, pair 1 and 2 were clustered together and were distinguished from tumors with no TP53 and CTNNB1 alterations. These data suggest genomic alterations might lead to gene expression profiles that are characteristic of tumors with specific driver gene mutations.

Integrative Analysis

Our observations, and previous reports in ACC (Zheng et al. 2016; Fojo et al. 2020), have indicated a significant subset of ACC may be associated with homologous recombination (HR) repair deficiency, although the impact of mutations affecting several HR genes is unclear. To investigate the impact of HR mutations in our cohort, we performed mutation signature analysis and investigated homologous recombination repair using HRDetect that aggregates four single-nucleotide variant (SNV) and structural variant (SV) mutation signatures, homologous recombination deficiency (HRD) index, and the proportion of small deletions associated with microhomology (Fig. 3C; Davies et al. 2017; Zhao et al. 2017). Case 5 had a high proportion of mutations attributed to Catalog Of Somatic Mutations in Cancer (COSMIC) single base substitution 3 (SBS3) signature, which previously has been associated with HRD (Figs. 3A and 4A), and HRDetect results were consistent with this tumor being deficient in homologous recombination repair (Fig. 3C), despite the absence of a BRCA1 or BRCA2 pathogenic mutation (Davies et al. 2017). In addition to the structural variant in RAD52, copy losses of several DNA repair genes including BRCA2, CHEK2, and ATM were detected in this case (Fig. 3A). The patient received the PARP inhibitor olaparib and had disease stabilization for 8.5 mo. Case 2 was also predicted to be HRD by HRDetect and similarly lacked a BRCA1/2 alteration. Interestingly, Case 4 with the heterozygous BRCA2 mutation and strong SBS3 signature was predicted to be HR-proficient by HRDetect, in line with the presence of a functional wild-type allele. Apart from HRD-associated signatures, Case 3 and Case 6 had a high proportion of mutations attributed to SBS31, SBS35, and DBS5, all previously associated with prior platinum drug treatment (Fig. 4A,B). Case 3 received platinum treatment prior to metastatic disease (86 d), whereas Case 6, which exhibited very high SBS31 and SBS35 exposures (ranked second and first, respectively, compared against 570 cancer genomes of mixed origins), was treated with platinum for 163 d prior to biopsy. The presence of genome signatures consistent with HRD provides support for involvement of the HR pathway in ACC pathogenesis and may represent a compelling target in ACC as in other tumors of the endocrine system.

Figure 4.

The contribution of previously reported single base substitution (SBS) (A), doublet base substitution (DBS) (B), and indel (ID) (C) mutational signatures (Catalog of Somatic Mutations in Cancer [COSMIC v3.1]) in ACC tumors are compared against 570 cancer genomes of mixed origins (Pleasance et al. 2020).

Case 1 with the pathogenic germline MUTYH mutation was characterized by a high proportion of mutations attributed to SBS18 (ranked 20th compared against 570 cancer genomes), previously associated with MUTYH deficiency (Fig. 4A; Thibodeau et al. 2019). Case 7 with deleterious mutations in MSH6 and MLH1 demonstrated a high mutation burden (whole-genome tumor mutation burden of 28.6 mutations/Mb), and microsatellite instability (39% of microsatellites were unstable; Niu et al. 2014). Mutation signature analysis also revealed a high proportion of mutations attributed to SBS26, SBS44, DBS7, ID2, and ID7, all previously associated with mismatch repair deficiency (MMRd) (Fig. 4; Drost et al. 2017). SBS44 (ranked eighth), DBS7 (ranked sixth), ID2 (ranked 11th), and ID7 (ranked fourth) were particularly high compared to a large, previously sequenced cancer cohort. Interestingly, Case 7 also harbored mutations consistent with signatures SBS7c and SBS38 (Fig. 4A), both of which are primarily observed in melanomas (Saini et al. 2016). We and others have reported SBS38 in non-UV-associated melanoma subtypes, raising the possibility of a common underlying mutational process in ACC and a subset of non-UV-exposed melanomas (Newell et al. 2020; Pleasance et al. 2020). After progressing on EDP chemotherapy and mitotane, the patient was initiated on avelumab, but unfortunately died before completing one cycle because of rapidly progressing disease.

DISCUSSION

There are relatively few studies of comprehensive genomic characterizations of ACC (Zheng et al. 2016; Fojo et al. 2020). Current exome sequencing studies have revealed a few recurrently altered genes in ACC, with examples of variants of uncertain significance impacting multiple genes aligned to several cancer associated pathways. Here we describe a cohort of patients with metastatic ACC profiled using whole-genome and transcriptome analysis (WGTA) and provide detailed clinical annotation.

Our study is limited by the small sample size and heterogeneity of the cohort. Four out of seven samples were obtained from primary adrenal tumors, whereas the remaining three were obtained through biopsies of metastatic lesions. This introduces significant heterogeneity. The fact that four patients received systemic therapy prior to biopsy, including two patients who received cytotoxic chemotherapy, also creates significant bias. Our analysis highlights the genomic diversity of ACC, which makes inferring specific conclusions about the management of this rare disease challenging. The clinical–genomic correlates are presented as hypothesis-generating findings.

Whole-genome sequencing allows for detection of potentially disruptive structural variants, including the STK11 and RAD52 variants detected here, and evaluation of genome-wide mutation signature phenotypes including HRDetect and microsatellite instability. Transcriptome sequencing informs on the impact of observed genomic alterations, differentiation, and revealed distinct expression profiles that were correlated with driver gene mutation status.

Consistent with previous studies, somatic mutations in CTNNB1 and/or TP53 were observed in our cohort, including an exon level deletion affecting CTNNB1 similar to that described in 7% of ACC cases (Maharjan et al. 2018). Although no somatic mutations in other previously reported driver genes were detected in our cohort, we observed the activation of pathways implicated in ACCs. We observed copy-number gains of CCNE1 in six cases and detected homozygous loss of RB1 in two cases and homozygous loss of CDKN2A and CDKN2B in one case, suggesting the deregulation of the cell cycle and suggesting that targeting the cell cycle pathway may be beneficial in a subset of patients. We did not detect MEN1 mutations in our cohort; however, we did observe multiple copy gains and concurrent high expression of KDM5A, a demethylase that antagonizes MEN1-mediated methylation. A single TCGA ACC sample harbored a KDM5A amplification, and this case was also negative for MEN1 mutation (Zheng et al. 2016). Increased KDM5A expression may represent an alternate mechanism of epigenetic dysregulation compared to MEN1 loss of function and is the subject of active investigation as an oncology target (Yang et al. 2021).

In addition to alterations in ACC driver genes, alterations consistent with deregulation of epigenetic signaling and DNA repair processes were observed. We identified somatic mutations in SMARCA4, SMARCC2, and ARID1B, suggesting deregulation of the chromatin remodeling complex, SWI/SNF. Inactivating mutations in genes encoding SWI/SNF subunits are frequently observed in various cancer types and have been associated with deficiency in DNA repair as well as response to multiple targeted therapies (Mittal and Roberts 2020; Ribeiro-Silva et al. 2019), but the role of the SWI/SNF complex in ACC tumorigenesis and therapeutic response has not yet been investigated. In addition, somatic variants were detected in KMT2A and KMT2C histone methyltransferases and in ATRX and DAXX as previously reported in primary ACCs (Zheng et al. 2016). We also observed increased EZH2 expression, particularly in cases with RB1 loss. EZH2 inhibitors are approved for use in EZH2 mutant follicular lymphoma and in SMARCB1/INI1 deficient epitheliod sarcoma (Groisberg 2020) and may be therapeutically beneficial in a subset of SWI/SNF deficient, or EZH2 overexpressing ACCs.

Four out of the seven cases were associated with mutations in DNA repair genes, including somatic alterations in BRCA2, MLH1, MSH6, ATM, and a germline pathogenic mutation affecting MUTYH. Mutations in the DNA-mismatch repair pathway have been observed in primary ACCs (Zheng et al. 2016) and ACC has been shown to be a Lynch syndrome–associated cancer (Raymond et al. 2013). Germline MUTYH mutations are also known to be associated with ACC predisposition; however, germline alterations affecting HR-related genes have only been rarely attributed to ACC (El Ghorayeb et al. 2016; Xie et al. 2018), suggesting a potential subset of ACCs is characterized by HRD derived from somatic alterations as opposed to germline mutation. Two cases in our cohort demonstrating a strong HRD signature lacked mutations in BRCA1/2, suggesting either epigenomic dysregulation or alterations affecting other HR genes may contribute to the HRD in ACC. Profiling mutational signatures associated with HRD may reveal ACC patients that can benefit from platinum-based therapy or PARP inhibitors, as was observed in other tumor types (Zhao et al. 2017; Chopra et al. 2020; Golan et al. 2021).

Despite numerous shared alterations, this cohort also illustrates the high degree of heterogeneity that exists within ACC and represents a challenge to successfully propose targeted therapy for ACCs. One patient (Case 6) had a partial response to combination chemotherapy (etoposide, doxorubicin, and cisplatin), and another patient (Case 5) has had sustained benefit from mitotane and subsequently from olaparib, whereas other treatments only led to minimal clinical benefit. Therapy informed by WGTA analysis was given for Case 5, Case 6, and Case 7. Case 5 had evidence of HRD despite the absence of a BRCA1 or BRCA2 mutation and received the PARP inhibitor olaparib resulting in disease stabilization for 8.5 mo. This patient had a potential disruptive structural variant in RAD52, which participates in HR. In mammals RAD52 deletion alone does not compromise HR repair; however, depletion of human RAD52 is synthetically lethal with mutations in other members of the HR pathways including BRCA1, BRCA2, PALB2, and RAD51 paralogs (Feng et al. 2011; Chun et al. 2013; Lok et al. 2013). In this case the structural variant in RAD52 together with the copy losses of other HR genes might have resulted in the HRD phenotype and partial response to PARPi. The patient showed mixed response to the subsequent carboplatin and paclitaxel combination with disease progression 5 mo after initiation of treatment. Case 6 was enrolled in a phase 1 clinical trial of a TTK inhibitor based on the presence of a CTNNB1 activating mutation and preclinical evidence (Zaman et al. 2017), but did not demonstrate clinical benefit. This patient subsequently received the PARP inhibitor olaparib based on the presence of an in-frame deletion affecting ATM with loss of heterozygosity and low ATM gene expression based on transcriptome and immunohistochemistry, although this sample was not classified as HRD by HRDetect. A phase 2 study in gastric cancer (Bang et al. 2015) suggested improved overall survival for patients treated with a PARP inhibitor and a taxane, when compared to a taxane alone. This benefit was preferentially seen in a cohort of patients with low ATM expression by immunohistochemistry, although a confirmatory phase 3 study was negative (Bang et al. 2017). In this case, the PARP inhibitor olaparib was combined with temozolomide based on the recommended dose and schedule derived from a phase 1 clinical trial in small-cell lung cancer (Farago et al. 2018), but unfortunately the patient did not benefit from this combination. Case 7 had evidence of microsatellite instability and mismatch repair deficiency; and based on these findings the patient was initiated on an immune checkpoint inhibitor (Le et al. 2017), but the patient's clinical status had rapidly deteriorated and she passed away within a few days of initiating treatment without a response assessment.

Overall, the findings of this cohort suggest that molecular profiling of ACC could be utilized to target specific treatments to certain alterations affecting several potentially actionable pathways including DNA repair, cell cycle, oncogenic signaling pathways, and genes involved in epigenetic regulation found in specific subpopulations. Comprehensive genomic and transcriptome analysis of ACC has the potential to reveal therapeutic strategies tailored to patients’ individual tumors. Although the lack of response seen in previous studies of targeted agents may be related to the poor activity of these drugs, it is also possible that the heterogeneity of the underlying disease contributed to the disappointing results observed. To better guide future studies, more information is required on the molecular makeup of metastatic ACC.

METHODS

Sample Collection and Preparation

Whole-genome sequencing of tumors and matched normal tissue was performed in seven cases; four of the biopsies were obtained from primary adrenal tumors, whereas the remaining three were biopsies of metastatic lesions, as outlined in Table 1. Whole-transcriptome sequencing was performed on six of the seven tumors. Patients with locally advanced or metastatic ACC with a life expectancy of greater than 6 mo were enrolled on the Personalized OncoGenomics (POG) project (ClinicalTrials.gov identified: NCT02155621) at BC Cancer, in Vancouver, Canada (Laskin et al. 2015; Jones et al. 2017).

Tumor sections were reviewed by a pathologist to determine tumor content and cellularity, and remaining sections were used for DNA and RNA preparation. Peripheral blood was collected for use as a normal tissue reference and pathogenic and likely pathogenic germline variants in 98 known cancer predisposition genes were analyzed as outlined in Dixon et al. (2020). DNA and RNA were coextracted for the construction of polymerase chain reaction (PCR)-free DNA libraries and poly(A) selected RNA libraries as described previously (Sheffield et al. 2015).

Sequencing and Bioinformatics

Tumor genomes sequencing was performed to 80× coverage and normal genomes to 40× coverage on Illumina HiSeq platform with 125- or 150-bp paired-end reads. Illumina HiSeq2500 or NextSeq500 were used for RNA sequencing with 150–200 million 75-bp paired-end reads. Sequence reads were aligned to the human reference genome (hg19) by the BWA tool (Li and Durbin 2010) and somatic SNVs and small indels were identified using SAMtools (v0.1.17) (Li et al. 2009) and Strelka (v1.0.6) (Saunders et al. 2012). Regions of copy-number alteration (CNA) and LOHs were determined using CNAseq (v0.0.6) (Jones et al. 2010) and APOLLOH (v0.1.1) (Ha et al. 2012), respectively. Amplification was defined as the ploidy-corrected copy-number gain being greater than the estimated tumor ploidy. Copy gains are defined as an increase in gene copies equivalent to or less than the estimated average ploidy. SVs in RNA-seq data were detected by ABySS v1.3.4 and TransABySS (v1.4.10) (Birol et al. 2009; Simpson et al. 2009) and Chimerascan (v0.4.5) (Iyer et al. 2011) and DeFUSE (v0.6.2) (McPherson et al. 2011); SVs in the DNA were identified using ABySS and TransABySS followed by Manta v1.0.0 (Chen et al. 2016) and Delly v0.7.3 (Rausch et al. 2012). SV calls from multiple algorithms were merged into a consensus caller MAVIS (v2.1.1) that performs a subsequent validation by local assembly (Reisle et al. 2019). Structural variants were filtered for events that were called by multiple tools and were supported by an assembled contig and were deposited for download (https://bcgsc.ca/downloads/POG_ACC/). The contribution of previously reported mutational signatures in COSMIC v3.1 (https://cancer.sanger.ac.uk/signatures; Alexandrov et al. 2020) was calculated using nonnegative least squares optimization. HRDetect scores were computed using a logistic regression model described previously in Davies et al. (2017). Six HRD-associated mutation signatures were used in the algorithm: (i) Signature 3, (ii) Signature 8, (iii) SV signature 3, (iv) SV signature 5, (v) the HRD index, and (vi) the fraction of deletions with microhomology. Somatic SNVs called by Strelka (v1.0.6) were used for single base substitution signature calculation. The contribution of previously reported mutational signatures in COSMIC v3.1 (https://cancer.sanger.ac.uk/cosmic/signatures) was calculated using Monte Carlo Markov chain (MCMC) sampling (https://github.com/eyzhao/SignIT). Somatic, DNA-derived SVs identified by the consensus caller MAVIS were used for SV signature calculation (Reisle et al. 2019). MAVIS calls that were detected by more than one tool and for which the contig could be assembled were included in the analysis, and the contribution of the previously reported SV mutational signatures was calculated using MCMC sampling (https://github.com/eyzhao/SignIT) (Nik-Zainal et al. 2016). HRD index was computed as the arithmetic sum of LOH, telomeric-allelic imbalance (TAI), and large scale transitions (LST) scores (Zhao et al. 2017). All signatures were normalized and log transformed, and a threshold of 0.7 was used as previously reported (Davies et al. 2017; Zhao et al. 2017).

Microsatellite instability (MSI) was calculated by MSIsensor (v0.2) (Niu et al. 2014).

Gene Expression Analysis

A transcriptome of an average depth of 200 million reads was generated from the tumor sample. The RNA sequencing reads were analyzed with JAGuaR (Butterfield et al. 2014) to include alignments to a database of exon junction sequences and subsequent repositioning onto the genomic reference. The RNA sequencing data were processed using the Genome Sciences Centre's wtss (whole-transcriptome shotgun sequencing) pipeline coverage analysis (version 1.1) with the “stranded” option to determine gene and exon read counts and normalized expression level. The level of expression of each gene was determined as the number of RPKM.

Gene expression level was evaluated by converting RPKM of select genes into percentile ranks against transcriptomes of ACC tumors from The Cancer Genome Atlas project (https://portal.gdc.cancer.gov/). Expression subtypes classification shown in Figure 3A (Exp subtype track) was inferred by calculating spearman correlation factor between individual TCGA samples and POG samples. Hierarchical clustering shown in Figure 3D was performed after selecting the top 2500 genes according to the magnitude (from largest to smallest in absolute values) of their coefficients contributing to the first principal component analysis (PCA) component and after conversion of the RPKM values to z-scores using the expression data of six patients with RNA-seq results. Agglomerative clustering method was used to obtain clusters, computing the distance matrix with the Euclidean distance. The adrenocortical differentiation score (ADS) shown in Figure 3A (ADS score track) was calculated by using 25 gene expression as described by Zheng et al. (2016). To obtain the expression subtypes described in Zheng et al. (2016), kmeans clustering was performed on the publicly available mRNA data of 76 TCGA ACC samples with expression subtypes classification (downloaded from https://pancanatlas.xenahubs.net) and the POG mRNA data that was calculated using the TOIL RNA-seq pipeline (Vivian et al. 2017) to minimize computational batch effects between POG samples and recomputed TCGA samples. The TOIL pipeline uses CutAdapt for adapter trimming and STAR (Dobin et al. 2013) to generate alignments and performs quantification using RSEM (Li and Dewey 2011) and Kallisto (Bray et al. 2016).

ADDITIONAL INFORMATION

Data Deposition and Access

The data sets generated and analyzed during the current study have been deposited at the Michael Smith Genome Sciences Centre at BC Cancer (BCGSC; https://bcgsc.ca/downloads/POG_ACC/). The genomic and transcriptome sequencing data sets have been deposited at the European Genome-phenome Archive (EGA; https://ega-archive.org/) as part of the study EGAS00001001159 (Patient 1 EGAD00001002037, Patient 2 EGAD00001002049, Patient 3 EGAD00001003053, Patient 4 EGAD00001004905, Patient 5 EGAD00001003676, Patient 6 EGAD00001005844, Patient 7 EGAD00001005897).

Ethics Statement

The study was conducted in accordance with the principles of the Declaration of Helsinki and with approval of the Research Ethics Board at BC Cancer, and all participants provided written consent.

Acknowledgments

This work would not be possible without the participation of our patients and families, the POG team, Canada's Michael Smith Genome Sciences Centre technical platforms, the generous support of the BC Cancer Foundation and their donors and Genome British Columbia. The results published here are in part based upon analyses of data generated by the following projects and obtained from dbGaP (http://www.ncbi.nlm.nih.gov/gap): TCGA managed by the National Cancer Institute and National Human Genome Research Institute (http:// cancergenome.nih.gov). Data from the Pan-Cancer Analysis of Whole Genomes (PCAWG) managed by the International Cancer Genome Consortium was retrieved from https://dcc.icgc.org/pcawg.

Author Contributions

J.-M.L. and V.C. conceptualized the study. Analyses were performed by J.-M.L., V.C., L.M.W., L.C., and G.W. Provision of patient samples and curation of data was conducted by J.-M.L., C.K.K., D.J.R., and J.L. The original draft was written by J.-M.L., V.C., L.M.W., and M.A.M., J.L., C.K.K., D.J.R., and S.J.M.J. reviewed and edited the manuscript. Funding was acquired by M.M., J.L., and S.J.M.J.

Funding

This work was supported by the BC Cancer Foundation and Genome British Columbia (project B20POG).

Competing Interest Statement

The authors have declared no competing interest.