Enlibio,Your biological research partner!
ELISA
Enlibio Biotech Co.,Ltd,Wuhan,China 2008-2023
R&D 15th year
HomeHome
ProductsNav
NewsNav
Contact UsNav
About UsNav
DistributorsNav
SitemapNav
Human ELISA KitsMouse ELISA kitsRat ELISA kitsOther species ELISA kitsELISA Kit Assist ReagentsELISAsPrimary AntibodySecondary AntibodyIndustrial Antibody developmentAntibody Development ServicesAntibody Related ReagentsWestern blot ReagentsImmunohistochemistry ReagentsAntibodiesRecombinant ProteinsNatural ProteinsPeptidesProtein servicePeptide serviceImmunogensCell LineDownload ManualEnlibio NewsBiological abstractContact UsTechnical SupportMessagesAbout UsJob OpportunityCompany structurePicture show

Genomic and immune landscape Of metastatic pheochromocytoma and paraganglioma

Issuing time:2023-03-03 10:44

Abstract

The mechanisms triggering metastasis in pheochromocytoma/paraganglioma are unknown, hindering therapeutic options for patients with metastatic tumors (mPPGL). Herein we show by genomic profiling of a large cohort of mPPGLs that high mutational load, microsatellite instability and somatic copy-number alteration burden are associated withATRX/TERTalterations and are suitable prognostic markers. Transcriptomic analysis defines the signaling networks involved in the acquisition of metastatic competence and establishes a gene signature related to mPPGLs, highlighting CDK1 as an additional mPPGL marker. Immunogenomics accompanied by immunohistochemistry identifies a heterogeneous ecosystem at the tumor microenvironment level, linked to the genomic subtype and tumor behavior. Specifically, we define a general immunosuppressive microenvironment in mPPGLs, the exception being PD-L1 expressingMAML3-related tumors. Our study reveals canonical markers for risk of metastasis, and suggests the usefulness of including immune parameters in clinical management for PPGL prognostication and identification of patients who might benefit from immunotherapy.

Introduction

Pheochromocytomas and paragangliomas (PPGLs) are rare and highly heterogeneous neuroendocrine tumors, associated with mutations in one of the at least 20 driver genes related to the disease1. Up to 20% of PPGLs are metastatic (mPPGLs), with a heterogeneous 5-year overall survival rate that ranges from 40 to 77%2. Although some genotype-phenotype correlations have been defined for patients with PPGLs, disease management remains challenging due to a variable clinical behavior, lack of accurate markers of metastasis risk and inefficient therapeutic standards for metastatic stage3,4. Moreover, a poor understanding of the underlying mechanisms driving metastasis hinders the development of an efficient therapeutic strategy for mPPGLs.

The genomic landscape of PPGLs has been poorly described and little is known about the role that secondary events play in the progression of the disease. So far, only two large-scale multi-omic studies reporting the genomic characteristics of PPGLs have been published5,6. These studies established the molecular classification of PPGLs in three genomic subtypes: pseudohypoxic, kinase signaling, and Wnt-altered. However, these previous studies included few metastatic cases. Moreover, the tumor microenvironment (TME) has been shown to play a major role in the prognosis and therapeutic stratification of other tumor types7. However, the TME remains largely unexplored in PPGLs and may be pivotal for therapy success.

In this large-scale study, we performed a comprehensive characterization of the genomic landscape of mPPGLs and interrogated the immune compartment of these tumors. We compared the results from our series, enriched in metastatic cases, with the TCGA PPGL cohort5, which is mainly composed of non-metastatic tumors. Using this strategy, we defined the mutational architecture of mPPGL, identified transcriptional alterations linked to metastasis-related processes, proposed prognostic markers, and outlined the immune infiltration profile in PPGL.

Results

We performed genomic profiling of 156 PPGLs from 128 unrelated patients: whole-exome sequencing (WES) in 87 paired germline–tumor samples and RNA sequencing (RNA-Seq) in 114 tumor samples (hereinafter, CNIO cohort; Fig.1a; Supplementary Table1). Seventy-five of these patients had metastatic disease, detected at initial diagnosis or during follow-up, and accounted for 99 metastatic samples included in the study (either primary tumors, relapses or metastases). We also included in the study data of 175 PPGLs from 174 unrelated patients from the PPGL TCGA project: WES data from 174 paired germline–tumor samples and RNA-Seq data from 150 tumors (from now on, TCGA cohort). By a comprehensive analysis of these data, we deciphered the mutational, copy number alteration (CNA), transcriptomic and immune landscapes of mPPGLs (Fig.1b). Tumors were classified according to their genotype/ transcriptional profile into the previously defined PPGL genomic subtypes5.

aCharacteristics of the CNIO cohort used for WES and RNA-Seq analysis. Characteristics shown are: platform available (WES and RNA-Seq), preservation of material (FFPE or frozen), tumor type (non-metastatic, aggressive, metastatic, relapse or metastases), primary tumor and metastasis location, sex, age at diagnosis, genomic subtype and genotype. NA: information not available.bGeneral overview of the study. See also Supplementary Table1and Supplementary Fig.1.

Mutational landscape of metastatic pheochromocytoma and paraganglioma

We identified somatic single-nucleotide variations (SNVs) and small insertions/deletions (INDELs) in 87 paired germline-tumor samples with available WES data (mean depth of 105×; 12 tumors from 11 patients were discarded due to poor quality; see methods section). Data from 169 TCGA patients (174 samples) were added for the analyses (Supplementary Fig.1). To note, for 20 patients (16 from CNIO cohort and 4 from TCGA cohort), we had available more than one tumor sequenced by WES (range, 2–7).

WES data analysis revealed 5618 somatic variants across 261 samples. From those, 3641 had a variant allele frequency (VAF) > 10% and 2227 were missense, 140 truncating, 127 affected start/stop codons, 52 in-frame deletions or insertions, 66 altered canonical splice sites, and 1029 were classified as ‘other’ (including synonymous variants); 3.2% of the variants were classified as high impact consequence variants and affected cancer genes (from COSMIC Cancer Gene Census or OncoKB).

WES data was used to calculate the tumor mutational burden (TMB) and microsatellite instability (MSI) scores of each tumor (Fig.2a). The TMB (based on true-somatic coding events with VAF > 0.1; category 5 in Supplementary Fig. 2) varied between cases, with a median of 11 (range, 0–142), in agreement with previous PPGL studies5,6,8,9, and much lower than in other cancer types10. Two of the PPGL studies suggested an association between higher TMB and metastatic behavior5,9. In this study, we confirmed this observation and also detected that TMB increased from non-metastatic primary tumors (median number of variants: 8; range 2–25) to primaries of mPPGLs (17.5; 1–54) to metastases (24; 0–142) (Fig.2b). We also observed a significantly different TMB among genomic subtypes, with the Wnt-altered subtype having the highest values, regardless of tumor behavior (Supplementary Fig.3a, b). Although there is scarce information about MSI in PPGLs11,12,13, it is known that MSI contributes to the tumorigenic hypermutated phenotype in other cancer types14. In our study, we observed differences in the MSI score according to the clinical behavior (non-metastatic tumors median: 0.13, range 0.11–0.18; primary metastatic tumors: 0.16, range 0.10–0.22; and metastases median: 0.18, range 0.11–0.25) (Fig.2c). Moreover, the MSI score and the TMB exhibited a moderately significant correlation (Fig.2d), indicating that the higher TMB observed in metastatic tumors could be linked to a higher MSI status. Both features were associated with risk of metastasis development in univariate and multivariate analyses (including as covariates sex, mutations in Krebs cycle susceptibility genes15andATRX/TERTalterations), suggesting they are independent predictors of metastatic risk (Table1).

The data included correspond to 261 tumor-normal pairs from the CNIO and TCGA cohorts.aOverview of the number of variants, MSI score, and other clinical and genomic features of the whole series. Tumors have been ranked by TMB (calculated with category 5 variants). The legend in the bottom indicates the color code for each item. The filtering strategy for WES events is shown in Supplementary Fig.2. bTMB andcMSI score (extracted with MANTIS) across PPGL tumors (n= 256). P-values were calculated with a two-sided Mann–Whitney–Wilcoxon (MWW) test and significant ones are shown in the figures.dCorrelation between TMB and MSI score. Two-sided Pearson’s correlation coefficient is shown.e, fProgression-free survival analysis of patients according to primary tumors’ TMB and MSI score, respectively. Only primary tumors from non-metastatic and metastatic patients were analyzed.Higher TMBindicates tumors with values > than the third quartile (n= 37);lower TMBfor the remainder cases (n= 175).Higher MSIscore indicates cases with MSI score >0.15 (n= 40);lower MSIfor the remainder cases (n= 172). Kaplan–Meier plots of time to progression (time between the first PPGL diagnosis and the first documented metastasis) are shown together with P-values calculated using a log-rank test. Median progression time (± standard error) of each group is depicted in the corresponding color. Patients without evidence of metastases were censored at the date of the last follow-up.g–iTMB variation according to tumor volume (cm3) (n= 28), % Ki67 positive cells (n= 16) andMKI67mRNA expression (n= 191). Volumes and % Ki67 cells data, when available, were extracted from the pathological anatomy reports received with each specimen.High MKI67 mRNA expressionindicates expression levels above the 3rd quartile value of the whole cohort andlow MKI67 mRNA expressionwhen levels were beneath the 3rd quartile value.jTMB variation according to the age of the patient at surgery (n= 225). For (h)–(k) a two-sided MWW was applied to test for differences, and, except for (j), only primary tumors (non-metastatic, aggressive and metastatic) were considered.kFrequency ofATRX/TERT-alterations within genomic subtypes. Two-sided Freeman–Halton test was used to test for differences between genomic subtypes. Metastatic tumors include primary tumors and metastases; if paired primary-metastasis is available, only one tumor per patient is represented.lTMB andmMSI score inATRX/TERT-wild-type (WT) and inATRX/TERT-altered tumors (n= 256). Two-sided MWW was used to test for differences. For all box-plots in this figure: the median value is marked, and Tukey whiskers are represented. See also Supplementary Fig.3and Supplementary Table2. Source data are provided as a Source Data file.

To note, both higher TMB and MSI score were also associated with a shorter time to progression (TTP) (Table1; Fig.2e, f; Supplementary Fig.3c). In agreement, we observed a significantly higher TMB in larger tumors, and in those with >5% Ki67 positive cells and highMKI67expression levels (Fig.2g–i); no association was detected between the MSI score and the tumor volume or the % of Ki67 positive cells (Supplementary Fig.3d, e), probably due to the low number of tumors included in this analysis, but a significant association between the MSI score and theMKI67expression levels was identified (Supplementary Fig.3f).

TMB in PPGLs also correlated with the age of the patients at surgery (Fig.2j), feature that has already been shown for many cancers16. In fact, it has been described that the age-dependent increase in somatic variants that occurs naturally in tissues, increases the probability of mutating driver genes, making cancer a disease of ageing17.

Without considering the 10 main somatic PPGL gene drivers represented in our cohort (HRAS, FGFR1, CSDE1, MAML3, IDH1, RET, VHL, NF1, MAX,andEPAS1), ATRXwas found mutated in the tumors of 3.4% patients (Supplementary Table2) and was established as the only gene significantly mutated in metastatic tumors (Fisher’s test,P= 2.9 × 10−4), as previously reported8. SinceTERThas been linked to mPPGL18, we also exploredTERTalterations (overexpression, amplification, promoter mutations and hypermethylation; see methods section) in a subset of cases, and confirmed thatTERTalterations are events associated with mPPGLs (Supplementary Fig.3g; Fisher’s test,P= 8.4 × 10−7). We validated thatATRXandTERTalterations accumulate mainly inSDHB/FHpseudohypoxic tumors19, but also found a high frequency of these events in Wnt-altered tumors (Fig.2k). Moreover,ATRX/TERT-altered samples presented higher TMB and MSI score than non-altered ones (Fig.2l, m), a significantly higherMIK67expression (Supplementary Fig.3h) and a tendency to higher Ki67 positive cells and bigger tumors (Supplementary Fig.3i, j).

Somatic copy-number alteration profile in metastatic pheochromocytoma and paraganglioma

We analyzed the somatic CNA (SCNA) profiles of the CNIO and TCGA cohorts applying FACETS to the germline-tumor matched WES data (Fig.3a). Higher SCNA burden, measured by the number of SCNA events, was observed in more aggressive tumors (Fig.3b), with metastatic samples (primaries and metastases from metastatic patients) exhibiting the highest number of SCNA events (median: 15, range: 1–48 and median: 29, range: 5–65, respectively) in comparison to non-metastatic tumors (median: 8, range: 0–41). In addition, a tendency towards shorter time to progression was observed in patients with tumors with high SCNA burden (Supplementary Fig.4a). Similarly to the TMB, SCNA burden was also higher inATRX/TERT-altered tumors (Fig.3c). However, in the case of SCNA burden, no major differences were observed between tumors of different genomic subtypes (Supplementary Fig.4b).

The data included is from 234 tumor-normal pairs from the CNIO and TCGA cohorts.aOverview of the arm-level SCNA, and clinical, molecular and technical features of the whole series. Tumors have been ranked by SCNA burden (number of SCNA events). The legend in the bottom indicates the color code for each item. Top panel: % of samples with arm-level CNA in metastatic and non-metastatic primary tumors. Significant (Fisher test; FDR < 0.1) regions between both groups are annotated.bSCNA burden (number of SCNA events) across PPGLs (n= 232).Pvalues shown were calculated with a two-sided MWW test.cSCNA inATRX/TERT-wild-type (WT) and inATRX/TERT-altered tumors (n= 232). Two-sided MWW was used to test for differences.dProgression-free survival analysis of patients according to the presence of whole chromosome 5 gains (n= 17) or no whole chromosome gains (n= 177) in primary tumors. Only primary tumors from non-metastatic and metastatic patients included. Kaplan-Meier plot of time to progression (time between the first PPGL diagnosis and the first documented metastasis) and theP-value calculated using a log-rank test is exposed. Patients without evidence of metastases were censored at the date of the last follow-up. For all box-plots in this figure: the median value is marked, and Tukey whiskers are represented. See also Supplementary Fig.4. Source data are provided as a Source Data file.

Genome doubling was observed in 14.5% of primary metastatic tumors whereas only 5% of the non-metastatic tumors duplicated their whole genome (P= 0.023, two-tailed Fisher test). Contrary to that reported by Fishbein et al.5, this tendency was not linked to its genomic subtype (P= 0.38, two-tailed Freeman–Halton test).

The most frequent (>25% of all tumors) arm-level deletions affected 1p, 3p, 3q, 11p, 11q, 17p, 22p, and 22q chromosome arms. Arm-level chromosome gains were less common, mostly involving 1q and 7p, as previously reported for PPGLs6. Specifically, we observed 1p loss and 1q gain inSDHBtumors; 3p, 4q and the whole chromosome 11 loss inVHL-; 6q loss inHRAS-; chromosome 17 loss inNF1-related tumors. Whole chromosome gains were also detected, chromosome 4 inMAML3- and chromosome 10 inRET-cases. Although 11p deletion is a common event in PPGLs, this SCNA was not observed in tumors harboring oncogenicHRASmutations (FDR < 0.05). As already shown5,6, SDHB-, VHL- and NF1-related tumors deleted genomic regions where their specific locus is placed (1p, 3p, and 17q, respectively).MAML3- andRET-related tumors also exhibited gains in chromosomes 4 and 10, comprisingMAML3andRET, respectively.

We identified four altered arm-level SCNA differing between metastatic and non-metastatic primary tumors (Fig.3a). These included gains in 1q, 5p, and 5q, and deletion of 3q. Whole chromosome 5 gain was also associated with shorter TTP (Fig.3d), finding which is in line with the fact thatTERTis located at 5p15.33. This association remained statistically significant after multivariate Cox hazard regression analysis using sex, mutations in Krebs cycle susceptibility genes, TMB and MSI score as covariates (HR = 2.56, 95% CI = 1.15–5.70,P= 0.021), indicating its independent prognostic value.

Transcriptomic landscape of metastatic pheochromocytoma and paraganglioma

We selected a signature of 26 genes differentially expressed between primary tumors of patients with and without metastatic disease (n= 54 andn= 176, respectively) and independent of genomic subtypes (Fig.4a). After expression dichotomization and multivariate analysis using genomic subtype as a covariate, 25 out of the 26 genes were also associated with risk of metastasis (P < 0.005; Fig. 4b) and time to progression (P < 0.05; Fig. 4c). To validate these findings, we re-analyzed an independent cohort of PPGLs20and confirmed a significant differential expression between the two groups for 17 of the 25 selected genes (Supplementary Fig.5a), including 14 that demonstrated a value for stratifying patients according to metastatic risk (C5orf49, DNASE1L3, HMGB3, RRM2, CDT1, TTC9, CCNB2, IQGAP3, FAM83H, CDK1, PBK, NETO2, CSMD2andPMAIP1) (Supplementary Fig.5b). Functional enrichment analysis of this gene signature using STRING v1121indicated a significant overrepresentation (PPI enrichmentP= 3.08 × 10−7) of genes related to mitotic cell cycle, G1/S-specific transcription, p53 signaling and DNA damage response (Supplementary Fig.5c). CDK1, the only gene included in the 4 enriched gene sets, displayed the most significant association with metastatic risk in a combined analysis that included both cohorts (OR = 9.01, 95% CI = 5.17–15.71,P= 8.5 × 10−15; Supplementary Fig.5d; Fig.4d); furthermore, high expression ofCDK1showed a strong association with shorter TTP, regardless of the tumor genomic subtype (HR = 4.15, 95% CI = 2.32–7.24,P= 1.1 × 10−6; Fig.4e). A multivariate Cox hazard regression analysis including the proposed prognostic molecular markers (chromosome 5 gain, higher TMB and MSI score) as covariates, pinpointedCDK1high expression as an additional and independent indicator of prognosis (HR = 2.21, 95% CI = 1.11–4.79,P= 0.045). We also validated by immunohistochemistry (IHC) higher CDK1 protein levels in the nucleus of primary metastatic tumors (Fig.4f), which were correlated withCDK1gene expression (Supplementary Fig.5e).

aGene signature associated with mPPGL. mRNA expression levels of the 26 differential expressed selected genes, tumor behavior, genomic subtype and genotype are depicted in rows; primary tumors appear in columns (n= 230). Univariate (black) and multivariate (blue) logistic and Cox regression analysis of metastasis risk (b) and TTP (c), respectively. Gene expression was dichotomized as follows: for downregulated genes in mPPGLs, median expression was used as threshold (0 – below the median expression level; 1 – above the median expression level); for up-regulated genes in mPPGLs, high expression levels > than the 3rd quartile (0 – below the 3rd quartile threshold value; 1 – above the 3rd quartile threshold value). Multivariate analysis included as covariate genomic subtype. Only data from primary tumors from non-metastatic and metastatic patients were included.dBox plot ofCDK1expression of primary tumors included in this study united to those from an independent cohort20 (n= 417). Expression from both cohorts was z-score transformed (centered at the mean of non-metastatic group for each cohort). The P-value corresponds to a two-sided MWW test.eProgression-free survival analysis of patients according toCDK1expression. Kaplan–Meier plot of time to progression (time between the first PPGL diagnosis and the first documented metastasis) is shown together with theP-value calculated using a log-rank test. High levels (above the 3rd quartile value of the whole cohort) are represented in red (n= 71) and low expression (below the 3rd quartile) in blue (n= 172). Patients without evidence of metastases were censored at the date of the last follow-up.fRepresentative images (right) and quantification (left) of CDK1 IHC staining in a subset ofn= 41 PPGLs. Three PPGLs classified as aggressive were not included in this analysis. Scale bar in images = 100 μm. Unpaired two-sided t test was used to test differences between groups.gSignaling networks underlying mPPGLs. Plots show normalized enrichment score (NES) of significantly enriched gene sets (FDR < 0.05, GSEA) from MSigDB, grouped according to relevant biological processes. Each dot represents a gene set and is highlighted in the specific color for each process.hClose-up of gene sets gathered inimmune responseannotation and differentiation betweeninterferon signaling, T cells/activation differentiationandanti-tumor related cytokines. For all box-plots in this figure: the median value is marked, and Tukey whiskers are represented. See also Supplementary Figs.5, 6. Source data are provided as a Source Data file.

In order to identify altered processes in metastatic tumors, and not only specific genes, we performed a GSEA analysis comparing metastaticvsnon-metastatic primary tumors (Fig.4g; Supplementary Fig.6). This analysis revealed a positive normalized enrichment score (NES > 1.85; FDR < 0.05) not only in gene sets involved in cell cycle, DNA repair and p53 pathway, but also in those related to extracellular matrix (ECM) organization, translation/ ribosomes, ubiquitin/ proteasomes, and motility; this indicates that genes belonging to these categories are mainly upregulated in metastatic primary tumors. In contrast, gene networks associated with cilium formation, ion transport, Rho-GTPases, adhesion, nervous system (NS) development/ neuron projection and immune response displayed negative NES (NES < −1.85; FDR < 0.05), illustrating a downregulation of genes in these categories in metastatic primary tumors. Cilia loss and Rho pathway deregulation have already been reported in PPGLs9,22,23, reinforcing the importance of these signaling networks in PPGL malignant transformation. Signaling cascades related to the immune system and deregulated in metastaticvsnon-metastatic primary tumors (Fig.4h; Supplementary Fig.6o) included multiple sets involved in innate and adaptive immunity, cytokine secretion and inflammation that had not previously been associated with mPPGLs.

Metastasis risk classifier for pheochromocytoma and paraganglioma

We evaluated the classification power of the potential markers of mPPGL described in this study: high MSI score, high TMB, gain of chr 5 and highCDK1expression. We also aimed to assess their power in combination with the already described markers of mPPGL (germline mutation in Krebs cycle genes,MAML3-fusion,TERT-alt,ATRX-mutation and highMKI67expression5,8,15,19,24). The area under the receiver operating characteristic curves (AUC under ROC) of each event individually and all combinations possible (n= 511) were computed. The AUC, 95% confidence interval (CI) and sensitivity/specificity per each classifier are in Supplementary Data1. The best classifier is the one that takes into accountATRX-mutations, high MSI score, highCDK1expression andMAML3-fusions, showing a 100% sensitivity with the highest ability to predict mPPGL (AUC = 0.902, 95%CI = 0.855–0.948) (Supplementary Fig.7).

Functional enrichment analysis of mutated and CN-altered genes in metastatic pheochromocytoma and paraganglioma

Panther tool25was applied to detect significantly enriched gene sets among the list of genes with high impact consequence mutations in primary tumors from metastatic patients, and significantly differing focal SCNA between non-metastatic and metastatic primary tumors. The list of high impact mutated genes comprised 323 genes that grouped into 52 significantly enriched processes (fold-enrichment>1.8 and FDR < 0.05) classified into five different functional annotations:circulatory system development/morphogenesis, ECM organization, cell adhesion and motility, actin cytoskeleton, andNS development and neuron projection(Supplementary Fig.8a). Focal SCNA significantly differing between non-metastatic and metastatic primary tumors and resulting in altered gene expression was present in 911 genes (see methods section). After applying a functional enrichment analysis using the 911 gene list, we recognized 120 significantly enriched gene sets (fold-enrichment>1.8 and FDR < 0.05) (Supplementary Fig.9a).

Given the availability of RNA-Seq data, we executed GSEA with the 52 and 120 sets that resulted from the previous analysis in order to elucidate which categories were affected at the mRNA expression level. The analysis performed with the mutated gene sets (52) revealed differential expression (comparing metastaticvsnon-metastatic primary tumors) of genes involved in five described functional annotations (Supplementary Fig.8b), three of which were more relevant at the transcriptional level (>30% of sets of particular annotation with FDR < 0.05 in the GSEA analysis):ECM organization, cell adhesion and motility, andNS development and neuron projection(Fig.5a). Notably, ECM and adhesion (Supplementary Fig.6l, m) are terms classically linked to tumor progression and metastasis26,27. The weaker relative activity towards neuron differentiation (Supplementary Fig.6k) may indicate a dedifferentiated neuronal-like phenotype in mPPGLs. This feature could be related to a neuroendocrine-to-mesenchymal transition associated with aggressive traits28,29,30.

Functional enrichment analysis of mutated (a) and CN-altered (b) genes in mPPGLs.y-axis indicates the −log10(FDR) of each gene set identified in the functional enrichment analysis performed using the list of 323 genes mutated in the metastatic primary tumors (a) or the 911 genes with differing SCNA between metastatic and non-metastatic primary tumors (b). Only those gene sets with fold-change enrichment (FCe) > 1.8 and FDR < 0.05 by Fisher’s Exact test were considered.x-axis indicates −log10(FDR) × NES sign from the GSEA analysis performed using the gene sets recognized in the functional enrichment analysis and the RNA-Seq series of metastaticversusnon-metastatic primary tumors. The gene sets related to each functional annotation are shown in different colors as depicted in the legend. Dashed lines indicate the cut-off values for significance. The extended version of the figure, in which each gene set name is annotated, is depicted in Supplementary Figs.8, 9. Source data are provided as a Source Data file.

SCNA analysis identified gene sets enriched in mPPGLs, 98 exhibiting statistically significant differences between non-metastatic and metastatic primary tumors (Supplementary Fig.9b; Fig.5b). These could be classified into 5 different functional annotations:translation and ribosomes, DNA repair and TP53 pathway, splicing, cell cycle and chromatin regulationandubiquitin and proteasome. Ribosome biogenesis has been reported to be essential for tumorigenesis and sufficient to drive malignant transformation31. Ribosome biogenesis, by preventing p53 ubiquitination and degradation32,33, is also linked to the cell cycle34. Moreover, hyperactive rDNA transcription leads to DNA damage and genome instability, activating DNA damage responses31, that could also explain the higher TMB, MSI score and SCNA burden observed in mPPGLs (Figs.2, 3).

Therefore, integration of transcriptomic and genomic data (mutations and SCNA), suggests causal mechanisms underlying 7 out of the 15 functional annotations deregulated at the transcriptional level (Fig.4g).

Immune landscape of metastatic pheochromocytoma and paraganglioma

Using a combination of gene expression analyses and IHC scoring with a set of immune markers, we characterized the immune infiltration profile of each tumor. First, we used the data generated by Thorsson et al.35for the TCGA cohort and applied theCRI-iAtlas/ImmuneSubtypeClassifierR package to classify the tumors from the CNIO series according to the available RNA-Seq data; from this, we established the distribution of the six immune subtypes described by Thorsson et al. (C1–C6) across our collection (with 0.38% tumors in C1, 48.66% in C3, 45.21% in C4, and 5.75% in C5). C3 (inflammatory) and C4 (lymphocyte depleted) were the predominant immune subtypes in PPGLs, but interestingly they were represented differently according to tumor type; C3 was more abundant in non-metastatic tumors (52.3%), whereas C4 showed more prevalence among primary metastatic tumors and metastases (50.9 and 71.4%, respectively) (Fig.6a). In addition, different proportions were also observed among genomic subtypes: C3 was more abundant in the pseudohypoxic subtype (72.2%), while C4 was more frequent in the kinase signaling subtype (67.6%). Of note, regardless of genomic subtype, more metastatic primary tumors and metastases were classified as C4 (P < 1 × 10−4), and those patients with tumors categorized as C4 seemed to have a shorter TTP than those classified as C3 (Supplementary Fig.10a).

aClassification of immune subtypes within tumor tissue types and genomic subtype. The proportion (%) of samples of each immune subtype per tumor type and genomic subtype is shown. Aggressive tumors were excluded from pie charts due to the low number of samples available.bPlot showing the −log10(P) resulting from a two-sided MWW analysis to define differences between metastatic primaryversusnon-metastatic tumors in the different immune cell classes estimated with CIBERSORTx. The ‘sign(effect)’ indicates the direction of the fold-change between the proportions in both groups. Cell types with >85% of the samples with ‘0 s’ were excluded from the analysis. Columns that surpass the red dashed line haveP < 0.05. The top row summarizes univariate logistic regression analysis comparing immune cell proportions and metastatic risk. The color of the cell is relative to the −log10(OR), and * indicatesP < 0.05. cPercentage of CD8+ T cells infiltrated among tumor cells detected by immunohistochemistry (left) in a subset ofn= 39 PPGLs with different clinical behavior. Three PPGLs classified as aggressive were not included in this analysis, and for two cases IHC were not assessable. Median ± IQR is shown. Two-sided MWW was applied to test for differences. Representative images (right) of CD8 IHC. Scale bar in images = 200 μm.dTop panel: median enrichment score of the different Fges in metastatic and non-metastatic primary tumors. Two-sided MWW was applied to test for differences between metastatic (n= 55) and non-metastatic (n= 176) primary tumors; significantPvalues are shown. Univariate (black) and multivariate (blue) logistic and Cox regression analysis of metastasis risk (middle panel) and TTP (bottom panel), respectively. Enrichment scores were used as a continuous variable. Multivariate analysis included genomic subtype as covariate. Only data from primary tumors from non-metastatic and metastatic patients were included. Significant associations in Fges scores after multivariate analysis are shaded in blue.eKaplan–Meier plots of time to progression in patients according to different immune cell type levels found in primary tumors. Only primary tumors from non-metastatic and metastatic patients included. High levels (above the median level of the whole group) are represented in red (n= 92,n= 113, andn= 107, respectively for NK resting, NK activated and T regulatory) and low expression (below the median level) in blue (n= 130,n= 109 andn= 115, respectively for NK resting, NK activated and T regulatory).P-values shown inside the plots were calculated using a log-rank test. Patients without evidence of metastases were censored at the date of the last follow-up. See Supplementary Fig.10for the extended version. Source data are provided as a Source Data file.

Next, we used the immune score, extracted by ESTIMATE36, to evaluate the presence of immune cells in the tumors using RNA-Seq data. Although we did not observe different infiltration levels between the different tumor types (Supplementary Fig.10b), we discovered substantial variations in the proportion of immune cells extracted by deconvolution of the RNA-Seq data using CIBERSORTx37(Fig.6b). We observed a significantly higher number of resting and lower number of activated CD4+ T cells, lower proportion of CD8+ T cells, and a higher number of resting and lower quantity of activated NK cells in metastatic primaries compared to non-metastatic tumors. The lower levels of CD8+ T cells in metastatic tumors was validated by IHC (Fig.6c). The suppression of activated T cells is supported by the GSEA analysis with 24 terms related to activation/differentiation/proliferation identified among the negative NES (Fig.4h; Supplementary Fig.6o, blue arrows). This notion was reinforced by the differences between non-metastatic and metastatic primary tumors found in the functional gene expression signatures (Fges) enrichment scores, calculated by GSVA as in Bagaev et al.38, related to effector cell traffic and NK cells (Fig.6d, top panel). We also noted a higher representation of macrophages M2-like, known by their anti-inflammatory and pro-tumoral role39, in metastatic primary tumors. This is in agreement with the observation that a higher proportion of metastatic tumors were classified into C4 subtype, characterized by high M2-like macrophage response35.

Several gene sets (n= 9) related to interferon signaling also had negative NES in the GSEA analysis (Fig.4h; Supplementary Fig.6o, green arrows), and a significantly lowerTh1 signature- andanti-tumor cytokines-, and higherpro-tumor cytokines-Fges scores were found in metastatic primary tumors (Fig.6d, top panel), both supporting the immunosuppressive microenvironment encountered in the metastatic primary tumors. These differences increased when metastases were incorporated into the analysis (Supplementary Fig.10c). Other cytokines with anti-tumor immunity properties, such as TNF, IL6, IL2, and IL140, were also identified by the GSEA analysis (Fig.4h; Supplementary Fig.6o, purple arrows).

Interestingly, lowerNK cells-, effector cell traffic-, Th1 signature- andantitumor cytokines-Fges scores were associated with risk of metastasis development (multivariate analysis adjusting for genomic subtype; Fig.6d, middle panel), suggesting that TME is independent of genotype.

The latter Fges scores were also associated with TTP (Fig.6d, bottom panel). We also found that some infiltrating immune cell types estimated with CIBERSORTx were significantly associated with TTP; high levels of resting and low levels of activated NK cells were associated with shorter TTP (Fig.6e; Supplementary Fig.10d). Similarly, high levels of regulatory T cells (Tregs) were also associated with shorter TTP.

Immunogenomics as a theranostic tool in the immunotherapy contexture of metastatic pheochromocytoma and paraganglioma

Bagaev et al. recently proposed four Fges-based pan-cancer microenvironment subtypes that correlate with patient response to immunotherapy in multiple cancers38. To classify our tumors according to their TME, we applied k-means clustering using the 29 Fges scores. We defined four TME subtypes in our cohort (Fig.7a); among these, the immune-enriched non-fibrotic was the least represented (14.3% of the cases), and the subtype with fewer metastatic tumors (Supplementary Fig.11a) and with better prognosis (Fig.7b). TME subtypes were strongly correlated with the ImmuneScore and StromalScore, obtained from ESTIMATE36, and predicted levels of infiltrating immune and stromal cells, respectively (P= 4 × 10−38and 5.4 × 10−44, Krustal–Wallis test) (Fig.7a). We also observed a trend towards TME subtypes being associated to neoantigen load and TMB (Supplementary Fig.11b, c). TME subtype proportions also differed between immune subtypes (Supplementary Fig.11d) and PPGL genomic subtypes (Supplementary Fig.11e). The pseudohypoxic cluster was mainly enriched in the fibrotic subtype, characterized by high vascularity and an immunosuppressive phenotype38as recently reported by Celada et al.41; in contrast, the kinase signaling cluster was highly represented in the immune-enriched non-fibrotic subtype, reflecting the influence of the genomic subtypes on TME. These differences were also evident when comparing Fges scores between genomic subtypes (Supplementary Fig.11f). We observed higher representation of a Th1 signature and checkpoint molecules in the kinase signaling tumors, and the pseudohypoxic tumors were characterized by higher signature Fges scores related to angiogenesis and fibrosis (cancer-associated fibroblasts, matrix remodeling, endothelium), as already reported42,43.

aHeatmap of the 267 PPGL tumors profiled by RNA-Seq and classified into the four distinct TME subtypes described Bagaev et al. based on unsupervised k-means clustering of the 29 Fge enrichment scores. Genomic and clinical features are exhibited as depicted in the legend.bKaplan–Meier plot of time to progression in patients according to the primary tumor TME subtype (n= 33 for IE,n= 55 for F,n= 74 for IE/F andn= 62 for D). Only primary tumors from non-metastatic and metastatic patients included.P-value was calculated using a log-rank test. Patients without evidence of metastases were censored at the date of the last follow-up.cExpression (median normalized z-score expression) of main immunoregulators (IMs) according to the clinical behavior and genomic subtype. Two-sided MWW and two-sided Krustal-Wallis tests were applied to test for differences between metastatic (n= 55) and non-metastatic (n= 176) primary tumors and between the different genomic subtypes (pseudohypoxic,n= 33; kinase signaling,n= 16; Wnt-altered,n= 6) within metastatic primary tumors, respectively; both tests were applied on the normalized data. SignificantP-values are shown.dRepresentative images (top) and quantification (bottom) of PD-L1 IHC staining in a subset ofn= 44 PPGLs. Quantification is represented for each genomic subtype. Two-sided Freeman–Halton test was used to test for differences between groups. OnlyMAML3class exhibited statistically significant differences compared to the other groups. Scale bar in images = 150 μm.eNeoantigen load across genomic subtypes in primary tumors (n= 170). Only tumors with WES and RNA-Seq data available were included in the analysis; WT Wwnt-altered tumors have not been represented. The median value is marked and Tukey whiskers are represented.P-values shown in the figure correspond to two-sided MWW test. OnlyMAML3class exhibited statistically significant differences compared to the other groups.fPercentages of CD8+ T cells infiltrated among the tumor cells detected by immunohistochemistry in a subset ofn= 42 PPGLs with different genomic subtypes. For two cases, IHC was not assessable. Median ± IQR is shown. MWW was applied to test for differences between groups. See also Supplementary Figs.11–13. Source data are provided as a Source Data file.

In order to further explore the immunophenotyping of patients with mPPGL, we examined the expression of key immunomodulators (IMs) targeted in clinical oncology by several IM agonists and antagonists44. Gene expression of IMs (Fig.7c; Supplementary Fig.12a) varied between metastatic and non-metastatic tumors and across genomic subtypes, possibly suggestive of their role in delineating the immune infiltration of the tumors. Genes with the greatest differences (P < 0.01) between metastatic and non-metastatic tumors were CD274, CD276, TGFB1, VEGFA, CD40andTLR4. CD274, which encodes for programmed death-ligand 1 (PD-L1), showed the lowest levels in pseudohypoxic tumors and the highest in Wnt-altered subtype, including several metastaticMAML3-tumors. PD-L1 IHC in an independent series ofn= 44 PPGLs (Supplementary Table3) confirmed that 34% (15/44) of PPGLs, 7 harboring aMAML3-fusion, exhibit PD-L1 staining in tumor cell membranes. These represented 100% ofMAML3-tumors included in the series and, with the exception of oneSDHBtumor, were the only tumors exhibiting highly positive PD-L1 staining (Fig.7d); no association was observed according to clinical behavior (Supplementary Fig.12b). A higher neoantigen load was also observed inMAML3-tumors compared to all other tumors (MWW test:P= 0.0067; Fig.7e).

Given the differences observed in the aforementioned key immune players, further investigation is needed in order to elucidate the differences in terms of immune infiltrates between the tumors classified in the different genomic subtypes. We observed differences between molecular clusters according to Fge scores, immune and TME classifications (Supplementary Fig.10c, d; Fig.6a), and, as anticipated, we also noticed differences according to relative levels of some immune cell types among genomic subtypes (Supplementary Fig.13, top). One difference was the greater infiltration of CD8+ T cells inMAML3-tumors (Fig.7f). On the other hand, we show that metastatic tumors present a similar pattern in the proportion of major classes of immune cells across the different genomic subtypes (Supplementary Fig.13, bottom), reinforcing the notion that some immune cell populations may be partially linked to the molecular genomic subtype, but others to the tumor behavior.

Discussion

The comprehensive genomic analysis of mPPGL described here highlights not only the heterogeneous landscape of secondary genomic alterations, but also the presence of shared TME features within genomic subtypes and/or according to tumor behavior.

In this study, we reveal that TMB, MSI status and SCNA burden of PPGLs increase with disease progression, suggesting their use as markers of metastatic risk at diagnosis of a primary tumor. We confirm differential TMB in PPGLs as previously reported5,9, and report the association between changes in MSI status and clinical behavior of PPGLs. Although we did not identify any PPGL tumor surpassing the cutoff value used for MSI-high classification in tumors with high TMB10,13, we propose that MSI thresholds may not be applied similarly to tumors with lower TMB, such as for PPGL or acute myeloid leukemia45,46. Frequent gain of chromosomal material had previously been observed in mPPGL47, but here we have linked an increased SCNA burden to mPPGL, mimicking observations in other tumors48. Importantly, TMB, MSI status and SCNA burden were associated withATRX/TERTalterations, reinforcing the notion that alterations in immortalization-related genes contribute to genomic instability49. In fact, different copy number signatures could underlie this heterogeneous genomic instability, as reported for ovarian cancer50.

ATRX/TERTalterations are validated as the only recurrent secondary prognostic-associated events identified in PPGLs, accumulating not only in pseudohypoxic tumors8,19,51,52, but also in Wnt-altered tumors. These findings bolster the concept that alterations contributing to immortalization mechanisms are those that confer a poor prognosis for patients with PPGL. Except forATRXmutations, TERTalterations and chromosome 5 gains, no other genetic alterations could be identified as metastasis priming events. Therefore, and even though PPGL tumor formation is clearly driven by germline or somatic mutations, our data suggest that the molecular switches involved in metastatic progression are likely variable and may involve non-genetic events. However, of course, future efforts should also be directed towards performing WGS to capture potential alterations in non-coding regions or large rearrangements53.

At the transcriptional level, we identified elevated CDK1 as a prognostic factor that, independently of genomic subtype, adds to the higher TMB, higher MSI score and chromosome 5 gains just described herein. Moreover, given that CDK1 is indispensable for driving the cell cycle54, and because the cell cycle is one of the major deregulated signaling networks in mPPGLs (described herein and recently by Zethoven et al.55), there may be potential for developing CDK1-specific inhibitors for treatment56.

Given what has been reported to date57, it will be difficult to identify a single molecular marker capable of accurately predicting the metastatic behavior of these tumors. Therefore, it is of particular interest that our study describes a classifier that takes into account only 4 events (ATRX-mutations, high MSI score, high CDK1 expression and MAML3-fusions) and achieves a sensitivity of 100%. Although further validation will be necessary, we have defined a classifier that could identify all patients at risk of metastasis at the time of diagnosis of the primary tumor.

The immune contexture is being evaluated as a prognostic tool in many tumor types58, but little has been reported in PPGLs41,55,59,60. Here, we identified a high number of downregulated gene sets related to the immune response. By interrogating the immune content, we also showed a complex and heterogeneous ecosystem at the TME level, characterized by the presence of CD8+ T cells, NK cells and a Th1 profile; this was associated with a good clinical outcome, as expected, due to their anti-tumor immune characteristics61. Additionally, we identified an association between high levels of macrophages M2 and Tregs, both cell types known by their anti-inflammatory and pro-tumoral role39,62, and metastatic behavior and poorer prognosis. The excessive catecholamine levels present in metastatic tumors could contribute to the establishment of the immunosuppressive TME63, by diminishing T and NK cell activation64 and by driving the M1- to M2-like macrophage polarization65. Our data demonstrate the potential usefulness of including immune parameters as prognostic classification tools for PPGL.


Cancer immunotherapy is in the spotlight for the treatment of numerous cancers and the TME is known to be important for tumor progression and anti-cancer therapeutic responses66. Therefore, it is possible that classification of PPGLs according to their TME subtype38 could help identify patients who might benefit from immunotherapeutic agents. We established that tumors with immune-enriched non-fibrotic subtype, showing an immune-active TME, exhibited the best prognosis and included mainly non-mPPGLs. In melanoma, bladder and gastric cancer, responders to immune checkpoint inhibitor have the highest proportion of immune-enriched non-fibrotic subtype38. However, what we observe in mPPGL is discouraging. Specifically, the fibrotic subtype, which in our cohort was enriched in mPPGLs, seems to be the least responsive to the aforementioned therapies38.

Strategies to individualize immunotherapy based on patient- and tumor-specific characteristics are a current focus of the scientific community67. With this in mind, we identified PD-L1 positivity, higher TMB and neoantigen load, and CD8+ T cell infiltration in MAML3-tumors compared to other PPGLs. This is noteworthy because all these traits are associated with improved responses to PD1 axis blockade in diverse tumors68,69; this suggests that PD-1/PD-L1 inhibitors may offer a potential treatment for patients with MAML3-tumors. Pembrolizumab, a PD-1 inhibitor displayed modest therapeutic efficacy in a clinical trial of mPPGLs70; although no association with genotype was observed, MAML3 status was not considered, so it remains possible that this could be used to stratify patients. This concept is supported by findings that gene fusions are a source of immunogenic neo-antigens that can mediate responses to immunotherapy, even with a low TMB and minimal immune infiltration71. Of relevance, it is known that WNT/β-catenin signaling, upregulated in MAML3-tumors5, correlates with immune exclusion72. The latter is caused by an increased PD-L1 transcription after MYC or β-catenin binding to its promoter73,74. Therefore, these observations suggest that WNT inhibitors may synergize with PD1 axis blockade in MAML3-tumors.

We could speculate that CD40 and TLR4 low expression identified in mPPGLs could pave the way to propose other approaches to expand the range of targetable immune checkpoints, as in situ vaccination with double CD40-TLR4 stimulation75 to turn “cold tumors” into “hot tumors” and restore PD-1 sensitivity. Notably, CD40 and/or various TLR agonists, currently in clinical trials, have recently been tested in a PPGL mouse model with promising results76.

The notion that the tumor immune portrait is orchestrated by the genomic subtype and that is modulated by the tumor behavior is supported by reported examples that show (i) a link between tumor genotype and immunological architecture, and (ii) the dynamic state of TME during tumor progression77.

In summary, this study serves as a valuable resource for future investigations in the field of PPGL genomics due to the large series of tumors with WES and RNA-Seq data associated with pathological, IHC and curated clinical information. Moreover, the present report contributes with significant data to advance the genomic and immune understanding of mPPGLs, and thereby potentially increase the therapeutic strategies needed for a bench-to-bedside clinical success.










Article classification: Biological abstract
Share to:
Tel:+86-027-87610298
Tel:+86-027-87610297
mail: info@enlibio.com   mail: sale@enlibio.com
mail: alex@enlibio.com   mail: order@enlibio.com
Add:Room A11-329, 1st Floor, No.1, SBI Venture Street, Optics Valley, East Lake
New Technology Development Zone, Wuhan, China.
Certificate NO.:U18Q28010569R0S
备案图标鄂ICP备18027482号-2
©Copyright 2008-2019. All Rights Reserved. Xinqidi,Your biological research partner!
鄂ICP备18027482号-2 | 本站支持  --