1. Sample cohort and clinical data
The CNIO series dataset comprises 114 samples from 105 patients with PPGL, previously reported by Calsina et al., 2023, supplemented with 81 samples obtained from 74 additional patients with PPGL available clinical information. For four patients, more than one tumor was available in the new sample set. Only tumors from patients who had provided informed consent were included in the study, in compliance with the principles of the Declaration of Helsinki. The protocol was approved by the following Ethics Committees: Hospital Universitario 12 de Octubre (15/024), Madrid, Spain; Universitäts Spital Zurich (2017-00771), Zurich, Switzerland; Klinikum der Universität (379-10), Munich, Germany; University of Würzburg (88/11), Würzburg, Germany; Azienda Ospedaliera Universitaria Careggi (Prot. N. 2011/0020149), Florence, Italy; Berlin Chamber of Physicians (Eth-S-R/14), Berlin, Germany; Radboud University Medical Center (9803-0060), Nijmegen, The Netherlands; University Hospital Carl Gustav Carus at TU Dresden (EK210052017 and EK189062010), Dresden, Germany. The annotated clinical data included tumor behavior, sex, age at first tumor diagnosis and presence of single or multiple PPGL from the 74 new patients. Tumors were classified as aggressive if they exhibited capsular or adipose tissue invasion, or vascular infiltration, or showed multiple recurrences without definitive evidence of metastatic disease (Table 1). In addition, we used clinical and omic data from publicly-available studies for subsequent analyses (Supplementary Figure 1).
Table 1. Clinical characteristics of 74 new patients incorporated into the CNIO series
Patient´s clinical characteristics
|
New CNIO patients (n=74)
|
Primary tumor localization
|
|
PCC
|
57% (42)
|
PGL
|
26% (19)
|
Multiple PGL
|
3% (2)
|
PCC+PGL
|
7% (5)
|
Metastasis
|
5% (4)
|
NA
|
3% (2)
|
Sex
|
|
Female
|
49% (36)
|
Male
|
47% (35)
|
NA
|
4% (3)
|
Age at initial diagnosis of PCC/PGL;
median (range) in years
|
48 (14-80)
|
Clinical behavior
|
|
Non-metastatic disease
|
68% (50)
|
Aggressive disease
|
7% (5)
|
Metastatic disease
|
26% (19)
|
2. Nucleic Acid Extraction and RNA-Seq Processing in CNIO series
The 81 new tumor samples used in the study were either formalin-fixed paraffin-embedded (FFPE) or fresh-frozen (FF). Tumor selection, ensuring the presence of more than 80% cancer cells, was performed by an expert pathologist (E.C.) based on hematoxylin and eosin-stained slides. DNA and RNA isolation and quality control were performed as previously described in Calsina et al., 2023(10).
RNA-Seq (n=81) library preparation was performed as previously reported by Calsina et al., 2023(10) using QuantSeq 3' mRNA-Seq Library Prep Kit FWD for Illumina (Lexogen, 015). The prepared libraries were subsequently sequenced on NovaSeq™ system (Illumina). To combine the new CNIO samples (n=81) with the CNIO+TCGA dataset generated in Calsina et al., 2023, RNA data were processed using ComBat-seq as described in that same article.
3. Molecular Characterization of PPGL series: Driver Pathogenic Variants and Secondary Events Involving ATRX and TERT
Germline and somatic PV in PPGL driver genes were characterized using NGS AmpliSeq Custom DNA Panel (Illumina, San Diego, CA, USA), and variants were confirmed by Sanger sequencing. Among the 81 new samples included, eight were negative for all tested drivers, hereafter referred to as wild-type (WT). TERT promoter- (C228T and C250T) and ATRX-PV were analyzed in all CNIO tumor samples using a customized NGS panel (12). While ATRX-PV data were available from TCGA and COMETE-derived series (Fishbein et al., 2017; Job et al., 2019), TERT promoter mutational status was only described in the latter series (14). However, to acquire complete information on TERT status, TERT expression, TERT promoter hypermethylation and/or copy number gains were also analyzed in the whole series (12–15).
4. MAML3-fusion Identification
Identification of tumors harboring a MAML3-fusion was performed by a combination of different approaches: multi-omic (transcriptome (n=765) and methylome (n=568)), break-apart fluorescence in situ hybridization (FISH), and PD-L1 immunohistochemistry (IHC). Sample workflow is summarized in Supplementary Figure 1. Globally, the combined series included data of 779 unique patients.
Transcriptomic analysis
Transcriptomic analysis to identify fusion cases in the whole CNIO series was conducted by combining the RNA-Seq from Calsina et al., 2023 and 81 non-published CNIO samples, together with data from several publicly available data sets: RNA-Seq from the TCGA series (3), microarray data from E-MTAB-733 (16), GSE 67066 (13) and CNIO microarray series (GSE19422) (17) (Supplementary Figure 1). These six series were combined by Z-score, transforming the expression values for each gene for each cohort (centered at the mean expression), selecting the genes described in the TCGA transcriptomic profile (DVL3, MYC, WNT4, NKD1, FZD7, WNT5B, WNT7B, MAML3, PTCH1, GLLI2, FZD3, FZD8 and miR-375) as associated with MAML3-fusion (3). MicroRNA (miRNA) data was extracted from Calsina et al., 2019.
DNA methylation analysis
Methylation analysis was carried out to identify MAML3-fusion cases based on a hypo-methylated probe signature derived from Fishbein et al. 2017 data set. We identified a total of ≈3000 differentially methylated probes included in the Illumina Infinium HumanMethylation450 BeadChip among the MAML3-tumors versus other genotypes (log Fold Change (|LogFC|)>1 and False Discovery Rate (FDR) < 0.05). Of these ≈3000 probes, 243 were included in the Illumina Infinium HumanMethylation27 BeadChip. The beta-value of each of the 243 probes per sample was Z-score transformed (centered at the mean beta-value of each probe for each series). The CNIO (19–22), TCGA(3) and GSE43298 (9) series were pooled in a single matrix for sample clustering.
Break-apart FISH assay for MAML3-fusion
To detect MAML3-fusions a dual-color break-apart FISH probe was developed using six bacterial artificial chromosome (BAC) clones from BACPAC Genomics (Supplementary Figure 2a). Analysis was performed in FFPE-whole slides or tissue microarrays (TMA) cores from available samples.
Three BAC clones (RP11-6A22, RP11-615B12, and RP11-903H21) located in MAML3 5' region were labeled with Spectrum-Green using the Nick translation assay (Sigma Aldrich). The remaining three clones (RP11-625H13, RP11-876B4, and RP11-28L10) located in the 3' region were labeled with Spectrum-Orange.
FFPE sections were mounted on positively charged slides (SuperFrost, Thermo Scientific) following the Histology FISH Accessory Kit (DAKO) instructions as previously described (23). Tissue sections were deparaffinized in xylene, rehydrated through a series of ethanol solutions, and treated with 2-[N-morpholino] ethanesulphonic acid (MES) and pepsin for protein digestion. Then, samples were denatured in the presence of the FISH probe and allowed to hybridize overnight in a humid chamber. Subsequently, the slides were washed and mounted with DAPI (Sigma Aldrich).
Two independent investigators (MCM and SR-P) manually scored the FISH signals by examining the number of nuclei with split signals, assessing presence of MAML3-fusions.
PD-L1 Quantification by Immunohistochemistry (IHC)
PD-L1 IHC positivity evaluation was performed according to the criteria described by Calsina et al., 2023. Each case was scored semi-quantitatively depending on the tumor proportion score (TPS), as tier 1 (0–15% positive tumor cells), tier 2 (15–30%), and tier 3 (30–100%). Analysis was performed in FFPE-whole slides or TMA cores from available samples.
5. MAML3-fusion Validation Techniques
PCR
Total RNA isolated from tumors (previously described) was reverse-transcribed and used as template for a subsequent PCR with specific primers: upstream UBTF portion (exon 17: 5′-GGAGCAGCAAAAGCAGTACA-3′ and exon 19: 5′-CCCAAACCCCCCAAATCCAG-3′) and downstream MAML3 portion of the fusion transcript (5′ TCTCCATTAAGTGGTGGTGGATCGAG 3′). Amplicons were Sanger sequenced to validate the fusion.
RNA Fusion NGS Panel for MAML3-fusion Validation
A customized Archer's FusionPlex panel (FusionPlex MAML3 Supplement 22729 concatenated with FusionPlex Lung v2 18090 v1.3 v1.0, Archer®) was used to validate the presence of MAML3-fusions in samples. The design incorporated a spike-in panel consisting of primers that covered the five MAML3 gene exons, both in 5' and 3' directions.
Library preparation was carried out using 300 ng of RNA input per sample and the Archer Universal RNA reagent Kit v2 and Archer Molecular Barcode (MBC) Adapters for Illumina, following the manufacturer’s instructions. Libraries were quantified using the KAPA Library Quantification Kit from Roche, and equimolar quantities of the libraries were pooled together.
Sequencing was performed on the MiSeq™ System from Illumina, ensuring a minimum depth of 1 million paired-end reads with a read length of 150 bp. The generated sequencing data were subsequently analyzed using Archer Analysis v7 software developed by ArcherDX. Samples were considered as positive when the fusion breakpoint was supported by a minimum of three unique start site reads.
6. Additional OMIC-analysis of the CNIO+TCGA Transcription Matrix
RNA-Seq Analyses and TME Profiling
The platforms used for the previous transcriptomic study are composed of microarrays (E-MTAB-733, GSE 67066 and GSE19422) and RNA-Seq (CNIO and TCGA) data. Due to the impossibility of unifying data from both technologies for whole-transcriptomic analyses, microarray matrices were excluded for further analysis. Thus, we used the CNIO+TCGA RNA-Seq data, previously described by Calsina et al., 2023 and we included 81 new samples, resulting in a total of 197 tumors in the CNIO series.
CNIO+TCGA matrix was composed by 349 tumors, hereafter referred to as “Series 2”. This series included a balanced number of non-metastatic and metastatic patient samples, used for subsequent differential expression analysis (DE) using DESeq2 (24), gene set enrichment analysis (GSEA) and immune cell-type estimation. Study workflow is summarized in Supplementary Figure 1.
Gene set enrichment analysis
GSEA was performed in Series 2 using software version 2.2.2. (RRID:SCR_003199), focusing on the Molecular Signature Database (MSigDB) collections, specifically the 'H: hallmark gene sets'. Analysis was conducted using default settings, with the exception of the collapse set, which was set to false. The scoring_scheme was set to classic, plot_top_x was set to 100, set_max was set to 1000, and set_min was set to 10. The comparisons were made between MAML3-tumors and other genotypes using the RNA-Seq matrix data.
Immune Cell Type Estimation
Absolute abundance of 22 different immune cell subtypes was calculated using CIBERSORTx (25) in each sample from Series 2, using the DESeq2 normalized expression matrix. CIBERSORTx analysis parameters were the following: (1) Gene signature matrix: LM22, containing gene expression profiles of 22 immune cell subsets; (2) 1000 permutations were performed to obtain robust and reliable results; (3) B-mode batch correction to account for any potential batch effects and (4) Quantile normalization was disabled during the analysis.
7. Immunohistochemistry and image data quantification
IHC study was performed on a representative FFPE-sample subset with available material. Sections measuring 2-3 microns thick were prepared and dried overnight at 60°C. IHC staining was performed using the BOND-MAX Automated Immunohistochemistry Vision Biosystem from Leica Microsystems GmbH and AS-Link 48 from DAKO, following the standard protocols.
Sections were initially deparaffinized and pretreated with Epitope Retrieval Solution 2 (EDTA buffer pH 8.8) at 98°C for 20 minutes. After washing, peroxidase blocking was conducted for 10 minutes using the Bond Polymer Refine Detection Kit DC9800 from Leica Microsystems GmbH.
Tissues were washed and then incubated with the appropriate primary antibodies: mouse anti-CD31 (Ready to Use (RTU), mouse monoclonal [JC70A], Dako (IR610/IS610)), rabbit anti-PD-L1 (dilution: 1:150, rabbit monoclonal [E1L3N], Cell Signaling Technology (13684)), mouse anti-CD8 FLEX (RTU, mouse monoclonal [C8/144B], Dako (IR623), visualization system: EnVision FLEX + Mouse + Magenta), mouse anti-CD68 (RTU, mouse monoclonal [KP1]), mouse anti-CD163 (dilution: 1:200, mouse monoclonal [10D6], Leica (NCL-CD163)), mouse anti-perforin (dilution: 1:50, Mouse Monoclonal [5B10], Abcam (ab89821). In order to discriminate between CD8 lymphocytes and NK/Cytotoxic cells, which also express CD8, a double IHC using CD8/perforin was performed. Positive cells for CD8/perforin or perforin alone were considered NK/Cytotoxic cells.
Tissues were incubated with the anti-mouse or anti-rabbit secondary antibody with polymer or envision flex for 10 minutes and developed with DAB or magenta chromogen for 10 minutes. The slides were then stained with hematoxylin, dehydrated, cleared, and mounted with a permanent mounting medium for microscopic evaluation. Whole slide images were scanned using the AxioScan Z1 from Zeiss and captured with either Zen Blue software (ZEN Digital Imaging for Light Microscopy, Zeiss, RRID:SCR_013672) or AxioVision software (AxioVision Imaging System, RRID:SCR_002677).
Vascular patterns (CD31) in whole tissue sections were examined by an expert pathologist (E. C.). Blood vessels were analyzed using a workflow combining QuPath (26) and Fiji (Image J). Five regions of interest (ROIs) were selected for each tumor, and vessel detection was performed on each ROI using QuPath. The resulting images from each ROI were then analyzed in Fiji through skeletonization process with default settings. Skeletonization data were further processed using GraphPad Prism 9 and Microsoft Excel obtaining vessel length, vessel number and vessel branching data, subsequently combined to create the classification of the seven vascular phenotypes described in this article.
CD8+, CD68+, and CD163+ cell quantification was performed using a digital image analysis workflow developed in QuPath. This analysis was conducted in the same tumor region of each slide (ROI of 2.5x107 µm2) stained with the corresponding antibody. The pipeline automatically detected each cell, stroma, and vessels, and quantified positive cells for each mentioned antibody among the detected cells. Only tumor-infiltrating cells were considered, excluding regions with immune nodules and immune cells detected within or in the immediate periphery of blood vessels.
8. Quantification, Statistical Analysis and Sample Clustering
Statistical analyses were conducted using different software packages: GraphPad Prism 9, R version 3.2.2., and IBM SPSS Statistics v19.
Categories for p-values are the following: p-value <0.001: denoted as *** indicating highly significant results, p-value between 0.001 and 0.01, denoted as ** showing very significant results and p-value between 0.01 and 0.05, denoted as * indicating significant results. Type of test performed and p-values are described for each experiment.
Time to progression curves (based on Series 1) were generated, according to the time interval between primary tumor diagnosis and the appearance of a metastatic lesion and/or recurrence. The risk of developing early metastases or aggressive disease was assessed by univariate logistic regression model.
Morpheus software, available at https://software.broadinstitute.org/morpheus, was used for graphical visualization, classification, and sample clustering.