Genomic and immune landscape Of metastatic pheochromocytoma and paraganglioma

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).

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 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).

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.

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.

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.