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

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.

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

A Outer ring: Overview of the analytical methods applied to the cohort annotated with the number of tumours analysed. Inner ring: The number of tumours and patients (in parenthesis) analysed from each anatomical location. Where a metastasis was analysed, the location of the primary tumour is indicated. Created with components from BioRender. (Licence: Tothill, R. (2025), https://BioRender.com/y24d952). B A summary of the combination of analytical methods applied to each tumour with respect to anatomical primary location. The upper panel indicates the total number of samples from each anatomical location analysed with the respective combination of assays (lower panel). C Anatomical location (dot colour) and clinical behaviour (tile colour) of each tumour analysed. Paired samples from the same patient are joined by a line. D Summary of the germline SDHB mutations detected by WGS across the cohort and their respective position within the protein amino acid sequence. (i) Amino acid changes for single nucleotide variants and small insertion/deletion events (NP_002991.2). (ii) The total number of patients observed with each amino acid change where bar colour indicates the clinical disease course. (iii) Schematic of the consequence of splice donor mutations. (iv) Regions of the SDHB protein sequence deleted by large scale structural events. Bar colour indicates the clinical disease course and the number of patients affected is indicated within the bar.

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

UMAP dimensional reduction was used to cluster WTS (A, n = 315), small-RNA-seq (B, n = 416) and DNA methylation data (C, n = 413) with previously published data^,,,,,. UMAP clustering was also repeated for respective A5 data types in isolation (D–F, n = {D: 91, E: 90, F: 93}). G Differential expression profiling was performed between abdominal-thoracic PCPG and HN-PG. The heatmap shows CPM (log₂, z-score) expression values for each tumour (x-axis) for the top differentially expressed genes (limma moderated t-test Benjamini-Hochberg adjusted p-value < 0.05, ranked by log-fold-change, top and bottom 30 genes are shown)(y-axis). The top annotation bars indicate the suspected cell-of-origin based on UMAP clustering and the anatomical location of the tumour, respectively. The left annotation bar indicates whether the gene was also in a differentially methylated region for the same contrast. H, I Spatial distribution of adjusted p-values (y-axis, -log10, RNA-seq p-values from limma Benjamini-Hochberg corrected moderated t-test, methylation p-values derived using the RUV-inverse method from the missMethyl package) along chromosomes 7, 12, and 17 (x-axis) from differential expression (H) and probe-level differential methylation analysis (I) between sympathetic PCPG (n = 75) and HN-PG (n = 12). J Expression of catecholamine biosynthesis and processing pathway genes. Line colour indicates the anatomical location of the tumour and sub-panels segregate tumours based on which catecholamines were above upper normal limit during clinical testing. Expression data for the A5 cohort were combined with a larger compendium of publicly available data representing the different PCPG subtypes (n = 315) before values were normalized to Z-scores. Only A5 tumours are shown.

Fig. 3. Genomic and clinical features of PCPG.

A (i) Location of the primary tumour. An asterisk indicates whether the tumour had a non-chromaffin expression signature. (ii) Clinical behaviour. Primary tumours from patients with metastatic disease are annotated as to whether the tumour was able (filled circle) or unable (open circle = no metastasis sequenced, strike-through = confirmed unrelated) to be confirmed as the metastatic clone through sequencing of a paired metastasis. (iii) Identifier linking samples from the same patient. (iv) Indicates if the sample was resected after cytotoxic treatment. (v) Immunohistochemical scoring of Ki67; greater than 3% positive cells is indicated by a black square. (vi) Tumour mutation burden (mutations per megabase) and (vii) number of structural variants are indicated using a heatmap. (viii–x) Signature analysis was performed using the MutationalSignatures package. Signatures that contributed 15% of mutations and at least 500, 50, or 10 mutations for SBS, ID, and DBS signatures, respectively, in at least one tumour are shown. Signatures are ranked (top to bottom) based on the mean proportional contribution. Colour indicates the proportion of mutations attributable to each signature (xi) TERT/ATRX mutation status. B Aggregate and per-sample copy number aberrations in sympathetic (rows 2 and 3, respectively) and parasympathetic (rows 5 and 6, respectively) PCPG. Statistically significant recurrently deleted (downward arrow) or amplified (upwards arrow) regions are indicated (rows 1 and 4). Clinical behaviour and whole genome doubling are annotated to the left of per-sample data (rows 3 and 5). The percentage of genome altered and structural variants count are annotated on the right. Losses and gains as well as percent genome altered are relative to sample ploidy rather than diploid. C Genes that were recurrently altered and predicted to be drivers by Cancer Genome Interpreter across the A5 cohort (top) and the respective number of mutations observed in published datasets (bottom). D Schematic showing somatic mutations observed in EPAS1.

Fig. 4. TERT promoter mutations and structural alterations.

A sn-ATAC-seq covering the TERT gene and promoter (inset) region. Data is shown for tumour cells from two ATRX mutant tumours, three TERT mutant tumours, and one tumour with a TERT structural variant. Normal chromaffin, adrenocortical, fibroblast, and endothelial cells are shown for contrast. B (i) Circos plot describing somatic alterations detected by WGS in E123-M1. Outer ring indicates the chromosome, second outer ring marks SNVs (blue: C > A, black: C > G, red: C > T, grey: T > A, green: T > C, pink: T > G) and Indels (yellow: insertions, red: deletions). Third outer ring indicates total copy number, (green: >2, red: <2), the fourth outer ring indicates the copy number of the minor allele (blue: >1, orange: <1). The centre circle displays structural variants. (ii) Copy number segmentation data from the AFF4 and TERT gene regions in E123-M1 showing multiple short segmental copy-number changes due to chromothripsis. The break points contributing to an AFF4-TERT fusion are marked with red dashed lines. (iii) Schematic of AFF4-TERT fusion detected by WTS. C TERT expression (y-axis, log₂ CPM) versus gene copy number at the TERT locus (x-axis). D TERT expression (y-axis, log₂ CPM) versus the methylation status (x-axis) of probes in the TERT promoter and TERT hypermethylated oncological regions (rows 1 and 3). Density plot of probe methylation beta-values (rows 2 and 4). Point and line colour indicate TERT/ATRX mutation status.

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

A Genomic (top) and protein (bottom) position of alterations in ATRX. B Tumour/normal ratio of telomere read content (y-axis) in relation to ATRX^mut (n = 10), TERT^mut (n = 20) and wild-type (n = {Chromaffin: 48, Non-chromaffin:14}) status (x-axis). Presence (red) or absence (black) of C-circles is indicated by dot colour. Tumours from the same patient are joined by a line. Values were averaged across paired tumours and a two-sided Student’s t-test was used to test for differences. C Telomere variant repeats (TVRs) of the type NNNGGG were detected in WGS from tumour and matched normal using TelomereHunter. Count values were normalized against intratelomeric (green) or total (blue) read count. The tumour/normal ratio of normalized counts (y-axis) are shown with respect to the presence (n = 12) or absence (n = 76) of detected c-circles (x-axis). Only TVRs which were significantly different using a two-tailed Students t-test (p < 0.05) applied to values normalized to total reads are shown. TVRs that were significant after false-discovery correction (Benjamini & Hochberg, p < 0.1) are indicated with a double-asterisk. P-values before (p) and after (pBH) correction are shown. Data points are coloured to indicate a tumour/normal telomere content ratio (log₂) greater than (red) or less than (grey) 0.5. D Reads containing telomeric sequences were counted from WTS data using TelomereHunter. The telomeric content (y-axis) was computed as the number of unmapped reads containing telomeric sequence times 1,000,000 divided by the total number of reads with a GC content similar to telomeric repeats (n = {TERT^mut: 20, ATRX^mut: 10, Chromaffin^wildtype: 45, Non-chromaffin^wildtype: 14}). Data point colours, x-axis, and statistical testing are as described in (B). The hinges of each boxplot correspond to the first and third quartiles and the median value is marked. The whiskers extend to the largest and smallest value no greater than 1.5 times the interquartile range above or below the upper and lower hinges, respectively.

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

A A schematic of data availability and TERT/ATRX status in primary tumours reported as metastatic. B A metastatic C-circle positive and TERT-mutant polyclonal primary: (i) The lower-left panel shows the variant allele frequency (VAF, y-axis) for all somatic mutations in the metastatic primary (E143-P1) and paired metastases (E143-M1/2) (x-axis). Copy number status in the region of each variant is indicated by dot colour. Horizontal lines indicate tumour purity as a proportion (read from y-axis). Grey lines connect shared mutations while a pTERT mutation is highlighted with a red line. The top panel shows the tumour/normal telomere content ratio, the colour of each bar indicates C-circle status. The panel on the right shows the genomic copy number status along each chromosome (y-axis) for each tumour (x-axis). (ii) A schematic illustration of the clonal evolution of metastatic disease in patient E143. Cell colour indicates the presence of the ALT phenotype (green) or pTERT mutation (red). C The tumour mutation burden (mutations per megabase, y-axis) observed in each tumour with respect to ATRX^mut (n = 10), TERT^mut (n = 20) and wild type (n = {Chromaffin: 48, Non-chromaffin: 14}) status (x-axis). A one-tailed Student’s t-test was used to test for statistical significance. The y-axis has been truncated to accommodate an extreme outlier which was excluded during statistical testing. D The dimensions (centimetres, y-axis) of the largest primary tumour reported for each patient. Patients are stratified by the presence of a TERT/ATRX mutation and the presence (M) or absence (NM) of metastatic disease (n = {Chromaffin TERT(M): 13, ATRX(M): 5, WT(M): 9, WT(NM): 26; Non-chromaffin: WT(M): 4, WT(NM): 9}). The hinges of each boxplot correspond to the first and third quartiles and the median value is marked. The whiskers extend to the largest and smallest value no greater than 1.5 times the interquartile range above or below the upper and lower hinges, respectively. A one-tailed Student’s t-test was used to test for statistical significance.

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

Differential gene expression contrasting non-metastatic primary tumours with either (i) ATRX mutant tumours, (ii) TERT mutant tumours, or (iii) all metastases and metastatic primaries. A fourth contrast between TERT and ATRX mutant tumours was also performed. A An upset plot showing the intersection of genes that were significant (limma moderated t-test Benjamini-Hochberg adjusted p-value < 0.05, log-fold-change > 1) in each contrast. Bar colour indicates gene-ontology association for protein coding genes or gene biotype for non-protein-coding genes. B Differential gene-expression between non-metastatic primary and metastatic tumours showing fold change (log₂, y-axis) versus p-value (-log₁₀, x-axis, limma moderated t-test Benjamini-Hochberg adjusted). Genes that were also significant in non-metastatic primary vs TERT altered and non-metastatic primary versus ATRX altered contrasts are coloured according to gene ontology annotation. C Heatmap (centre panel) showing genes differentially expressed in both the ATRX^mut vs non-metastatic primary and ATRX^mut vs TERT^mut contrasts (n = {ATRX^mut: 9, TERT^mut: 15, non-metastatic primary: 21}). Annotation bars (left) indicate whether the gene was found in a differentially methylated region (ATRX^mut vs non-metastatic primary), and the correlation of the expression of each gene to that of MKI67 (Spearman correlation). The right panel shows snRNA-seq expression aggregated by cell type (right sub-panel) or ATRX/TERT mutation status (left sub-panel). Dot colour indicates mean expression while dot size indicates the fraction of cells expressing the gene (nCells = {ATRX^mut: 29533, TERT^mut: 19177, Tumour-wild-type: 14916; Tumour: 63626, Myeloid cells: 2198, Chromaffin cells: 3238, Adrenocortical cells: 3723; Schwann-cell-like cells (SCLCs): 727, Fibroblasts: 1855, Endothelial cells: 1849, Lymphocytes: 1200}. D Expression (log₂ CPM, y-axis) of genes differentially expressed between ATRX altered and non-metastatic primary tumours. Expression data is shown from the A5 cohort (top) and the TCGA/Flynn et al. cohorts (bottom). Data point colour indicates PCPG subtype. A two-sided Student’s t-test was used to test for differences (n = {A5 - WT: 81, ATRX^mut: 10, TCGA/Flynn - WT: 217, ATRX^mut: 6}) E Differentially expressed genes in both the TERT-altered vs non-metastatic primary and ATRX^mut vs TERT^mut contrasts (n = {ATRX^mut: 9, TERT^mut: 15, non-metastatic primary: 21}). See panel (C) description for panel elements. Box plots to be interpreted as per Fig. 6.

Fig. 8. Evolution of PCPG under treatment pressure.

A Clinical timeline of for patient E169. Bi Shared and private coding, non-coding/synonymous (NC/Syn.), and structural (SV) variants between paired metastases taken before and after CVD treatment. Bii Copy number status along each chromosome (y-axis) in paired metastases (x-axis). C Mutation signature analysis using COSMIC v3 SBS (i) and InDel (ii) signatures (y-axis) in paired metastases (x-axis). Heatmap colour indicates signature contribution. D Expression of MGMT in the A5 cohort. E–G Patient E167, see description for panels (B, C). H Expression of MLH1 in the A5 cohort. I Total mutation counts (y-axis) for PCPG tumours (x-axis) in the Project GENIE data registry. Bar colour indicates the presence of a mutation in the mismatch repair pathway. J Trinucleotide context for mutations observed in E167-M1 (top), E167-M2 (second from top), the highest mutation load tumour from the GENIE dataset (second from bottom), and ctDNA derived from a patient with metastatic SHDB-related PGL treated with Temozolomide. K¹⁸F-FDG-PET imaging for Temozolomide treated patient at time of blood draw for cfDNA analysis. L IchorCNA analysis of ctDNA derived from Temozolomide treated patient (M) Variant allele frequencies (y-axis) for somatic variants (x-axis) observed in ctDNA derived from Temozolomide treated patient. Datapoints are coloured to indicate transition/transversion or insertion/deletion type.