Multi-omic analysis of SDHB-deficient pheochromocytomas and paragangliomas identifies metastasis and treatment-related molecular profiles

Flynn, Aidan; Pattison, Andrew D.; Balachander, Shiva; Boehm, Emma; Bowen, Blake; Dwight, Trisha; Rossello, Fernando J.; Hofmann, Oliver; Martelotto, Luciano; Zethoven, Maia; Kirschner, Lawrence S.; Else, Tobias; Fishbein, Lauren; Gill, Anthony J.; Tischler, Arthur S.; Giordano, Thomas; Prodanov, Tamara; Noble, Jane R.; Reddel, Roger R.; Trainer, Alison H.; Ghayee, Hans Kumar; Bourdeau, Isabelle; Elston, Marianne; Ishak, Diana; Ngeow Yuen Yie, Joanne; Hicks, Rodney J.; Crona, Joakim; Åkerström, Tobias; Stålberg, Peter; Dahia, Patricia; Grimmond, Sean; Clifton-Bligh, Roderick; Pacak, Karel; Tothill, Richard W.

doi:10.1038/s41467-025-57595-y

Download PDF

Article
Open access
Published: 17 March 2025

Multi-omic analysis of SDHB-deficient pheochromocytomas and paragangliomas identifies metastasis and treatment-related molecular profiles

Nature Communications volume 16, Article number: 2632 (2025) Cite this article

6061 Accesses
8 Citations
4 Altmetric
Metrics details

Subjects

Abstract

Hereditary SDHB-mutant pheochromocytomas (PC) and paragangliomas (PG) are rare tumours with a high propensity to metastasize although their clinical behaviour is unpredictable. To characterize the genomic landscape of these tumours and identify metastasis biomarkers, we perform multi-omic analysis on 94 tumours from 79 patients using seven molecular methods. Sympathetic (chromaffin cell) and parasympathetic (non-chromaffin cell) PCPG have distinct molecular profiles reflecting their cell-of-origin and biochemical profile. TERT and ATRX-alterations are associated with metastatic PCPG and these tumours have an increased mutation load, and distinct transcriptional and telomeric features. Most PCPG have quiet genomes with some rare co-operative driver events, including EPAS1/HIF-2α mutations. Two mechanisms of acquired resistance to DNA alkylating chemotherapies are identifiable; MGMT overexpression and mismatch repair-deficiency causing hypermutation. Our comprehensive multi-omic analysis of SDHB-mutant PCPG therefore identifies features of metastatic disease and treatment response, expanding our understanding of these rare neuroendocrine tumours.

Single-nuclei and bulk-tissue gene-expression analysis of pheochromocytoma and paraganglioma links disease subtypes with tumor microenvironment

Article Open access 21 October 2022

Cytoplasmic HIF-2α as tissue biomarker to identify metastatic sympathetic paraganglioma

Article Open access 18 July 2023

Single-cell chromosome and bulk transcriptome analysis as a diagnostic tool to differentiate between localized and metastatic pheochromocytoma and sympathetic paraganglioma

Article 03 May 2025

Introduction

Pheochromocytomas (PC) and paragangliomas (PG) are heritable neuroendocrine neoplasms arising from the sympathetic and parasympathetic nervous system. They are clinically remarkable for catecholamine hypersecretion (e.g., dopamine, noradrenaline, adrenaline) causing morbidity and occasionally death from cardiovascular sequelae^1,2. PCPG are typically slow growing, but metastases develop in 10–20% of patients³. Current biomarkers are inadequate for predicting progression to metastatic disease, which is currently incurable.

Metastatic risk, which guides surveillance recommendations^4,5,6, is highest in patients carrying constitutional pathogenic variants in SDHB. SDHB-deficiency and cellular succinate accumulation inhibits prolyl hydroxylase domain proteins (PHD1/2) leading to oxygen-independent stabilization of hypoxia inducible factor alpha (HIF-α) and HIF target gene activation - a phenomenon known as pseudohypoxia^7,8. Succinate excess also inhibits other 2-oxoglutarate-dependent dioxygenases, including TET2 and jumonji histone lysine demethylases, causing genome-wide hypermethylation and epigenetic reprogramming^8,9. While the causative role of SDHB-deficiency in metastatic PCPG is well founded, other genomic changes also likely contribute, such as ATRX and TERT mutations^{10,11,12,13,14,15,16}. However, most PCPG genomic studies to date have used whole-exome sequencing, limiting detection of non-coding variants, structural alterations, telomeric features and DNA mutation patterns.

Herein, we conduct an integrated multi-omic analysis of an internationally sourced and clinically well-annotated cohort of SDHB-mutant PCPG including parasympathetic head and neck PG (HN-PG) - a relatively understudied group to date. Importantly, paired primary and metastatic tumours are analysed in a subset of cases. We make important observations with respect to cell-of-origin, treatment response and clinical outcome and confirmed the presence of TERT/ATRX mutations and other features in metastatic cases.

Results

Whole genome and multi-omic analysis of SDHB-associated PCPG

To characterize the genomic landscape of SDHB-mutant PCPG, we performed multi-omic analysis of an international (A5 consortium) cohort representing 94 primary and/or metastatic PCPG tumours from 79 patients (Fig. 1A) (see Supplementary Data 1 and 2 for clinical details). Analysis included whole-genome sequencing (WGS) of tumour and matched blood (n = 94), whole transcriptome sequencing (WTS) (n = 91), small-RNA seq (n = 90), DNA methylation profiling (n = 93) and C-circle analysis (n = 89). Droplet-based (10x) single nuclei (sn)RNA-seq and snATAC-seq was applied in a subset of cases (Fig. 1B).

**Fig. 1: Multi-omic analysis of SDHB-mutated PCPG.**

Review of histopathology and clinical data confirmed disease diagnosis in all cases. Paired synchronous/metachronous primaries were analysed in four cases (E128, E136, E159, E229) and clonal independence confirmed by discordant somatic profiles (Supplementary Fig. 1 and 2). With respect to clinical behaviour, 40 tumours from 37 patients were ostensibly non-metastatic with at least twelve-months of clinical follow-up (median: 60 months, range 12–456). An additional eight primary PCPG from eight patients had less than 12 months of clinical follow-up without metastases reported (termed: short clinical follow-up). Thirty-four of 79 patients had confirmed metastatic disease. Twenty-eight primary PCPG were analysed from 26 patients who developed metastases and in six patients at least one paired metastasis was available, confirming the clonal link between the primary and metastatic tumours (termed: metastatic primary) (Fig. 1C). For 21 primary PCPG from 20 patients, metastasis was reported but metastatic tissue was not available (termed: primary - metastasis reported) (Fig. 1C). Finally, in eight patients, one or more metastases were analysed but a primary tumour was not analysed.

WGS confirmed pathogenic germline SDHB mutations consisting of missense (n = 41), nonsense (n = 21), frameshift insertions/deletions (n = 6), large deletion events (n = 7), and donor splice-site mutations (n = 4) (Fig. 1D). Two recurrent large deletion events were identified in two or more unrelated patients (chr1:g.17052136-17054668_del, chr1:g.17048756-17064432_del). WGS showed variable length homologous DNA sequence at the breakpoint ends of a subset of the large deletions, implicating homologous recombination as a potential cause of the germline deletions (Supplementary Fig. 3). Somatic loss-of-heterozygosity (LOH) at 1p36.13 with or without somatic copy-loss was detected in all tumours and SDHB immunohistochemistry confirmed tumoural SDH-deficiency.

Molecular profiling confirms genotype-subtype and cell-of-origin transcriptional profile

Molecular subtyping is known to reflect the genotypic features of PCPG^{13,17,18,19,20,21}. To confirm genotype-phenotype relationships among A5 tumours we applied UMAP clustering to WTS (n = 91), small-RNA-seq (n = 90) and DNA methylation (n = 93) datasets. Published RNA-seq data across a spectrum of PCPG genotypes was included for comparison^13,14,20. Here, we used the C1/C2 annotation based on seven PCPG gene-expression subtypes previously defined by single-nuclei and bulk gene-expression analysis²². As expected, all A5 SDHB-mutant tumours clustered among C1A (SDHx) tumours by WTS (Fig. 2A). UMAP of the small-RNA-seq data showed a similar clustering pattern; however, four A5 tumours from two patients clustered with an outlier group consisting of unrelated genotypes (Fig. 2B). This small-RNA-seq outlier group included C2Bi (MAX) PCPG that frequently exhibit loss-of-heterozygosity of chr14q and silencing of the imprinted DLK1-MEG3 miRNA cluster located on chr14q31-32²⁰ (Supplementary Fig. 4). By DNA methylation profiling, SDHB-mutant PCPG were distinguished from other PCPG genotypes, with the exception of one tumour (E229-P1) (Fig. 2C), and SDHB-mutant tumours showed genome-wide hypermethylation (Supplementary Fig. 5), consistent with a prior report⁹. Unsupervised clustering therefore confirmed the expected molecular profile of the A5 SDHB-mutant tumours as well as the pathogenicity of constitutional SDHB variants.

**Fig. 2: Genomic profiling separates sympathetic and parasympathetic PGL.**

UMAP clustering of only the A5 tumours showed clear separation of parasympathetic (non-chromaffin cell) HN-PG from sympathetic (chromaffin cell) PCPG (Fig. 2D–F). Two mediastinal dopamine-secreting PG (E148-P1, E155-P1) clustered with the HN-PG group, suggesting these tumours may be non-chromaffin cell paragangliomas²³ while one HN-PG (E185-P1) and one abdominal-thoracic PG (E128-P2) showed variable clustering. Supervised differential expression analysis between HN-PG and sympathetic PCPG identified 4615 genes and 36 microRNAs (adjusted p-value < 0.05, log fold change > 1) (Supplementary Data 3 and 4). HN-PG and the suspected non-chromaffin mediastinal tumours had low expression of the chromaffin cell marker CARTPT²⁴, neural transcriptional regulators (TFAP2B, TOX3, GATA3, POU4F, NEUROG2, PAX2), HOX genes (HOXA1-10, HOXB4, HOXB6-9, HOXC4-HOXC13), and the long non-coding RNA HOTAIR (Fig. 2H, Supplementary Fig. 6). Similar to the HOX gene clusters, clustered genes on chr12q24.32 including TMEM132C and lncRNAs were also lowly expressed. Concordant DNA methylation patterns were also observed (Fig. 2G, I, Supplementary Fig. 7).

HN-PG are predominantly biochemically silent and accordingly catecholamine biosynthesis genes including TH and DBH were lowly expressed (Fig. 2J). Although most sympathetic PCPG are biochemically active, a small number of A5 cases had normal catecholamine levels and such tumours also had low TH expression. Some patients with HN-PG/non-chromaffin cell tumours and sympathetic PCPG had elevated plasma 3-methoxytyramine (3-MT) (a dopamine metabolite^25,26) and these tumours had low DBH expression - DBH being necessary to convert dopamine to norepinephrine. Low TH and DBH expression corresponded with gene-promoter methylation (Supplementary Fig. 8). Meanwhile, variable expression of the norepinephrine transporter (SLC6A2) was observed with gene-promoter methylation corresponding with low expression in sympathetic PCPG (Supplementary Fig. 9). Low SLC6A2 (NET) expression is clinically relevant since NET is required for uptake of theranostic agents e.g.,^123/124/131I-meta-iodobenzylguanidine (MIBG)^27,28.

With respect to neurotransmitters, TPH1 (tryptophan hydroxylase), essential for serotonin (5-HT) synthesis, was overexpressed in the HN-PG/non-chromaffin cell tumours (Fig. 2G, Supplementary Fig. 10). TPH is known to be expressed in glomus cells of the carotid body (presumed to be the HN-PG cell-of-origin)²⁹, and 5-HT storage vesicles have been observed in some HN-PG, although urinary 5-HIAA is not elevated in most HN-PG patients³⁰. Meanwhile, while immunohistochemical detection of choline acetyltransferase (CHAT) has been reported in HN-PG³¹, CHAT mRNA was not overexpressed in these tumours (Supplementary Fig. 11).

Genome-wide mutational features of SDHB-deficient PCPG

Consistent with prior studies, therapy naive PCPG (n = 84) had a low mutation burden (median =0.32, range 0.03–4.24 mut/Mb)^11,13,14,20. Single-base substitution (SBS) signature analysis (COSMICv3) showed predominantly age-related or clock-like mutational signatures (SBS1, SBS5) (Fig. 3A). With respect to doublet base substitutions (DBS) the predominant patterns were DBS 2, 4, 9 and insertion/deletion (ID) signatures ID 5, 8 and 9. With respect to systemic cytotoxic treatments, 10/94 (12%) of tumours were resected post-treatment (Supplementary Fig. 12) and chemotherapy-related mutational signatures were observed in some cases including SBS11 (DNA alkylation) (discussed below).

**Fig. 3: Genomic and clinical features of PCPG.**

PCPG tumours showed recurrent patterns of chromosomal changes similar to observations in previous studies^13,14,18 including loss of chr1p, chr2q, chr3, chr8p, chr11, chr17p and chr22 (GISTIC q < 0.05) (Fig. 3B, Supplementary Fig. 13). Significant copy-number gains were observed on chr1q in sympathetic (chromaffin cell) PCPG, while parasympathetic (HN-PG/non-chromaffin cell) tumours had frequent loss of chr1q. A minority of tumours had chromothripsis (n = 6) or genome-doubling (n = 9), as previously reported in PCPG^13,14,20. Notably, we identified significant chr7 gains in the metastatic group (GISTIC q < 0.05). Sub-clonal copy number events were observed in several samples by WGS and copy number inference from sn-RNA-seq (Supplementary Fig. 14A–D), with primary tumours exhibiting greater clonal diversity than metastases (One-sided Student’s t-test, p < 0.01) (Supplementary Fig. 14E).

Few genes were recurrently altered and only 15 cancer-related genes had > 1 predicted pathogenic mutation using Cancer Genome Interpreter³² (Fig. 3C, Supplementary Data 5). Consistent with prior studies, TERT and ATRX were the most frequently mutated genes^11,13,15 (discussed in detail below). ATRX mutations were significant drivers using the dN/dS³³ method while TERT promoter (pTERT) mutations were significant using OncoDriveFML³⁴. Two tumours had somatic EPAS1 (HIF-2α) pathogenic variants near codon 531 (prolyl) targeted by PHD2 (Fig. 3D). One EPAS1 mutation (E230-P1, NP_001421.2:p.Tyr532Cys) was previously described in a PCPG patient with polycythemia^35,36, while a second EPAS1 mutation (E183-P1, NP_001421.2:p.Asp539Gly) involved the same codon as a previously reported EPAS1 pathogenic variant³⁷. Other recurrently mutated cancer genes included chromatin modifiers KMT2C (n = 2) and BCOR (n = 2). TP53 somatic mutations were found in three metastatic tumours consistent with a known low frequency in PCPG^13,14,20. Finally, a MYCN hotspot mutation was found in one tumour (E129-P1, NP_005369.2:p.Pro44Leu), as previously reported in another study³⁸.

Somatic TERT mutations, structural rearrangement and promoter methylation

Hotspot pTERT mutations (chr5:1295113 G > A, chr5:1295135 G > A, GRCh38) were identified in 14 metastatic cases (17.7%). Targeted deep sequencing of pTERT excluded potential sub-clonal mutations in the other cases. snATAC-seq data showed a promoter peak central to the de novo ETS binding site³⁹ in three pTERT mutants tested and all snATAC-seq reads had the somatic variant, confirming open chromatin (Fig. 4A). Chromothripsis of chromosome 5 created an AFF4-TERT fusion in one case (E123) (Fig. 4B). WGS supported a DNA breakpoint in AFF4 3-prime of untranslated exon 1 and DNA sequence 390 bp upstream of the TERT transcription start site (Fig. 4Bii). RNA-seq split-reads supported a fusion transcript involving AFF4 exon 1 spliced to TERT exon 3 (Fig. 4Biii). This validates our prior report of TERT structural rearrangements in PCPG¹⁶.

**Fig. 4: TERT promoter mutations and structural alterations.**

All TERT altered cases had elevated TERT expression with the highest expression observed in the AFF4:TERT fusion tumour. TERT copy-number gain did not correlate with TERT overexpression in the absence of TERT mutation or fusion events (Fig. 4C). We found TERT promoter hypermethylation including the THOR domain (probes cg11625005, cg10767223) in TERT altered PCPG, as previously described¹⁵ (Fig. 4D). TERT promoter hypermethylation and concordant overexpression was also observed in a single metastatic TERT wild-type tumour (E171-M1) (Fig. 4D).

ATRX alterations are associated with the ALT phenotype

Pathogenic ATRX mutations included nonsense and splice-site mutations, small frameshift insertions or deletions, and a large deletion in one case (Fig. 5A). One missense ATRX mutation (unknown significance) was detected (NM_000489.6:c.6778 C > T, NP_000480.3:p.H2260Y). ATRX expression was low in most ATRX altered cases, presumably due to nonsense-mediated decay (Supplementary Fig. 15).

**Fig. 5: Loss of *ATRX* leads to telomere dysregulation.**

Alternative lengthening of telomeres (ALT) is a known feature of ATRX altered PCPG and other cancer types¹². C-circle analysis was done to predict ALT activity⁴⁰ and C-circles were detected in all ATRX altered and one TERT altered primary tumour (E143-P1) (discussed below)(Fig. 5B). As expected, the telomeric tumour:blood ratio, a proxy for tumour telomere length, was higher in ATRX altered compared to TERT altered or TERT/ATRX wild-type tumours (Two-sided Student’s t-test p-value < 0.00025) with no significant difference in telomeric DNA content observed between the latter two groups. Long telomeres in the absence of C-circles were detected in four ATRX/TERT wild-type PCPG, including two HN-PG.

ALT features can include aberrant telomeric variant repeat (TVR) usage, intrachromosomal insertion of telomeric DNA and expression of telomeric repeat-containing RNA (TERRA)⁴¹. Relative TVR usage was altered in C-circle positive versus negative tumours (Fig. 5C). When normalized to total read count (rather than telomeric read count) only three TVRs (TTAGGG, GTTGGG, GTAGGG) remained significant after false-discovery rate correction (Benjamini & Hochberg, p < 0.1). When normalized to the telomeric reads the canonical TVR (TTAGGG) was proportionally more abundant in C-circle positive compared to negative tumours. Only one intrachromosomal telomeric insertion was detected across the cohort, occurring in ATRX altered case E198-M1. TERRA expression was elevated in ATRX altered tumours; however, the difference did not reach statistical significance (Two-sided Student’s t-test p > 0.05)(Fig. 5D). It should be noted that the tag counting method of telomere quantitation cannot differentiate reads originating from centromeric and telomeric regions.

TERT/ATRX alterations are associated with metastatic progression and late somatic events

TERT/ATRX alterations were found in 21/34 (61%) metastatic cases and significantly associated with metastasis (Pearson’s Chi-squared p-value 9.349e-08). Only one non-metastatic tumour (E156-P1) had a pTERT mutation; however, this case had a short clinical follow-up (Fig. 6A). Evidence for late acquisition of pTERT mutations was found in two cases: 1) E146 where the primary was TERT wild-type but the paired metastasis was pTERT mutant and 2) case E143 where a subclonal pTERT mutation was detected in the primary (VAF = 10.3%) but higher VAF was detected in two paired metastases (E143-M1 VAF = 33%, E143-M2 VAF = 28%), indicating all metastatic cells were pTERT mutant (Fig. 6B). Interestingly, the subclonal pTERT mutant primary (E143) also had detectable C-circles and longer telomeres, but these ALT features were absent in both metastases. The data therefore suggests the pTERT mutation and ALT features had co-existed but were likely within different cell populations, while only the pTERT mutant cells metastasized.

**Fig. 6: *TERT*/*ATRX*-alterations and their association with metastatic progression.**

Ten primary tumours from patients with metastatic disease did not have TERT/ATRX alterations, but the late acquisition of TERT/ATRX alterations and presence in metastatic cells could not be excluded given that the metastatic tissue was not sequenced. However, metastasis in the absence of TERT/ATRX alterations was also evident in three metastatic tumours for two patients (E158, E159). These metastases had no other evidence of a TLM (i.e., they did not overexpress TERT, and had neither detectable C-circles nor an increased telomere length ratio) and had no recurrent cancer gene alterations that may explain metastatic progression. Interestingly, we found the tumour mutation burden was significantly elevated in TERT/ATRX-altered PCPG (one-sided Student’s t-test p < 0.001) consistent with a previous study¹¹, but metastases from the TERT/ATRX wild-type cases had a low mutation burden (Fig. 6C). Interestingly, while it is has been previously shown that a larger primary tumour is associated with metastatic risk and TERT/ATRX mutations¹¹, only patients with TERT altered tumours in our series had significantly larger primaries when contrasting to patients without TERT/ATRX alterations (two-sided students t-test, p < 0.05)(Fig. 6D). Finally, survival analysis showed a significant association with poorer outcome in the TERT/ATRX altered group, which is also consistent with prior studies^11,15 (Supplementary Table 1, Supplementary Fig. 16).

Transcriptional patterns in TERT/ATRX altered and metastatic PCPG

We next assessed transcriptional changes in TERT and ATRX altered PCPG, respectively, by contrasting each group independently to non-metastatic PCPG and directly to each other. Furthermore, we contrasted all metastatic and non-metastatic PCPG in a separate analysis (Fig. 7A). All HN-PG/non-chromaffin tumours (defined by transcriptional profile) were excluded given the absence of TERT/ATRX alterations in these tumours. Furthermore, we excluded low tumour purity cases ( < 50% estimated by WGS data) and a subclonal TERT-mutant case. Contrasting the TERT and ATRX altered tumours to non-metastatic tumours we found 1152 and 1448 differentially expressed genes, respectively (Supplementary Data 6). Some of these genes were common to both analyses (n = 273), as well as when contrasting all metastatic and non-metastatic tumours (n = 180).

**Fig. 7: Differential gene expression between *TERT* and *ATRX*-altered and non-metastatic tumours.**

Consistent with prior studies, genes overexpressed in metastatic PCPG were enriched for cell cycle and proliferation gene ontologies^11,22. MKI67 was among the top overexpressed genes in keeping with elevated Ki67 IHC staining (Fig. 7B). Other genes included cell cycle regulators (e.g., TOP2A, BUB1, FOXM1, AURK2B, CDK1), that were highly correlated with MKI67 expression (r > 0.88). Other genes overexpressed in metastatic PCPG included the polycomb repressor complex 2 gene EZH2, as previously described²², and the transcription factors OTX1 and TPX1.

Long non-coding RNAs and pseudogenes were enriched in genes differentially expressed in ATRX altered tumours (Pearson’s Chi-squared p-value = 9.9e-14)(Supplementary Fig. 17). The serine/threonine protein kinase RIPK4 was also overexpressed, while genes downregulated in ATRX altered tumours included the p53-dependent G2M cell cycle regulator RPRM and the GTP-binding gene DRG2 – the latter gene also implicated in G2M cell cycle checkpoint control (Fig. 7C)⁴². Other genes down-regulated in ATRX altered tumours included USE1, SULT4A1 and the nuclear-encoded mitochondrial gene COX17. Many differentially expressed genes were also within differentially DNA methylated regions (Fig. 7C). Gene-expression in neoplastic cells was confirmed using the snRNA-seq data and differential expression in ATRX mutant tumours was validated using independent bulk RNA-seq datasets^13,14 (Fig. 7D). Meanwhile, few genes were uniquely expressed in TERT altered tumours; amongst these were IRX3, SDK1 and TRIP13 (Fig. 7E).

Differential expression of small-RNA from TERT altered and ATRX altered tumours with non-metastatic primary tumours revealed down-regulation (adj. p < 0.05, log fold-change < −1.5) of four (hsa-miR-485-5p, hsa-miR-433-3p, hsa-miR-539-3p, hsa-miR-127-5p) and three (hsa-miR-3065-3p, hsa-miR-3065-5p, hsa-miR-148a-3p) miRNAs, respectively. In concordance with a previous study⁴³, three miRNA (hsa-miR-96-5p, hsa-miR-183-5p, hsa-miR-182-5p) were found to be upregulated (adj. p < 0.05, log fold-change > 1.5) in both TERT altered and ATRX altered tumours relative to non-metastatic primaries. These three miRNAs are collocated on chr7q32.2 and share targets including Forkhead Box O1 transcription factor (FOXO1), snail family transcriptional repressor 2 (SNAI2), and RECK, a membrane-anchored matrix metalloproteinase regulator and potential inhibitor of metastasis^44,45,46 (Supplementary Fig. 18, Supplementary Data 7).

To explore which transcription factors may be driving changes in transcriptional programming in TERT and ATRX altered tumours we performed binding motif enrichment analysis using single-nuclei ATAC-sequencing. When compared to normal chromaffin cells, TERT altered cells demonstrated an enrichment of binding sites associated with the Activator Protein-1 (AP-1) transcription factor family (JUNB, JUND, FOSL1, FOSL2), as well as the SWI/SNF related BAF chromatin remodelling complex subunit C1 (SMARCC1), and the Wilms’ tumour 1 (WT1) transcription factor. ATRX-mutant cells also showed an enrichment of peaks containing WT1 binding motifs, as well as the Sp1 and Sp2 transcription factors, however the greatest enrichment was seen for CTCF binding motifs. Interestingly, normal chromaffin cells showed an enrichment (i.e. depleted in TERT/ATRX altered cells) of the glucocorticoid (NR3C1) and androgen (AR) receptors (Supplementary Fig. 19).

Treatment-related mutagenesis and resistance to DNA alkylating chemotherapy

DNA alkylating chemotherapies (e.g., dacarbazine and temozolomide) are common treatments for metastatic PCPG⁴⁷. Two patients had paired tumours taken before and after cyclophosphamide, vincristine and dacarbazine (CVD) treatment. Case E169 had a para-aortic abdominal-thoracic PGL with distant lung, liver, and bone metastases (Fig. 8A). Tumours resected before (lung, E169-M1) and after (spinal disease, E169-M2) 23 cycles of CVD were used for genomic analysis. Paired E169 tumours had a high number of shared SNVs (n = 3134) including mutations in TERT and BCOR; however, the number of SNVs doubled in the post-CVD (n = 12607) tumour including mutations in TP53, RPL5 and POLE (DNA Polymerase Epsilon) (Fig. 8B). SBS signatures also changed in the post-treatment sample with increased representation of SBS8, SBS12, SBS16 and SBS21 as well as indel signature ID1 (Fig. 8C). However, elevated SBS11 associated with DNA alkylating chemotherapy and SBS14 associated with DNA polymerase epsilon-deficiency was absent. Importantly, MGMT was overexpressed in the post-treatment tumour (E169-M2) (Fig. 8D). MGMT overexpression is a known temozolomide resistance mechanism in glioblastoma cells⁴⁸ and therefore potentially an acquired dacarbazine resistance mechanism in PCPG.

**Fig. 8: Evolution of PCPG under treatment pressure.**

A second case (E167) had a left adrenal primary PC and then 14 years later developed HN-PGL and metastatic disease in the hip (Fig. 8E). The patient progressed with development of liver metastases after ¹³¹I-MIBG therapy and was subsequently treated with 22 cycles of CVD. Further progression occurred with new intracranial and osseous disease. Liver and intracranial metastases represented disease pre- and post-CVD treatment respectively, with clonal relatedness confirmed by common SNVs (n = 2404) including a pTERT mutation (Fig. 8F). MGMT expression was low in both tumours (Fig. 8D). The post-CVD metastasis acquired more than a million additional SNVs corresponding to a strong SBS11 signature (Fig. 8G). Interestingly, the post-treatment intracranial metastasis (E167-M2) had a somatic NRAS (NP_002515.1:p.G12A) mutation and predicted to have elevated MAP kinase pathway activity by WTS (Supplementary Fig. 20). Importantly, a protein truncating MLH1 splice-site mutation (chr3:g.37091976 G > A) was detected in the post-CVD metastasis with concordant loss-of-heterozygosity and low MLH1 expression (Fig. 8F, H). Canonical SBS mismatch-repair (MMR) mutational signatures (e.g. SBS6, SBS14, SBS21) were not detected in the post-CVD treated tumour; however, insertion/deletion signature ID2 was elevated and insertion/deletion events at microsatellite sites revealed a high MSI-score (11.7 microsatellite indels per megabase) in the post-CVD tumour (E167-M2), which is above the prescribed microsatellite instability threshold ( > 4)(Supplementary Fig. 21).

Pathogenic variants in MMR genes (MLH1, MSH2, MSH6, PMS2) as well as treatment-related hypermutation and SBS11 have been described in temozolomide-treated glioblastomas^49,50 and more recently also pancreatic neuroendocrine tumours⁵¹. To search for evidence of MMR-deficiency in additional PCPG we analysed variant data from AACR Project GENIE (release 13), representing comprehensive cancer panel data for 75 PCPG. We identified one hypermutated case with a dominant C > T mutational pattern and an MSH6 variant predicted to be a driver by Cancer Genome Interpreter (GENIE-DFCI-001077-11300, NP_000170.1:p.G162E) (Fig. 8I, J). Furthermore, we performed 523-gene panel sequencing of plasma cell-free (cf) DNA from a patient with SDHB-associated PCPG with high disease volume as shown by ¹⁸F-FDG-PET (Fig. 8K). Blood was taken after three cycles of temozolomide combined with ¹⁷⁷Lutetium-DOTA-octreotate. The patient had high fraction of circulating tumour (ct) DNA (45%) (Fig. 8L). Hypermutation was evident in cfDNA with a predominance of C > T transitions corresponding to SBS11 (Fig. 8J), although most C > T mutations had low VAF ( ~ 0.5%) in contrast to truncal mutations involving pTERT, MYCN, NUTM1 and MAP2K1 (VAF > 10%) (Fig. 8M). A mutation in an MMR gene was not detected; however, the hypermutated clone likely represented a small fraction of total disease volume therefore detection of an MMR gene mutation may have fallen below the detection limit of the cfDNA assay (VAF > 0.5%).

Discussion

We performed a comprehensive genomic analysis of SDHB-associated PCPG thereby creating a data resource to understand development of metastatic disease in these patients. Consistent with prior studies, most PCPG had stable genomes with few recurrently mutated cancer genes^11,13,14,20. Importantly, we confirm SDHB mutant PCPG are molecularly distinguishable from other PCPG genotypes, while parasympathetic (non-chromaffin) tumours including HN-PG had a very distinct molecular profile. Our data confirms TERT and ATRX alterations are the most common secondary driver events in PCPG and are enriched in metastatic disease^11,13,15. However, a small number of metastatic PCPG lack TERT/ATRX alterations, suggesting yet unidentified mechanisms of metastatic progression in these cases.

Mutual exclusivity of TERT/ATRX alterations in PCPG suggests these events have a common and potentially redundant role in metastasis. Both genes are involved in telomere maintenance - a hallmark of cancer⁵². However, given the majority of PCPG, including a subset of metastatic PCPG, have no obvious telomere lengthening mechanism (TLM), how, or indeed whether, these tumours maintain their telomeres is unknown. While methods such as the Telomere Repeat Amplification Protocol (TRAP)⁵³ and the C-circle assay are required to exclude the presence of a known TLM in tumour cells, absence of TLM has also been reported in other neuroepithelial cancers^54,55. Meanwhile, while ATRX/TERT alterations were strongly associated with metastatic progression, the presence of long telomeres in the absence of TLM appears unlikely to be sufficient to drive metastatic progression¹⁵.

ATRX is a SWI-SNF chromatin remodelling protein with multiple functions beyond telomeric regulation, while TERT is also thought to have non-canonical functions⁵⁶. We observed that increased cell cycle activity appears to be common to ATRX/TERT altered PCPG, which is consistent with prior studies^11,22 and in keeping with increased Ki67 staining within metastatic PCPG⁵⁷. Experimental models implicate ATRX and TERT in regulating cell cycle activity. For instance, ATRX binds to regulatory elements of cell cycle genes including CHEK1, as previously reported in a glioblastoma model⁵⁸. Intriguingly, we found repression of RPRM and DRG2 in ATRX altered PCPG, which are known to regulate the G2M cell cycle checkpoint^42,59. It remains to be determined if RPRM and DRG2 are ATRX target genes. TERT can also promote cellular proliferation and has a pro-survival function. For instance, cellular proliferation can be inhibited by wild-type reversion of a pTERT mutation in a glioblastoma cell line⁶⁰. Furthermore, catalytically inactive TERT can rescue TERT-null cells, indicating a pro-survival function that is independent of telomerase activity⁶¹. How TERT promotes cellular proliferation is poorly understood and may be cell type and context dependent. Prior studies have shown TERT occupancy at gene-promoters including tRNAs⁶² while TERT has also been shown to interact with oncogenic pathways involving MYC, WNT/B-catenin or NFkB⁵⁶. Given few differentially expressed genes were restricted to TERT altered PCPG, experimental models will be required to understand any non-canonical TERT functions in PCPG.

ATRX/TERT alterations have potential implications for treatment and surveillance. An increased sensitivity of ATRX altered tumour cells to DNA damage response (DDR) inhibitors targeting ATM⁵⁸, WEE1⁶³, ATR and PARP^64,65 has been reported in other neuroectodermal cancers, and therefore might be effective in PCPG. Sensitivity to radiosensitizing DDR-targeting drugs in combination with radionuclide therapies (e.g., ¹⁷⁷Lutetium-DOTA-octreotate or ¹³¹I-MIBG) also warrants investigation in the context of ATRX/TERT alterations^66,67. Although we confirmed TERT/ATRX alterations are associated with metastatic disease, the utility of these mutations as predictive biomarkers may be limited by the late timing of these events. Detection of TERT/ATRX alterations by cfDNA analysis may be a viable alternative to tissue biopsies, although the broader utility of cfDNA detection in PCPG still requires further validation.

DNA alkylating agents are commonly used for PCPG but there are currently no validated biomarkers to predict treatment response. We identified hypermutation together with MMR mutations (e.g., MLH1, MSH6) or MGMT overexpression in dacarbazine or temozolomide-treated PCPG. These observations are potentially clinically important. Firstly, monitoring PCPG patients for development of MMR mutations, SBS11 or potentially MGMT overexpression may indicate development of treatment-resistant cells. Secondly, MMR mutations may sensitize tumour cells to other treatments. For instance, hypermutation can theoretically generate neoantigens causing T-cell recognition in tumours, making them more immunogenic and potentially more sensitive to immune checkpoint inhibition (ICI), although this theory is tempered by observations in glioblastoma where temozolomide-induced hypermutation did not translate to ICI-response⁵⁰. Synthetic lethality could be another approach in MMR-deficient tumours, including use of WRN helicase inhibitors⁶⁸. Further investigation is required to determine whether these treatments are efficacious in MMR-deficient PCPG.

Although therapeutically actionable mutations are rare in SDHB mutant PCPG, co-operative drivers can still represent potential treatment opportunities in some patients. In particular, somatic EPAS1 variants in SDHB-mutant tumours may have translational relevance. It is currently unknown whether HIF2α inhibitors (e.g., belzutifan) are effective for SDHB mutant PCPG, but it is plausible that co-operative EPAS1 mutations could make these tumours more sensitive⁶⁹. Screening for co-operative somatic EPAS1 mutations in SDHB-mutation carriers should therefore be employed for future HIF2α inhibitor clinical trials. Although we have reported a comprehensive analysis of SDHB mutant PCPG, future studies will be needed to validate our findings and identify additional low frequency genomic features that may be used as clinical biomarkers and future treatment targets.

Methods

Patient samples and human research ethics

The research was conducted under an approved protocol of the human research ethics committee at Peter MacCallum Cancer Centre under the guidelines of the National Health and Medical Research Council and in accordance with the Helsinki Declaration of 1975, as revised in 1983. Patients were recruited at 11 sites under protocols approved by their respective institutional review boards with written informed consent without any compensation. Patient recruitment sites included those forming the American-Australian-Asian-Adrenal-Alliance (A5) International Research Consortium including the Peter MacCallum Cancer Centre (Australia), the Kolling Institute Neuroendocrine Tumour Bank under a protocol approved at North Sydney Local Health District (Australia), National Institute of Health (USA), University of Colorado (USA), University of Texas Health Science Center at San Antonio (USA), University of Florida (USA), University of Michigan (USA), Centre hospitalier de l’université de Montréal (Canada), as well as four non-A5 sites at Tufts Medical Centre (USA), Waikato Hospital (New Zealand), the National Cancer Centre (Singapore), and Uppsala University (Sweden). Patient sex was determined by genomic profiling. The findings of this study apply equally to both sexes.

Tissue sample review and nucleic acid extraction

Tissue samples were assessed for tumour cell content by a single pathologist (AG) by haematoxylin and eosin staining of 4uM cryosections. Ten serial 10uM sections were taken for nucleic acid extraction adjacent to the H&E reviewed sections. If required, microdissection of tissue was performed on ethanol fixed and stained fresh tissue sections to dissect tumour regions from the sectioned material. DNA was extracted from tissue samples and matched patient whole blood using the DNeasy® Blood & Tissue kit (Qiagen, Germany) according to manufacturer’s protocol. Total RNA was extracted from tissue using the miRNeasy mini kit (Qiagen, Germany) according to the manufacturer’s protocol. Nucleic acids were quality control checked by Agilent TapeStation microfluidic electrophoresis (Agilent, USA) and NanoDrop spectrophotometer (Thermo Fisher Scientific, USA).

Targeted amplicon sequencing

To validate mutations of interest, targeted amplicon sequencing was done using custom primers compatible with universal indexed primers for Illumina cluster generation. PCR amplification for the region of interest was done in a single well (for each patient sample) and the expected PCR amplicon product then verified by agarose gel electrophoresis. Multiple amplicons for each patient sample were pooled, purified, and then barcoded with the index primers, quantified by fluorometer (Qubit dsDNA High-sensitivity, USA) and then sequenced on a MiSeq using a V3 300 cycle kit (Illumina, USA). FASTQ files were aligned to human genome (GRCh38) using BWA-mem and BAM files were visualized in IGV⁷⁰ to confirm variant of interest and record the variant allele frequency.

Whole genome sequencing

Library preparation

Libraries were prepared using the Illumina® TruSeq™ DNA Nano library preparation method according to the manufacturer’s instructions at The University of Melbourne Centre for Cancer Research (UMCCR) using 200 ng input DNA and a 550 base pair insert size. Samples were sequenced in separate batches on the Illumina® Nova-Seq 6000 according to manufacturer’s instructions (Illumina, USA).

WGS alignment and small variant calling

BCL files were demultiplexed and converted to FASTQ files with Illumina® bcl2fastq (version 2.20.0.422) with default settings (including adapter trimming). Alignment and variant calling steps were run using the bcbio-nextgen cancer somatic variant calling pipeline (version 1.2.4-76d5c4ba) (https://github.com/bcbio/bcbio-nextgen). The versions of all programs used by bcbio-nextgen are additionally shown in Supplementary Table 2. Briefly, the tumour and normal sequencing data in FASTQ format were first processed by Atropos (version 1.1.28)⁷¹ to clip homopolymers (minimum 8 bases) from the ends of reads. Trimmed FASTQ files were aligned to the human genome (GRCh38) with BWA-mem (version 0.7.17)⁷² with predominantly default settings. The two exceptions were that the -M flag was enabled to mark shorter split read hits as a secondary alignment and seeds with > 250 occurrences were skipped (reduced from the default 500).

Germline variants were called by vardict (version 1.8.2)⁷³, GATK HaplotypeCaller (GATK version 4.1.8.1)^74,75 and Strelka2 (version 2.9.10)⁷⁶ variant callers. Somatic variants were called by vardict, Mutect2 (see GATK version)⁷⁷ and Strelka2 (version 2.9.10) variant callers. In both instances a variant was accepted as true if it was detected by 2 of the 3 callers. Somatic variants were further post-processed by umccrise (version 1.2.0-rc.1, https://github.com/umccr/umccrise). To further reduce spurious variant calls, two variant blacklists were created. The first was generated by counting supporting reads for all detected somatic variants across all normal BAM files using SAMtools⁷⁸ (version 1.12) mpileup. Any somatic variants that were supported by three or more reads in more than two normal controls were excluded. The second blacklist was created by extracting all observations of pass-filter somatic variants from unfiltered variant call files (i.e., including those that were rejected by the caller) for all three somatic callers. The number of observations of a given variant across all samples and callers were grouped by pass or reject status and summed. Any variant with a reject to pass ratio of greater than one was placed on the blacklist.

WGS copy number and structural variant calling

Structural variants (SV) were called with GRIDSS (version 2.8.0)⁷⁹. GRIDSS output was annotated for repeat sequences using RepeatMasker (version 4.1.5) and filtered for somatic SVs using GRIPSS (version 1.11). To ensure high quality somatic calls, all SVs passing GRIPSS filtering were examined in the raw output of GRIDSS across all samples. Any SV (based on exact genomic coordinates) that was observed (whether ultimately rejected or passed) in more than two unrelated samples was removed. Finally, SVs were required to have either a minimum of three supporting split reads or paired reads, or a combination of at least one split and one paired read. Somatic copy number profiles were generated using PURPLE (version 3.1, https://github.com/hartwigmedical/hmftools) from input generated by COBALT (version 1.11), AMBER (version 3.5), and GRIDSS. Structural variants were further annotated with Linx (version 1.16) using the output from PURPLE.

Telomere length estimation and telomere variant repeat detection

Telomere content was estimated from WGS for each tumour and normal pair using TelomereHunter (version 1.1.0)⁸⁰ with default parameters. Telomere content was computed from the number of intratelomeric reads normalized by the total number of reads with a GC composition like that of telomeres. Telomeric enrichment in tumours was defined as a log₂ ratio of tumour telomere content to blood telomere content greater than 0.5. TelomereHunter was also applied to WTS for each tumour to estimate TERRA content.

Mutational signature analysis

For each sample, the single base substitution (SBS), doublet base substitution (DBS), and small insertions and deletions (ID) mutation count matrices were extracted from variant call format (VCF) files containing somatic variant calls using the MutationalPatterns (version 3.10.0)⁸¹ package for R. The optimal linear combination of the COSMIC V3⁸² SBS, DBS, and ID mutational signatures that best reconstructed their respective profiles was then estimated using the fit_to_signatures function.

Microsatellite instability assessment

Insertions and deletions in microsatellite regions were assessed using PURPLE (version 3.1). Samples with more than four microsatellite Indels per megabase were considered to have microsatellite instability as per the software default setting.

Identification of recurrent copy number alterations

The Genomic Identification of Significant Targets in Cancer (GISTIC, version 2.0)⁸³ algorithm was used to identify recurrent copy alterations. The input data for GISTIC was modified from the somatic segmentation output from PURPLE by taking log₂ of the segmental copy number incremented by a small offset (0.01) and subtracting one to centre the diploid state around zero (log₂(CN + 0.01)-1). GISTIC was run in broad mode with segments containing less than 50 markers undergoing joining. The broad length cutoff was set to 75% of a chromosome and arm-peel was activated. The analysis was applied independently to sample subsets including HN-PG/parasympathetic PG, non-metastatic and metastatic sympathetic PCPG. In cases where a patient had more than one primary tumour, one was chosen at random, in cases where a paired primary and metastasis were present, the metastasis was used.

Identification of driver genes

The dNdScv R package³³ was used to assess somatic mutations to detect genes under positive selection. Before running the analysis, quality filtered somatic variant calls were lifted over from the GRCh38 to HG19 genome build using the rtracklayer package for R with a liftover chain obtained from UCSC (hg38ToHg19.over.chain). The resulting mutation set was analysed with dndscv function using default parameters. Additionally, the HG19 mutation set was used to run the OncodriveFML³⁴ driver detection tool via the web-based application interface (https://bbglab.irbbarcelona.org/oncodrivefml/home, accessed 5th February 2024). The algorithm was run using genomic region subsets encompassing coding sequences or promotor regions with Combined Annotation Dependent Depletion (CADD) scores version 1.0.

Clonal evolution analysis

For each patient with paired samples, variants detected in any sample from a given patient were recalled using SAMtools (version 1.12) mpileup in all samples for that patient. The resulting allele counts were then used as input for Pairtree (version 1.0.1)⁸⁴. Variants that were detected in between zero and three reads in any sample were marked for exclusion from all samples for a given patient as it was difficult to discern sub-clonality from sequencing artifacts. Variant read probabilities (VRP) were computed based on copy number segment data from PURPLE. In diploid regions the VRP was set to 0.5. In regions of single copy loss or copy-neutral loss of heterozygosity where the variant allele fraction was greater than half the sample purity the VRP was set to 1. In regions of copy-neutral loss of heterozygosity where the variant allele fraction was less than half the sample purity the VRP was set to 1 divided by the major allele copy number. For all other region types the VRP was set to 1 divided by the total copy number. Clustering was performed using the clustervars script. The resulting clusters were manually refined and integrated with copy number data to produce evolutionary trees.

RNA-seq

Library preparation

Samples were prepared using the NEB-Next directional RNA-Seq kit (NEB, USA) and underwent 150 bp paired-end sequencing on the Illumina NovaSeq 6000 (Illumina, USA) according to manufacturer’s instructions.

Read quantification

RNA sequencing reads were aligned with STAR⁸⁵ (version 2.7.9a) using a two pass approach. First, all data was aligned using a STAR reference built with the Gencode human reference genome with decoy HLA sequences (GRCh38) and Gencode gene transcript annotation file (version 36, primary assembly annotation). Following alignment, the splice junction output from all samples was merged. Supporting read evidence for each junction was summed and junctions with less than 10 unique supporting reads were removed. The resulting junction file was used to build a new STAR reference and the alignment was repeated. Gene level feature counts were obtained with HTSeq-count (version 0.11.3) using the previously described Gencode transcript annotation file. All parameters were default except Mode which was set to intersection-nonempty.

Differential gene expression

To contrast gene expression between parasympathetic and sympathetic (including pheochromocytoma) paraganglioma, samples were assigned to one of two groups based on clinical annotation and UMAP clustering. Samples where UMAP clustering did not reflect the clinical annotation or tumour purity was below 50% were excluded from the contrast. To perform differential expression analysis, transcript counts obtained from HTSeq-count were converted to a digital gene expression object using the edgeR (version 3.42.4)⁸⁶ package for R. Lowly expressed genes were removed using the filterByExpr function with default parameters. Library scaling factors were then computed using the calcNormFactors function. A design matrix was constructed without an intercept to include the differential group and patient sex as a covariate. Normalization factors were applied and correction weights for count-dependent variance were computed using the voom function from the limma⁸⁷ (version 3.56.2) package for R. To account for multiple samples from the same patient, the inter-duplicate correlation was computed using the duplicateCorrelation function with the patient identifier as the blocking factor. Linear model fitting was performed using the lmFit function supplied with the duplicate correlation computed above and the patient identifier as the blocking variable. Contrasts were fitted using the contrasts.fit function and test-statistics were computed using the eBayes function.

To contrast metastatic, TERT-altered, and ATRX-altered tumours with non-metastatic primary tumours samples were assigned to groups based on clinical behaviour. The clinical behaviour groups were non-metastatic primary tumour, primary tumour where metastatic disease was present, metastatic lesion, or a tumour with short clinical follow up. These groups were then further segregated based on the presence or absence of a TERT or ATRX mutation. Parasympathetic tumours and tumours with low purity ( < 50%) were excluded from all the contrasts. Samples with an ATRX/TERT alterations where the allele fraction was less than 25% were excluded from genotype-centric contrasts, and samples with short clinical follow-up were excluded from contrasts involving clinical behaviour. See Supplementary Data 2 for a full list of contrast participation. Differential expression was performed as described above.

RNA fusion detection

Whole transcriptome sequencing data aligned with STAR aligner (described above) was analysed with the Arriba (version 2.1.0) gene fusion detection tool using the supplied blacklist, known fusion, and protein domain files with default parameters.

Harmonization with published RNA-seq datasets

Read counts previously published by Fishbein et al.¹³ were obtained from the Genomic Data Commons using the TCGAbiolinks⁸⁸ (version 2.29.1) package for R with these parameters: project = “TCGA-PCPG”, data.category = “Transcriptome Profiling”, data.type = “Gene Expression Quantification”, workflow.type = “STAR - Counts”. Additional data previously published by Flynn et al.¹⁴ was obtained from the European Genome-Phenome Archive (EGAS00001005861). Raw counts from the publicly available data were merged with our dataset retaining genes present in all three datasets. Count matrices were converted to a digital gene expression object with edgeR and lowly expressed genes were removed using the filterByExpr function with default parameters. Library scaling factors were then computed using the calcNormFactors function and counts were transformed to log counts per million using the cpm function. Possible technical variance was ameliorated using the removeBatchEffect function from the limma package where the dataset identifier was supplied as the batch and the primary driver genotype (e.g., RET, SHDx, VHL) as the factor to preserve.

Small RNA-seq

Library preparation

Samples were prepared using the NEXTFLEX® Small RNA-Seq Kit v3 (Bioo Scientific). Samples were sequenced at the molecular genomics core (Peter MacCallum Cancer Centre) using 50 bp single end sequencing on the Illumina NextSeq 500 (Illumina, USA).

Read quantification

FASTQ files were pre-processed with two rounds of cutadapt (version 1.6.3) to first remove adapter sequences (-a TGGAATTCTCGGGTGCCAAGG, --minimum-length 23) and then remove the first and last four bases (-u 4, -u -4). Following trimming, sequences were aligned to GRCh38 with subread (version 1.6.3) and quantified using featurecounts (version 1.6.3) in stranded mode against the miRBase (version 20) miRNA transcripts⁸⁹. Following read quantification, read-count processing was performed in an identical manner to that described under RNA-Seq analysis.

Harmonization with published small-RNA-seq datasets

Small-RNA sequencing data previously published by Castro-Vega et al.²⁰ were obtained from ArrayExpress (E-MTAB-2833). Additional data published by Fishbein et al.¹³ were downloaded from the Genomic Data Commons using the GDC Data Transfer Tool in BAM format and converted to FASTQ format using SAMtools (1.10). Alignment and read quantification were performed as described above under ‘Read quantification’. Read count harmonization was performed as described for RNA-seq dataset under ‘Harmonization with published RNA-seq datasets’.

Differential expression of small-RNA dataset

Differential expression analysis was performed as outlined for whole transcriptome RNA-sequencing. The same contrast groups were used, however, samples that were determined to be strong outliers by UMAP clustering of small-RNA expression data were excluded from the analysis. To ensure biological relevance for miRNA activity, genes were required to have at least 100 counts per million in more than one sample and a total of more than 1000 counts per million across all samples to be included in the analysis.

Methylation analysis

Methylation array processing

Methylation profiling was performed by the Australian Genome Research Facility (AGRF). Tumour DNA samples were normalized to approximately 200 ng of DNA in 45 µL. Bisulfite conversion was performed with the EZ-96 DNA Methylation MagPrep (Zymo Research, California, USA). Bisulfite converted DNA was processed on the Illumina Human Methylation EPIC BeadChip (Illumina, California, USA) following the Illumina Infinium HD Methylation assay protocol (AGRF NATA scope GGTMN00263).

Methylation array data pre-processing and normalization

The minfi (version 1.46.0) package for R was used to read in a targets manifest and red/green channel IDAT files using the read.metharray.sheet and read.metharray.exp functions, respectively. Probe detection p-values were determined with the detectionP function. Raw green/red values were converted to methylation signal and functionally normalized using the preprocessFunnorm function from minfi. Any probes with a detection p-value below 0.01 in any sample or those that fell on a SNP position (using dropLociWithSnps) were removed from the dataset. Furthermore, the maxprobes package for R was used to exclude any cross-reactive probes. Illumina probe annotation version 1.0 with GRCh38 coordinate information was obtained from: https://github.com/achilleasNP/IlluminaHumanMethylationEPICanno.ilm10b5.hg38.

Probe level differential methylation analysis

Probe level differential methylation analysis was performed using the Remove Unwanted Variation workflow with the RUV package via the MissMethyl (version 1.34.0) package for R^90,91. Firstly, functionally normalized methylation values were converted to M-values using the getM function from the minfi package and negative control probes were obtained using the MissMethyl’s getINCs function. For each contrast, two rounds of variation adjustment were performed. In the first round, an RUV model was computed with RUV-inverse algorithm using the RUVfit function applied using the Illumina negative control probes as controls. The rescaled variances were then applied using the RUVadj function and p-values for differential methylation of each probe were obtained with the topRUV function. Probes with a p-value greater than 0.5 were selected for use as controls in the second round of variation adjustment using the RUVfit function and the resulting model was applied using the RUVadj function.

Differentially methylated regions

Identification of differentially methylated regions was performed with the DMRcate⁹² package (version 2.14.1) for R. Firstly, functionally normalized methylation values were converted to M-values using the getM function from the minfi package. Contrast and design matrices were generated using the same process described for differential gene expression. For each contrast, individual probes were then annotated using the cpg.annotate function with analysis type set to differential and the default annotation set for Illumina EPIC arrays (ilm10b2.hg19). Differentially methylated regions were then identified using the dmrcate function with a lambda of 1000 and bandwidth scaling factor of 2. Significant DMRs were those with a Fisher multiple comparison statistic below 0.05.

Harmonization with previously published methylation array datasets

Illumina Infinium HumanMethylation450 methylation array data previously made available by Fishbein et al.¹³ were obtained from the Genomic Data Commons using the TCGAbiolinks⁸⁸ (version 2.29.1) package for R with these parameters: project = “TCGA-PCPG”, data.category = “DNA Methylation”, data.format = “IDAT”, legacy = False. The array data were processed as described under ‘Methylation array data pre-processing and normalization’ with the exception of the annotation file, which was substituted with the supplied Illumina Infinium HumanMethylation450 annotation. Additional data made available by Letouzé et al.⁹ were obtained from the Gene Expression Omnibus using the GEOquery package for R (accession GSE39198 and GSE43293). The GSE39198/GSE43293 data was only available as beta values and was not reprocessed. The beta-values from all three datasets were converted to m-values and merged, retaining only probes available on the Illumina Infinium HumanMethylation27 array design. A design matrix was generated with the model.matrix function from the limma package using the primary driver genotype (e.g., RET, SHDx, VHL, Unknown). Batch effects were then adjusted using the removeBatchEffect function from limma where the dataset was supplied as the first batch annotation and patient gender as the second batch annotation.

Single nuclei RNA-sequencing

Library preparation with 10x chromium

snRNA-seq was performed using the ‘Frankenstein’ protocol (https://doi.org/10.17504/protocols.io.bqxymxpw). Briefly, nuclei were extracted from frozen tissues and subjected to fluorescence-activated nuclei sorting (FANS) using 4′,6-diamidino-2-phenylindole (DAPI). Sorting was performed using a BD FACSaria 2 instrument, sorting between 3000 and 10,000 nuclei per sample, capturing both diploid and tetraploid populations. FAN-sorted nuclei were immediately processed using either the 10x Chromium Single Cell 5’ (PN-1000006, 4 samples) or 3’ (PN-1000075, 4 samples) Library & Gel Bead Kit (10x Genomics, USA). Once processed, snRNA-seq libraries were sequenced on the Illumina Nova-Seq 6000 (Illumina, USA) using 150 bp paired-end sequencing.

Data processing and analysis

Additional snRNA-seq data for two normal adrenal medulla samples (E240, E243) previously published by Zethoven et al.²² were obtained from the European-genome-phenome archive (accession EGAS00001005861). Single-nuclei RNA sequencing reads in FASTQ format were aligned and counted using CellRanger (version 7.2.0). Read-one data for libraries processed with the 3’ protocol was restricted to the first 26 bases using the r1-length flag. Feature counting was performed using the 10x supplied GRCh38 transcriptome (version refdata-gex-GRCh38-2020-A) with reads aligning to introns contributing to feature counts. Single-Cell Remover of Doublets (SCRUBLET, version 0.2.3) was used to identify gel-beads containing more than one cell. CellRanger count matrices were loaded using the Read10X and CreateSeuratObject functions from the Seurat (version 4.3.0.1) package for R in conjunction with the SCRUBLET output. For each sample, features present in less than three cells and cells containing less than 500 features or 750 unique RNA molecules, or more than 10 percent mitochondrial RNA were excluded. Cell cycle scoring was performed on normalized and scaled counts using the S-phase and G2M-phase marker gene lists provided with the Seurat package. The data was scaled again using ScaleData function with S-phase score, G2M-phase score, unique RNA molecule count, and percent mitochondrial RNA supplied as regression variables. Cell type clusters were found using the FindClusters function and cell types were assigned manually using marker genes for fibroblast (FAP, PDGFRB, PDGFRA, ACTA2, COL1A1), chromaffin (PNMT, TH, DBH, CHGA, CHGB), adrenocortical (STAR, CYP11B1, CYP11A1), endothelial (FLT1, EPAS1), Schwann-cell-like cells (CDH19, SOX10, S100B, VIM), monocyte/macrophage (MSR1, CD163, CCL3), and lymphocyte (CD2, CD3E, MS4A1) cell types.

Single nuclei ATAC-sequencing

Library preparation with 10x Chromium

snATAC-seq was conducted using the “Van Helsing” protocol (https://doi.org/10.17504/protocols.io.bw52pg8e)⁹³. Initially, nuclei were isolated following the steps outlined in the previously mentioned “Frankenstein” protocol (https://doi.org/10.17504/protocols.io.bqxymxpw)⁹⁴. Sorted nuclei were collected into a round-bottom 96-well plate containing 100 µL of ice-cold Wash Buffer. The sorted nuclei were then transferred to 0.2 mL PCR tubes; 50 µL of ATAC Wash Buffer-Dig was added to the well to ensure transfer of any remaining nuclei. The nuclei were then centrifuged at 500 x g for 5 min at 4 °C. The supernatant was carefully decanted, leaving behind approximately 10 µL. Without resuspending the pellet, 100 µL of ice-cold Diluted Nuclei Buffer was added. The sample was centrifuged again under the same conditions, and 100 µL of the supernatant was gently removed in two steps—first 90 µL, then the remaining 10 µL. Another 100 µL of Diluted Nuclei Buffer was added gently. This washing step was repeated, removing the supernatant in two steps to leave behind approximately 7–10 µL. The nuclei were resuspended in this small volume, ensuring thorough mixing by carefully washing the walls of the tube. For quantification, 1–2 µL of the nuclei suspension was diluted 1:5 with Diluted Nuclei Buffer and mixed 1:1 with Trypan Blue. Nuclei were counted using a Countess II FL Automated Cell Counter or a hematocytometer to estimate the expected number of nuclei based on the recovery factor. Finally, 5 µL of the nuclei suspension in Diluted Nuclei Buffer were used directly in the Chromium Single Cell ATAC Reagent Kits protocol (CG000168 Rev A), following the specific recommendations for low input samples to optimize the transposition reaction and subsequent library preparation. Once processed, snATAC-seq libraries were sequenced on the Illumina NextSeq 500 (Illumina, USA) using 50 bp paired-end, duel-indexed sequencing.

Data preprocessing

Raw sequencing data was demultiplexed using cellranger-atac mkfastq (V2.0.0). Cellranger-atac count was used to combine replicate sequencing runs, align single nuclei ATAC reads to the GRCh38 reference genome (GRCh38 2020 A 2.0.0, 10x Genomics), and call initial peak sets for individual samples. Cellranger-atac aggr was then used to combine samples into a single peak-barcode matrix.

Quality control

Signac⁹⁵ (version 1.12.9004) and Seurat (version 4.3.0.1) were used for quality control and downstream analysis of snATAC-seq data. Filters were applied to QC metrics that were independent of peaks, to avoid bias toward under-represented cell types. Transcription start site enrichment scores were calculated using the Signac TSSEnrichment function. Outlier nuclei with over 5000 unique reads were assumed to be doublets or nuclei clumps and removed. Nuclei with TSS enrichment scores below 2 were removed due to having a low signal-to-noise ratio.

Cell type identification

Cell type identification was performed via label transfer in Signac. A gene-by-cell accessibility matrix was computed by calculating the sum of reads aligning to the gene body and promoter regions of protein-coding genes for each cell. Log normalisation was performed on the gene activities with Seurat NormalizeData. Integration was performed to link snATAC-seq nuclei clusters to defined cell types within the snRNA-seq dataset generated from matched and similar samples⁹⁶. Integration was performed with the FindTransferAnchors Seurat function. Briefly, Canonical Correlation Analysis (CCA) was used to identify shared sources of variation between the snATAC gene activity and snRNA gene expression, and project the two datasets into the same reduced-dimensional space. Mutual nearest neighbours (cell pairs from each dataset within each other’s neighbourhood) were identified and used as integration anchors between the two datasets. These anchors were then used to transfer cell-type labels from the snRNA reference dataset to the snATAC-seq dataset, using the Seurat TransferData function, enabling prediction of snATAC-seq nuclei cell types. Cell type identification was confirmed by gene activity of canonical cell-type-specific marker genes.

Peak calling and visualisation

Open genomic regions were defined using pooled reads for each cell type and tumour sample with > 100 total nuclei individually using the MACS2⁹⁷ peak caller (version 2.2.9.1) implemented in the Signac CallPeaks function. Resultant fragment pileup bedgraph files for each cell type were visualized using the SparK command line tool⁹⁸.

Motif enrichment analysis from sn-RNA-seq data

Motif analysis was performed using ArchR (version 1.02). Fragment counts output by Cellranger-atac were converted to arrow format using the createArrowFiles function (minTSS = 4, minFrags = 1000). Doublet scoring was performed using the addDoubletScores function. As some samples were reported to be homotypic, additional doublet detection was performed using AMULET (version 1.1). Doublets were then marked using a modified filterDoublets function which applied the standard ArchR filters and removed any cells considered to be doublets by AMULET. Cell type annotation was then transferred from the single nuclei RNA-seq data using the addGeneIntegrationMatrix function. A subset of the data containing only tumour and normal chromaffin cells was used to generate pseudo-bulk replicates using the addGroupCoverages (minRep = 3, maxRep = 10, minCells = 100, maxCells = 500, sampleRatio = 0.5) and addReproduciblePeakSet functions. The cell groups were defined as Tumour-TERT-mutant, Tumour-ATRX-mutant, and Normal-chromaffin-wild-type. Marker features were identified with the getMarkerFeatures function using the Wilcoxon test and applying bias adjustment for transcription start site enrichment scores and number of fragments in each cell. Enriched motifs in each contrast were determined using the peakAnnoEnrichment function with a false discovery rate cutoff of less than or equal to 0.05 and an absolute log fold-change greater than or equal to one.

Cell-Free DNA panel sequencing

Sample collection

Patient blood was collected in Cell Free DNA BCT CE tubes (catalogue number 218996, Streck, USA). Within two hours of collection, whole blood was centrifuged at 1900 g (brake off) for ten minutes to separate the plasma. Plasma was aspirated and aliquots were centrifuged at 20000 g to pellet residual cells or debris. The plasma was stored at −80 °C until DNA extraction. Cell free DNA (cfDNA) was extracted from 5 mL of plasma using the Qiagen QIAamp circulating nucleic acid kit (catalogue number 55114, Qiagen, The Netherlands) according to manufacturer’s instructions.

Library preparation

A sequencing library was generated with 30 ng of input cfDNA using the Illumina TruSightOncology 500 ctDNA library preparation kit (catalogue number 20039252, Illumina, USA). Input quantity was determined based on the mononucleosome fraction (75-300 bp) using the Agilent TapeStation4200. Library preparation was performed as per manufacturer’s instructions with the addition of a 1.2x bead clean-up of the amplified enriched library. The fragment distribution of the enriched library was assessed on the TapeStation4200, and library quantification was performed using the KAPA Illumina library quantification kit (catalogue number 07960140001, Roche, Switzerland) and analysed with QuantStudio6 (ThermoFisher, USA).

Sequencing and data processing

Sequencing was performed at the University of Melbourne Centre for Cancer Research Centre using the Illumina NovaSeq 6000 (Illumina, USA) to achieve a depth of coverage greater than 35000x. Data was analysed with Illumina DRAGEN TruSight Oncology 500 ctDNA Analysis Software v1.2 on the Illumina Connected Analytics platform (ICA v1). PierianDx Clinical Genomics Workspace (v6.27.0.1) was used for variant interpretation.

Immunohistochemistry

Immunohistochemistry was performed and interpreted on formalin fixed paraffin embedded tissue sections using methods that have been described in detail elsewhere^99,100. Briefly, for both SDHB (ABCAM, Cambridge, UK: ab14714, clone 21A11, dilution 1:100) and Ki67 (clone M7240 Dako, Carpenteria, CA USA, dilution 1:50) commercially available mouse monoclonal antibodies were used.

C-circle assay

Isothermal rolling circle amplification and detection of C-circles was performed as described^40,101. For each tumour sample, 20 ng DNA was incubated with 10 μl 0.2 mg/ml bovine serum albumin, 0.1% Tween, 4 mM DTT, 1 mM each dATP, dGTP, dCTP and dTTP, 1× Φ29 Buffer, with or without Φ29 polymerase 7.5 U (New England Biolabs) at 30 °C for 8 hours, then at 65 °C for 20 min. Reaction products were diluted to 100 μl with 2× SSC and dot-blotted onto a 2× SSC-soaked Biodyne B nylon membrane (Pall Corporation). The membrane was air-dried, twice crosslinked using 1200 J UV (Stratagene) and hybridized overnight at 37 °C with end-labelled 32P-(CCCTAA) and PerfectHyb Plus hybridization buffer (Sigma). The membrane was then washed at 37 °C three times in 0.5x SSC/0.1% SDS and exposed to a storage phosphor screen, which was imaged using an Amersham Typhoon laser scanner (Cytiva) and analysed by densitometry using ImageQuant software (Molecular Dynamics). C-Circle levels for each tumour sample were reported as the difference in the dot-blot intensities resulting from reaction of the sample’s DNA with and without Φ29 polymerase.

Statistics & reproducibility

Student’s t-tests and Pearson’s Chi-squared tests for count data were performed using the t.test and chisq.test functions, respectively, from the base stats package for R (version 4.3.1). No statistical method was used to predetermine sample size. Samples and assays failing quality control measures for tumour content, nucleic acid quality, or assay specific controls were excluded. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

Raw FASTQ files for WGS, WTS, small-RNA-sequencing, sn-ATAC-sequencing, sn-RNA-sequencing, as well as GRCh38 aligned CRAM files for WGS (accession EGAS50000000346), and IDAT files for Illumina EPIC methylation array profiling (accession EGAS00001007844) are available from the European-genome-phenome (EGA) archive. As genomic data represent potentially identifying information, these data are under restricted access controlled by the Data Access Committee at the University of Melbourne Centre for Cancer Research. Access to sequence read-level data or germline variant data will require evidence of institutional human research ethics committee approval. Access can be obtained by application through the EGA portal and applications will be processed within four to six weeks. If approved, data access will be granted in perpetuity. Unrestricted access to a data package containing non-identifiable data from tertiary analysis including somatic small and structural variant calls, copy-number analysis, RNA-fusion calls, and gene expression matrices from bulk and single nuclei analysis is available from Figshare (https://doi.org/10.6084/m9.figshare.25792479). RNA-seq, Illumina Infinium HumanMethylation450 arrays, and small-RNA-seq data previously made available by Fishbein et al.¹³ was obtained from the National Computational Infrastructure’s Genomic Data Commons. Bulk RNA-seq data previously published by Flynn et al.¹⁴ and snRNA data published by Zethoven et al.²² was obtained from EGA (EGAS00001005861 [https://ega-archive.org/studies/EGAS00001005861]). Small-RNA-seq previously published by Castro-Vega et al. were obtained from ArrayExpress (E-MTAB-2833 [https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2833])²⁰. Illumina Infinium HumanMethylation450 and HumanMethylation27 arrays previously made available by Letouzé et al. were obtained from the Gene Expression Omnibus (accession GSE39198, and GSE43293)⁹. The remaining data are available within the Article, Supplementary Information, or Source Data file. Source data are provided with this paper.

Code availability

All code used in the generation and analysis of data contained within this publication is available from GitHub: https://github.com/UMCCR-RADIO-Lab/a5_sdhb_pcpg (https://doi.org/10.5281/zenodo.14890696)¹⁰².

References

Nazari, M. A. et al. Catecholamine-induced hypertensive crises: current insights and management. Lancet Diab. Endocrinol. 11, 942–954 (2023).
Article CAS MATH Google Scholar
Lenders, J. W., Eisenhofer, G., Mannelli, M. & Pacak, K. Phaeochromocytoma. Lancet 366, 665–675 (2005).
Article PubMed Google Scholar
Ayala-Ramirez, M. et al. Clinical risk factors for malignancy and overall survival in patients with pheochromocytomas and sympathetic paragangliomas: primary tumor size and primary tumor location as prognostic indicators. J. Clin. Endocrinol. Metab. 96, 717–725 (2011).
Article CAS PubMed Google Scholar
Lenders, J. W. M. et al. Genetics, diagnosis, management and future directions of research of phaeochromocytoma and paraganglioma: a position statement and consensus of the Working Group on Endocrine Hypertension of the European Society of Hypertension. J. Hypertens.38, 1443–1456 (2020).
Article CAS PubMed PubMed Central MATH Google Scholar
Amar, L. et al. International consensus on initial screening and follow-up of asymptomatic SDHx mutation carriers. Nat. Rev. Endocrinol. 17, 435–444 (2021).
Article CAS PubMed PubMed Central MATH Google Scholar
Davidoff, D. F. et al. Surveillance Improves Outcomes for Carriers of SDHB Pathogenic Variants: A Multicenter Study. J. Clin. Endocrinol. Metab. 107, e1907–e1916 (2022).
Article PubMed PubMed Central MATH Google Scholar
Selak, M. A. et al. Succinate links TCA cycle dysfunction to oncogenesis by inhibiting HIF-alpha prolyl hydroxylase. Cancer Cell 7, 77–85 (2005).
Article CAS PubMed Google Scholar
Xiao, M. et al. Inhibition of alpha-KG-dependent histone and DNA demethylases by fumarate and succinate that are accumulated in mutations of FH and SDH tumor suppressors. Genes Dev. 26, 1326–1338 (2012).
Article CAS PubMed PubMed Central MATH Google Scholar
Letouze, E. et al. SDH mutations establish a hypermethylator phenotype in paraganglioma. Cancer Cell 23, 739–752 (2013).
Article CAS PubMed Google Scholar
Liu, T. et al. The activating TERT promoter mutation C228T is recurrent in subsets of adrenal tumors. Endocr. Relat. Cancer 21, 427–434 (2014).
Article PubMed PubMed Central MATH Google Scholar
Calsina, B. et al. Genomic and immune landscape Of metastatic pheochromocytoma and paraganglioma. Nat. Commun. 14, 1122 (2023).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Fishbein, L. et al. Whole-exome sequencing identifies somatic ATRX mutations in pheochromocytomas and paragangliomas. Nat. Commun. 6, 6140 (2015).
Article ADS CAS PubMed MATH Google Scholar
Fishbein, L. et al. Comprehensive Molecular Characterization of Pheochromocytoma and Paraganglioma. Cancer Cell 31, 181–193 (2017).
Article CAS PubMed PubMed Central MATH Google Scholar
Flynn, A. et al. The genomic landscape of phaeochromocytoma. J. Pathol. 236, 78–89 (2015).
Article CAS PubMed MATH Google Scholar
Job, S. et al. Telomerase Activation and ATRX Mutations Are Independent Risk Factors for Metastatic Pheochromocytoma and Paraganglioma. Clin Cancer Res 25, 760–770 (2019).
Dwight, T. et al. TERT structural rearrangements in metastatic pheochromocytomas. Endocr. Relat. Cancer 25, 1–9 (2018).
Article CAS PubMed MATH Google Scholar
Dahia, P. L. et al. A HIF1alpha regulatory loop links hypoxia and mitochondrial signals in pheochromocytomas. PLoS Genet. 1, 72–80 (2005).
Article CAS PubMed MATH Google Scholar
Burnichon, N. et al. Integrative genomic analysis reveals somatic mutations in pheochromocytoma and paraganglioma. Hum. Mol. Genet. 20, 3974–3985 (2011).
Article CAS PubMed MATH Google Scholar
Castro-Vega, L. J. et al. Germline mutations in FH confer predisposition to malignant pheochromocytomas and paragangliomas. Hum. Mol. Genet. 23, 2440–2446 (2014).
Article CAS PubMed MATH Google Scholar
Castro-Vega, L. J. et al. Multi-omics analysis defines core genomic alterations in pheochromocytomas and paragangliomas. Nat. Commun. 6, 6044 (2015).
Article ADS CAS PubMed MATH Google Scholar
Comino-Mendez, I. et al. Exome sequencing identifies MAX mutations as a cause of hereditary pheochromocytoma. Nat. Genet 43, 663–667 (2011).
Article CAS PubMed Google Scholar
Zethoven, M. et al. Single-nuclei and bulk-tissue gene-expression analysis of pheochromocytoma and paraganglioma links disease subtypes with tumor microenvironment. Nat. Commun. 13, 6262 (2022).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Turchini, J., Cheung, V. K. Y., Tischler, A. S., De Krijger, R. R. & Gill, A. J. Pathology and genetics of phaeochromocytoma and paraganglioma. Histopathology 72, 97–105 (2018).
Article PubMed Google Scholar
Kameneva, P. et al. Single-cell transcriptomics of human embryos identifies multiple sympathoblast lineages with potential implications for neuroblastoma origin. Nat. Genet 53, 694–706 (2021).
Article CAS PubMed PubMed Central MATH Google Scholar
Eisenhofer, G. et al. Biochemical and clinical manifestations of dopamine-producing paragangliomas: utility of plasma methoxytyramine. J. Clin. Endocrinol. Metab. 90, 2068–2075 (2005).
Article CAS PubMed MATH Google Scholar
Richter, S. et al. Head/neck paragangliomas: focus on tumor location, mutational status and plasma methoxytyramine. Endocr. Relat. Cancer 29, 213–224 (2022).
Article CAS PubMed PubMed Central MATH Google Scholar
Fottner, C. et al. 6-18F-fluoro-L-dihydroxyphenylalanine positron emission tomography is superior to 123I-metaiodobenzyl-guanidine scintigraphy in the detection of extraadrenal and hereditary pheochromocytomas and paragangliomas: correlation with vesicular monoamine transporter expression. J. Clin. Endocrinol. Metab. 95, 2800–2810 (2010).
Article CAS PubMed Google Scholar
Taieb, D., Jha, A., Treglia, G. & Pacak, K. Molecular imaging and radionuclide therapy of pheochromocytoma and paraganglioma in the era of genomic characterization of disease subgroups. Endocr. Relat. Cancer 26, R627–R652 (2019).
Article CAS PubMed PubMed Central Google Scholar
Yokoyama, T., Misuzu, Y. Y. & Yamamoto, Y. Immunohistochemical localization of tryptophan hydroxylase and serotonin transporter in the carotid body of the rat. Histochem Cell Biol. 140, 147–155 (2013).
Article CAS PubMed MATH Google Scholar
van Hulsteijn, L. T., van Duinen, N., Romijn, J. A., Smit, J. W. & Corssmit, E. P. Urinary 5-HIAA excretion is not increased in patients with head and neck paragangliomas. Int J. Biol. Markers 27, e160–e163 (2012).
Article PubMed Google Scholar
Kimura, N. et al. Immunohistochemical Expression of Choline Acetyltransferase and Catecholamine-Synthesizing Enzymes in Head-and-Neck and Thoracoabdominal Paragangliomas and Pheochromocytomas. Endocr. Pathol. 32, 442–451 (2021).
Article CAS PubMed Google Scholar
Tamborero, D. et al. Cancer Genome Interpreter annotates the biological and clinical relevance of tumor alterations. Genome Med 10, 25 (2018).
Article PubMed PubMed Central MATH Google Scholar
Martincorena, I. et al. Universal Patterns of Selection in Cancer and Somatic Tissues. Cell 171, 1029–1041 e1021 (2017).
Article CAS PubMed PubMed Central MATH Google Scholar
Mularoni, L., Sabarinathan, R., Deu-Pons, J., Gonzalez-Perez, A. & Lopez-Bigas, N. OncodriveFML: a general framework to identify coding and non-coding regions with cancer driver mutations. Genome Biol. 17, 128 (2016).
Article PubMed PubMed Central Google Scholar
Yang, C. et al. Novel HIF2A mutations disrupt oxygen sensing, leading to polycythemia, paragangliomas, and somatostatinomas. Blood 121, 2563–2566 (2013).
Article CAS PubMed PubMed Central MATH Google Scholar
Welander, J. et al. Rare germline mutations identified by targeted next-generation sequencing of susceptibility genes in pheochromocytoma and paraganglioma. J. Clin. Endocrinol. Metab. 99, E1352–E1360 (2014).
Article CAS PubMed PubMed Central MATH Google Scholar
Pang, Y. et al. Nonmosaic somatic HIF2A mutations associated with late onset polycythemia-paraganglioma syndrome: Newly recognized subclass of polycythemia-paraganglioma syndrome. Cancer 125, 1258–1266 (2019).
Article CAS PubMed Google Scholar
Wilzen, A. et al. Malignant pheochromocytomas/paragangliomas harbor mutations in transport and cell adhesion genes. Int J. Cancer 138, 2201–2211 (2016).
Article CAS PubMed MATH Google Scholar
Bell, R. J. et al. Cancer. The transcription factor GABP selectively binds and activates the mutant TERT promoter in cancer. Science 348, 1036–1039 (2015).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
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).
Article CAS PubMed MATH Google Scholar
Sieverling, L. et al. Genomic footprints of activated telomere maintenance mechanisms in cancer. Nat. Commun. 11, 733 (2020).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Jang, S. H., Kim, A. R., Park, N. H., Park, J. W. & Han, I. S. DRG2 Regulates G2/M Progression via the Cyclin B1-Cdk1 Complex. Mol. Cells 39, 699–704 (2016).
Article CAS PubMed PubMed Central MATH Google Scholar
Calsina, B. et al. Integrative multi-omics analysis identifies a prognostic miRNA signature and a targetable miR-21-3p/TSC2/mTOR axis in metastatic pheochromocytoma/paraganglioma. Theranostics 9, 4946–4958 (2019).
Article CAS PubMed PubMed Central MATH Google Scholar
Chiang, C. H., Hou, M. F. & Hung, W. C. Up-regulation of miR-182 by beta-catenin in breast cancer increases tumorigenicity and invasiveness by targeting the matrix metalloproteinase inhibitor RECK. Biochim Biophys. Acta 1830, 3067–3076 (2013).
Article CAS PubMed Google Scholar
Guttilla, I. K. & White, B. A. Coordinate regulation of FOXO1 by miR-27a, miR-96, and miR-182 in breast cancer cells. J. Biol. Chem. 284, 23204–23216 (2009).
Article CAS PubMed PubMed Central Google Scholar
Qu, Y. et al. MiR-182 and miR-203 induce mesenchymal to epithelial transition and self-sufficiency of growth signals via repressing SNAI2 in prostate cells. Int J. Cancer 133, 544–555 (2013).
Article CAS PubMed MATH Google Scholar
Nolting, S. et al. Personalized Management of Pheochromocytoma and Paraganglioma. Endocr. Rev. 43, 199–239 (2022).
Article PubMed MATH Google Scholar
Kitange, G. J. et al. Induction of MGMT expression is associated with temozolomide resistance in glioblastoma xenografts. Neuro Oncol. 11, 281–291 (2009).
Article CAS PubMed PubMed Central Google Scholar
Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).
Article CAS PubMed PubMed Central MATH Google Scholar
Touat, M. et al. Mechanisms and therapeutic implications of hypermutation in gliomas. Nature 580, 517–523 (2020).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Backman, S. et al. The evolutionary history of metastatic pancreatic neuroendocrine tumours reveals a therapy driven route to high-grade transformation. J Pathol 264, 357–370 (2024).
Hanahan, D. & Weinberg, R. A. Hallmarks of cancer: the next generation. Cell 144, 646–674 (2011).
Article CAS PubMed MATH Google Scholar
Kim, N. W. et al. Specific association of human telomerase activity with immortal cells and cancer. Science 266, 2011–2015 (1994).
Article ADS CAS PubMed MATH Google Scholar
Dagg, R. A. et al. Extensive Proliferation of Human Cancer Cells with Ever-Shorter Telomeres. Cell Rep. 19, 2544–2556 (2017).
Article CAS PubMed MATH Google Scholar
Hiraga, S. et al. Telomerase activity and alterations in telomere length in human brain tumors. Cancer Res. 58, 2117–2125 (1998).
CAS PubMed MATH Google Scholar
Thompson, C. A. H. & Wong, J. M. Y. Non-canonical Functions of Telomerase Reverse Transcriptase: Emerging Roles and Biological Relevance. Curr. Top. Med Chem. 20, 498–507 (2020).
Article CAS PubMed MATH Google Scholar
Kimura, N. et al. Pathological grading for predicting metastasis in phaeochromocytoma and paraganglioma. Endocr. Relat. Cancer 21, 405–414 (2014).
Article PubMed MATH Google Scholar
Qin, T. et al. ATRX loss in glioma results in dysregulation of cell-cycle phase transition and ATM inhibitor radio-sensitization. Cell Rep. 38, 110216 (2022).
Article CAS PubMed PubMed Central MATH Google Scholar
Ohki, R. et al. Reprimo, a new candidate mediator of the p53-mediated cell cycle arrest at the G2 phase. J. Biol. Chem. 275, 22627–22630 (2000).
Article CAS PubMed MATH Google Scholar
Li, Y. et al. Non-canonical NF-kappaB signalling and ETS1/2 cooperatively drive C250T mutant TERT promoter activation. Nat. Cell Biol. 17, 1327–1338 (2015).
Article CAS PubMed PubMed Central MATH Google Scholar
Cao, Y., Li, H., Deb, S. & Liu, J. P. TERT regulates cell survival independent of telomerase enzymatic activity. Oncogene 21, 3130–3138 (2002).
Article CAS PubMed MATH Google Scholar
Khattar, E. et al. Telomerase reverse transcriptase promotes cancer cell proliferation by augmenting tRNA expression. J. Clin. Invest 126, 4045–4060 (2016).
Article PubMed PubMed Central MATH Google Scholar
Liang, J. et al. Genome-Wide CRISPR-Cas9 Screen Reveals Selective Vulnerability of ATRX-Mutant Cancers to WEE1 Inhibition. Cancer Res. 80, 510–523 (2020).
Article CAS PubMed MATH Google Scholar
Garbarino, J., Eckroate, J., Sundaram, R. K., Jensen, R. B. & Bindra, R. S. Loss of ATRX confers DNA repair defects and PARP inhibitor sensitivity. Transl. Oncol. 14, 101147 (2021).
Article CAS PubMed PubMed Central Google Scholar
George, S. L. et al. Therapeutic vulnerabilities in the DNA damage response for the treatment of ATRX mutant neuroblastoma. EBioMedicine 59, 102971 (2020).
Article CAS PubMed PubMed Central MATH Google Scholar
Cullinane, C. et al. Enhancing the anti-tumour activity of (177)Lu-DOTA-octreotate radionuclide therapy in somatostatin receptor-2 expressing tumour models by targeting PARP. Sci. Rep. 10, 10196 (2020).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Waldeck, K. et al. A genome-wide CRISPR/Cas9 screen identifies DNA-PK as a sensitiser to (177)Lutetium-DOTA-octreotate radionuclide therapy. Theranostics 13, 4745–4761 (2023).
Article CAS PubMed PubMed Central MATH Google Scholar
Chan, E. M. et al. WRN helicase is a synthetic lethal target in microsatellite unstable cancers. Nature 568, 551–556 (2019).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
Kamihara, J. et al. Belzutifan, a Potent HIF2alpha Inhibitor, in the Pacak-Zhuang Syndrome. N. Engl. J. Med. 385, 2059–2065 (2021).
Article CAS PubMed PubMed Central MATH Google Scholar
Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24–26 (2011).
Article CAS PubMed PubMed Central Google Scholar
Didion, J. P., Martin, M. & Collins, F. S. Atropos: specific, sensitive, and speedy trimming of sequencing reads. PeerJ 5, e3720 (2017).
Article PubMed PubMed Central MATH Google Scholar
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
Article CAS PubMed PubMed Central MATH Google Scholar
Lai, Z. et al. VarDict: a novel and versatile variant caller for next-generation sequencing in cancer research. Nucleic Acids Res. 44, e108 (2016).
Article PubMed PubMed Central MATH Google Scholar
DePristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet 43, 491–498 (2011).
Article CAS PubMed PubMed Central MATH Google Scholar
McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).
Article CAS PubMed PubMed Central MATH Google Scholar
Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat. methods 15, 591–594 (2018).
Article CAS PubMed MATH Google Scholar
Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219 (2013).
Article CAS PubMed PubMed Central MATH Google Scholar
Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
Article PubMed PubMed Central MATH Google Scholar
Cameron, D. L. et al. GRIDSS: sensitive and specific genomic rearrangement detection using positional de Bruijn graph assembly. Genome Res. 27, 2050–2060 (2017).
Article CAS PubMed PubMed Central MATH Google Scholar
Feuerbach, L. et al. TelomereHunter - in silico estimation of telomere content and composition from cancer genomes. BMC Bioinforma. 20, 272 (2019).
Article Google Scholar
Blokzijl, F., Janssen, R., van Boxtel, R. & Cuppen, E. MutationalPatterns: comprehensive genome-wide analysis of mutational processes. Genome Med. 10, 33 (2018).
Article PubMed PubMed Central Google Scholar
Alexandrov, L. B. et al. The repertoire of mutational signatures in human cancer. Nature 578, 94–101 (2020).
Article ADS CAS PubMed PubMed Central MATH Google Scholar
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).
Article PubMed PubMed Central MATH Google Scholar
Wintersinger, J. A. et al. Reconstructing Complex Cancer Evolutionary Histories from Multiple Bulk DNA Samples Using Pairtree. Blood Cancer Discov. 3, 208–219 (2022).
Article CAS PubMed PubMed Central Google Scholar
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Article CAS PubMed Google Scholar
Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
Article CAS PubMed MATH Google Scholar
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015).
Article PubMed PubMed Central MATH Google Scholar
Colaprico, A. et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 44, e71 (2016).
Article PubMed MATH Google Scholar
Kozomara, A., Birgaoanu, M. & Griffiths-Jones, S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 47, D155–D162 (2019).
Article CAS PubMed Google Scholar
Maksimovic, J., Gagnon-Bartsch, J. A., Speed, T. P. & Oshlack, A. Removing unwanted variation in a differential methylation analysis of Illumina HumanMethylation450 array data. Nucleic Acids Res 43, e106 (2015).
Article PubMed PubMed Central Google Scholar
Phipson, B., Maksimovic, J. & Oshlack, A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics 32, 286–288 (2016).
Article CAS PubMed MATH Google Scholar
Peters, T. J. et al. Calling differentially methylated regions from whole genome bisulphite sequencing with DMRcate. Nucleic Acids Res. 49, e109 (2021).
Article CAS PubMed PubMed Central MATH Google Scholar
Martelotto, L. G. Protocol for nuclei isolation from fresh and frozen tissues for parallel snRNA-Seq and snATAC-Seq on 10x ChromiumTM platform using the same nuclei preparation (updated version) V.5, https://doi.org/10.17504/protocols.io.bw52pg8e (2021).
Martelotto, L. G. Frankenstein’ protocol for nuclei isolation from fresh and frozen tissue for snRNAseq V.3, https://doi.org/10.17504/protocols.io.bqxymxpw (2020).
Stuart, T., Srivastava, A., Madad, S., Lareau, C. A. & Satija, R. Single-cell chromatin state analysis with Signac. Nat. methods 18, 1333–1341 (2021).
Article CAS PubMed PubMed Central Google Scholar
Stuart, T. et al. Comprehensive Integration of Single-Cell Data. Cell 177, 1888–1902 e1821 (2019).
Article CAS PubMed PubMed Central MATH Google Scholar
Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).
Article PubMed PubMed Central MATH Google Scholar
Kurtenbach, S. & Harbour, J. W. SparK: A Publication-quality NGS Visualization Tool. bioRxiv, 845529 (2019).
Fuchs, T. L. et al. Expanding the clinicopathological spectrum of succinate dehydrogenase-deficient renal cell carcinoma with a focus on variant morphologies: a study of 62 new tumors in 59 patients. Mod. Pathol. 35, 836–849 (2022).
Article CAS PubMed Google Scholar
Fuchs, T. L. et al. A Proposed Grading Scheme for Medullary Thyroid Carcinoma Based on Proliferative Activity (Ki-67 and Mitotic Count) and Coagulative Necrosis. Am. J. Surg. Pathol. 44, 1419–1428 (2020).
Article PubMed PubMed Central Google Scholar
Henson, J. D. et al. The C-Circle Assay for alternative-lengthening-of-telomeres activity. Methods 114, 74–84 (2017).
Article CAS PubMed MATH Google Scholar
Flynn, A. Multi-omic analysis of SDHB-deficient pheochromocytomas and paragangliomas identifies metastasis and treatment-related molecular profiles, github.com/UMCCR-RADIO-Lab/a5_sdhb_pcpg, https://doi.org/10.5281/zenodo.14890696 (2025).

Download references

Acknowledgements

The Clinical Genomics Platform at University of Melbourne Centre for Cancer Research for their support in generating primary data. Leah Meuter for her crucial work in curating clinical data for patients from the NIH. This work was supported by funds from the National Health and Medical Research Project Grant (APP1108032), Pheo Para Alliance (USA), Paradifference Foundation (Sweden) and The SDHB Coalition (USA). RWT was supported by a Victorian Cancer Agency Mid-Career Fellowship (TP828750). A.F. and A.Pa were supported by the Joseph Herman Trust at the University of Melbourne. K.P. was supported by the Intramural Research Program of the NICHD, NIH. P.L.M.D. is a Robert Tucker Hayes Distinguished Chair in Oncology and was supported by funds the NIH/NIGMS (GM114102), NIH/NCI (CA264248), Neuroendocrine Tumour Research Foundation, VHL Alliance Investigator Award, and Paradifference Foundation. L.F. was supported in part by ACS MRSG-15-063-01. T.D. and R.C.B. were supported by grants from Perpetual (Hillcrest) Foundation. The Novo Nordisk Foundation Center for Stem Cell Medicine, reNEW, is supported by a Novo Nordisk Foundation grant number NNF21CC0073729.

Author information

These authors contributed equally: Roderick Clifton-Bligh, Karel Pacak, Richard W. Tothill.

Authors and Affiliations

Centre for Cancer Research and Department of Clinical Pathology, University of Melbourne, VIC, Australia
Aidan Flynn, Andrew D. Pattison, Shiva Balachander, Emma Boehm, Blake Bowen, Fernando J. Rossello, Oliver Hofmann, Luciano Martelotto, Sean Grimmond & Richard W. Tothill
Kolling Institute of Medical Research, Royal North Shore Hospital St Leonards NSW, Melbourne, Australia
Trisha Dwight & Roderick Clifton-Bligh
Murdoch Children’s Research Institute, The Royal Children’s Hospital, Melbourne, VIC, 3052, Australia
Fernando J. Rossello
Novo Nordisk Foundation Centre for Stem Cell Medicine, Murdoch Children’s Research Institute, Melbourne, VIC, 3052, Australia
Fernando J. Rossello
Australian Regenerative Medicine Institute, Monash University, Victoria, Australia
Fernando J. Rossello
Peter MacCallum Cancer Centre, Melbourne, Australia
Maia Zethoven & Alison H. Trainer
Division of Endocrinology, Diabetes and Metabolism, Department of Internal Medicine, The Ohio State University, Columbus, OH, USA
Lawrence S. Kirschner
University of Michigan, Ann Arbor, MI, USA
Tobias Else & Thomas Giordano
Department of Medicine, Division of Endocrinology, Metabolism, Diabetes, University of Colorado, Aurora, CO, USA
Lauren Fishbein
Sydney Medical School, University of Sydney, Sydney, NSW, Australia
Anthony J. Gill & Roderick Clifton-Bligh
NSW Health Pathology, Department of Anatomical Pathology, Royal North Shore Hospital, St Leonards NSW, Sydney, Australia
Anthony J. Gill
Tufts Medical Center, Boston, MA, USA
Arthur S. Tischler
Eunice Kennedy Shriver National Institute of Child Health and Human Development, Bethesda, MD, USA
Tamara Prodanov & Karel Pacak
Children’s Medical Research Institute, Faculty of Medicine and Health, The University of Sydney, Westmead, Australia
Jane R. Noble & Roger R. Reddel
Sir Peter MacCallum Department of Oncology, University of Melbourne, Melbourne, VIC, Australia
Alison H. Trainer & Richard W. Tothill
University of Florida and Malcom Randall VA Medical Center, Gainesville, FL, USA
Hans Kumar Ghayee
Division of endocrinology and Research Center, Center hospitalier de l’Université de Montréal, Montreal, Canada
Isabelle Bourdeau
Waikato Clinical Campus, University of Auckland, Hamilton, New Zealand
Marianne Elston
Cancer Genetics Service, National Cancer Center Singapore, Singapore, Singapore
Diana Ishak & Joanne Ngeow Yuen Yie
Lee Kong Chian School of Medicine, Nanyang Technological University Singapore, Singapore, Singapore
Joanne Ngeow Yuen Yie
St Vincent’s Dept of Medicine, University of Melbourne, Melbourne, VIC, Australia
Rodney J. Hicks
Department of Medical Sciences, Uppsala University, Uppsala, Sweden
Joakim Crona
Department of Surgical Sciences, Uppsala University, Uppsala, Sweden
Tobias Åkerström & Peter Stålberg
Div. Hematology and Medical Oncology, Department of Medicine, Mays Cancer Center, University of Texas Health Science Center at San Antonio (UTHSCSA), San Antonio, TX, USA
Patricia Dahia

Authors

Aidan Flynn
View author publications
Search author on:PubMed Google Scholar
Andrew D. Pattison
View author publications
Search author on:PubMed Google Scholar
Shiva Balachander
View author publications
Search author on:PubMed Google Scholar
Emma Boehm
View author publications
Search author on:PubMed Google Scholar
Blake Bowen
View author publications
Search author on:PubMed Google Scholar
Trisha Dwight
View author publications
Search author on:PubMed Google Scholar
Fernando J. Rossello
View author publications
Search author on:PubMed Google Scholar
Oliver Hofmann
View author publications
Search author on:PubMed Google Scholar
Luciano Martelotto
View author publications
Search author on:PubMed Google Scholar
Maia Zethoven
View author publications
Search author on:PubMed Google Scholar
Lawrence S. Kirschner
View author publications
Search author on:PubMed Google Scholar
Tobias Else
View author publications
Search author on:PubMed Google Scholar
Lauren Fishbein
View author publications
Search author on:PubMed Google Scholar
Anthony J. Gill
View author publications
Search author on:PubMed Google Scholar
Arthur S. Tischler
View author publications
Search author on:PubMed Google Scholar
Thomas Giordano
View author publications
Search author on:PubMed Google Scholar
Tamara Prodanov
View author publications
Search author on:PubMed Google Scholar
Jane R. Noble
View author publications
Search author on:PubMed Google Scholar
Roger R. Reddel
View author publications
Search author on:PubMed Google Scholar
Alison H. Trainer
View author publications
Search author on:PubMed Google Scholar
Hans Kumar Ghayee
View author publications
Search author on:PubMed Google Scholar
Isabelle Bourdeau
View author publications
Search author on:PubMed Google Scholar
Marianne Elston
View author publications
Search author on:PubMed Google Scholar
Diana Ishak
View author publications
Search author on:PubMed Google Scholar
Joanne Ngeow Yuen Yie
View author publications
Search author on:PubMed Google Scholar
Rodney J. Hicks
View author publications
Search author on:PubMed Google Scholar
Joakim Crona
View author publications
Search author on:PubMed Google Scholar
Tobias Åkerström
View author publications
Search author on:PubMed Google Scholar
Peter Stålberg
View author publications
Search author on:PubMed Google Scholar
Patricia Dahia
View author publications
Search author on:PubMed Google Scholar
Sean Grimmond
View author publications
Search author on:PubMed Google Scholar
Roderick Clifton-Bligh
View author publications
Search author on:PubMed Google Scholar
Karel Pacak
View author publications
Search author on:PubMed Google Scholar
Richard W. Tothill
View author publications
Search author on:PubMed Google Scholar

Contributions

R.W.T., K.P., and R.C.B. conceived the study. A.F., A.D.P., S.B., E.B., B.B., T.D., F.R., L.M., M.Z. and J.R.N. performed experiments and data analysis. A.S.T. and A.J.G. performed pathology review and analysis. L.S.K., T.E., D.I., L.F., A.J.G., O.H., A.S.T., T.G., R.R.R., H.K.G., A.H.T., I.B., M.S.E., J.N.Y.Y., R.J.H., J.C., T.A., T.P., S.G., P.S., P.D., R.C.B., and K.P. provided material support, recruited patients, or curated clinical data. A.F., R.C.B., K.P., and R.W.T., wrote the manuscript. All authors edited and approved the final manuscript.

Corresponding authors

Correspondence to Roderick Clifton-Bligh, Karel Pacak or Richard W. Tothill.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature Communications thanks Eric Letouzé, Jong-Eun Park, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Description of Additional Supplementary Information

Supplementary Data 1

Supplementary Data 2

Supplementary Data 3

Supplementary Data 4

Supplementary Data 5

Supplementary Data 6

Supplementary Data 7

Reporting Summary

Transparent Peer Review file

Source data

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Flynn, A., Pattison, A.D., Balachander, S. et al. Multi-omic analysis of SDHB-deficient pheochromocytomas and paragangliomas identifies metastasis and treatment-related molecular profiles. Nat Commun 16, 2632 (2025). https://doi.org/10.1038/s41467-025-57595-y

Download citation

Received: 15 June 2024
Accepted: 26 February 2025
Published: 17 March 2025
Version of record: 17 March 2025
DOI: https://doi.org/10.1038/s41467-025-57595-y

Subjects

Abstract

Similar content being viewed by others

Introduction

Results

Whole genome and multi-omic analysis of SDHB-associated PCPG

Molecular profiling confirms genotype-subtype and cell-of-origin transcriptional profile

Genome-wide mutational features of SDHB-deficient PCPG

Somatic TERT mutations, structural rearrangement and promoter methylation

ATRX alterations are associated with the ALT phenotype

TERT/ATRX alterations are associated with metastatic progression and late somatic events

Transcriptional patterns in TERT/ATRX altered and metastatic PCPG

Treatment-related mutagenesis and resistance to DNA alkylating chemotherapy

Discussion

Methods

Patient samples and human research ethics

Tissue sample review and nucleic acid extraction

Targeted amplicon sequencing

Whole genome sequencing

Library preparation

WGS alignment and small variant calling

WGS copy number and structural variant calling

Telomere length estimation and telomere variant repeat detection

Mutational signature analysis

Microsatellite instability assessment

Identification of recurrent copy number alterations

Identification of driver genes

Clonal evolution analysis

RNA-seq

Library preparation

Read quantification

Differential gene expression

RNA fusion detection

Harmonization with published RNA-seq datasets

Small RNA-seq

Library preparation

Read quantification

Harmonization with published small-RNA-seq datasets

Differential expression of small-RNA dataset

Methylation analysis

Methylation array processing

Methylation array data pre-processing and normalization

Probe level differential methylation analysis

Differentially methylated regions

Harmonization with previously published methylation array datasets

Single nuclei RNA-sequencing

Library preparation with 10x chromium

Data processing and analysis

Single nuclei ATAC-sequencing

Library preparation with 10x Chromium

Data preprocessing

Quality control

Cell type identification

Peak calling and visualisation

Motif enrichment analysis from sn-RNA-seq data

Cell-Free DNA panel sequencing

Sample collection

Library preparation

Sequencing and data processing

Immunohistochemistry

C-circle assay

Statistics & reproducibility

Reporting summary

Data availability

Code availability

References

Acknowledgements

Author information

Authors and Affiliations

Contributions

Corresponding authors

Ethics declarations

Competing interests

Peer review

Peer review information

Additional information

Supplementary information

Source data

Rights and permissions

About this article