Introduction

Pheochromocytomas and paragangliomas (PPGLs) are rare, highly heterogeneous neuroendocrine tumors that pose significant challenges in prognosis assessment due to their complex molecular landscape and variable clinical behavior [1, 2]. Postoperative recurrence and metastasis rates for PPGLs vary widely, ranging from 6.0 to 17.4% [3], complicating treatment and management strategies. An accurate assessment of the risk for recurrence and metastasis following surgery is essential for optimal management, as it directly informs decisions regarding postoperative surveillance and therapy.

Current genomic subtypes, such as pseudohypoxic, kinase signaling, and Wnt-altered subtypes, have shown promise in identifying the molecular characteristics of PPGLs [4]. However, these approaches have limited applicability, as the overall germline mutation rate in PPGL patients is only 30-40% [5, 6] leaving a substantial portion of PPGLs without actionable genomic markers. This creates a gap in our ability to predict recurrence and metastasis based solely on genetic profiles. While some genotype-phenotype correlations have been established [7-9], molecular subclassifications have advanced our understanding of PPGLs biology, the role of secondary molecular events, particularly in the tumor microenvironment (TME), remains underexplored. The TME has been shown to influence prognosis and response to treatment in other tumor types, and emerging evidence suggests that it may play a critical role in the development and progression of PPGLs. However, its precise contribution to disease progression in PPGLs has yet to be fully elucidated.

This study aims to bridge these gaps by analyzing multiple PPGLs cohorts with a focus on the transcriptional landscape and immune microenvironment. Using weighted gene co-expression network analysis (WGCNA) [10], we seek to define the transcriptional subtypes of PPGLs, explore how these subtypes are linked to immune infiltration patterns, and identify key transcriptional alterations associated with recurrence and metastasis. We also aim to propose potential prognostic markers and assess the impact of immune microenvironment profiles on disease progression. This comprehensive analysis will provide new insights into the biology of PPGLs and improve the predictive accuracy of postoperative recurrence and metastasis risk.

Results

Characterization of Three Transcriptional Subtypes in PPGLs via WGCNA

In this study, we employed WGCNA to explore the transcriptional landscape of PPGLs, successfully identifying three distinct transcriptional subtypes that are associated with varying clinical outcomes and may affect disease progression through mechanisms such as immune infiltration. Our analysis encompassed a comprehensive cohort of 87 PPGLs patients from the First hospital of Jilin University (PPGLs cohort), supplemented by an extensive dataset from The Cancer Genome Atlas (TCGA) that included 179 samples (TCGA cohort), which allowed us to validate our findings on a broader scale. We utilized Immunohistochemistry (IHC), RNA sequencing, and whole-genome sequencing (WGS) to analyze the samples, integrating molecular profiles with clinical data to predict recurrence and metastasis based on the identified transcriptional subtypes (Fig. 1A). WGCNA revealed three major transcriptional subtypes within the PPGL cohorts, labeled as Cluster 1 (C1), Cluster 2 (C2), and Cluster 3 (C3). Functional enrichment analysis (Fig. 1B) further validated the biological significance of these transcriptional subtypes. Subtype C1 is enriched in pathways associated with hypoxia and HIF-1 signaling, which are vital for tumor survival and aggressiveness [11]. Subtype C3’s enrichment in metabolic pathways, such as fatty acid metabolism and parathyroid metabolism, reflects its adaptation to metabolic stress, potentially revealing targets for metabolic intervention.

Figure 1.

A detailed analysis of gene expression within each subtype revealed several critical marker genes that characterize each group (Fig. 1C). Subtype C1 is marked by genes involved in HIF-1 signaling and angiogenesis, such as ANGPT2 and CAV1 [12, 13]. In contrast, subtype C2 is defined by genes, including PCSK1N and FTL, while subtype C3 shows a significant upregulation of genes like DYNLL2 and GPX3. The neuroendocrine (NE) features, assessed through a NE score, reveal distinct patterns among the three PPGL subtypes (Fig. 1D). Notably, subtype C2 demonstrates significantly elevated NE score in comparison to subtypes C1 and C3, suggesting a strong neuroendocrine phenotype. Subtype C1 exhibited higher cell cycle scores (Fig. 1E), epithelial-mesenchymal transition (EMT) scores (Fig. 1F), and hypoxia-inducible factor (HIF) signaling scores (Fig. 1G) compared to subtypes C2 and C3. This suggests that C1 may promote rapid tumor growth, metastasis, and potentially lead to a poorer prognosis. Subtype C2, on the other hand, shows an increased inflamed score compared to subtypes C1 and C3, suggesting a more inflamed tumor microenvironment that may be more responsive to immunotherapy (Fig. 1H).

In summary, our WGCNA-based analysis of PPGLs reveals distinct transcriptional subtypes with unique molecular features, highlighting the disease’s heterogeneity and providing insights into the underlying mechanisms that may drive varied clinical behaviors.

Genetic Mutation Patterns and Pathway Alterations Across Transcriptional Subtypes of PPGLs

To elucidate the relationship between transcriptional subtypes and previous genetic mutation classifications, we utilized genetic variation data from 179 PPGLs patients in TCGA cohort. The mutation profile for each subtype highlights distinct genetic alterations in PPGLs (Fig. 2A). Subtype C1 exhibits a high frequency of mutations in genes associated with the pseudohypoxia pathway, such as VHL and SDHB [14], suggesting a unique adaptation to hypoxic conditions that may drive tumor progression. In contrast, subtype C2 predominantly features mutations in genes involved in kinase signaling pathways, including NF1 and RET [15], which imply active signaling networks that could influence tumor growth and response to therapies. Subtype C3 does not exhibit a clearly defined pathway but presents a diverse mutational landscape, indicating the complexity of tumor biology in this group.

Figure 2.

When analyzing the proportion of genetic alterations in key pathways (Fig. 2B), we found that more than half of the cases in each subtype lacked detectable driver mutations. Among the patients with genetic alterations, subtypes C1 and C3 exhibited a significantly higher proportion of pseudohypoxic signaling alterations compared to C2. In contrast, subtype C2 was more strongly enriched for kinase pathway mutations and WNT alterations. A closer examination of pseudohypoxic mutations revealed a significant concentration of SDHB and VHL mutations in subtype C1, while subtype C3 exhibited a notable concentration of SDHA and SDHD mutations (Fig. 2C). Notably, subtype C2 patients exhibited the highest proportion of kinase pathway involvement, along with the greatest variety of kinase pathway mutations (Fig. 2D).

Mutation features, such as microsatellite instability (MSI), tumor mutational burden (TMB), and mutations in pathways like MYC-associated factor X (MAX), which are commonly used to assess the efficacy of clinical immunotherapy, do not show significant differences among subtypes (Fig. 2E). Additionally, principal component analysis (PCA) effectively distinguished transcriptional subtypes within the PPGLs cohort, further supporting the validity of the subtype classification [16]. In both the Magnus [8] and TCGA cohorts, PCA similarly confirmed that transcriptional subtypes could significantly differentiate PPGLs subgroups, and that transcriptional subtyping outperformed traditional genetic mutation-based classification in distinguishing PPGLs (Fig. 2F). In summary, our study offers a thorough molecular characterization of PPGLs, identifying three distinct subtypes with unique transcriptional features.

Immune Infiltration Profiles Across Transcriptional Subtypes of PPGLs

To elucidate the immune-infiltration characteristics of different transcriptomic subtypes, we applied the xCell [17]algorithm to analyze the RNA-seq data from the PPGLs and TCGA cohorts, revealing tumor microenvironment infiltration features at the transcriptional level. The heatmap illustrates the distribution and abundance of different immune cells across the three subtypes (Fig. 3A). Subtype C2, which showed significantly higher inflamed scores compared to subtypes C1 and C3 (Fig. 1H), exhibited increased levels of CD4⁺ Th1 cells and cytotoxic lymphocytes, such as T cells and natural killer (NK) cells. In contrast, hematopoietic stem cells (HSCs) and endothelial cells (ECs) were found to be reduced. Immunohistochemical (IHC) staining analysis of infiltrated CD8⁺ T cells further supports the quantitative findings from the inflamed score, demonstrating that subtype C2 has more pronounced CD8⁺ T cell infiltration at the tumor site compared to subtypes C1 and C3 (Fig. 3B). These findings reinforce the heatmap data and suggest an active cytotoxic immune response in C2. In further statistical analysis, both the PPGLs cohort and the TCGA cohort showed that the C2 subtype had significantly higher infiltration of cytotoxic lymphocytes and CD4⁺ Th1 cells, while exhibiting low infiltration of immunosuppressive cells (such as HSC), reflecting an activated immune microenvironment. In contrast, the C1 subtype was enriched in HSC but showed the lowest levels of cytotoxic lymphocyte and CD4⁺ Th1 cell infiltration. The immune microenvironment infiltration profile of the C3 subtype was intermediate between the two (Fig. 3C-D). This validation in the TCGA cohort, a larger, independent dataset, underscores the robustness of our findings and their potential relevance across different PPGL populations.

Figure 3.

To further investigate the relationship between transcriptomic subtypes and the presence of tumor-infiltrating cytotoxic lymphocytes (TILs), we applied multiple linear regression models. After adjusting for various potential factors that could influence TIL infiltration, we aimed to determine whether transcriptomic subtypes still have an impact on the degree of TIL infiltration. In the PPGLs cohort, regardless of whether the crude model (Model 1), the age, gender, and primary site-adjusted Model 2, or Model 3 (which further adjusted for Ki-67, SSTR2, and SDHB expression based on Model 2) was used, the results showed that compared to the subtype C1, the subtype C2 exhibited the most significant TILs infiltration, followed by the subtype C3 (Fig. 3E). This conclusion was also validated in the TCGA cohort (Fig. 3F). In summary, our extensive immunological profiling across the transcriptional subtypes of PPGLs reveals distinct immune landscapes, with subtype C2 exhibiting an active immunological profile.

Transcriptional Diversity and Immune Microenvironment Characteristics in PPGL Subtypes

To further clarify the tumor microenvironment characteristics of each transcriptional subtype, we performed single-nucleus RNA sequencing (snRNA-seq) analysis on 32 PPGL patients in Magnus’s cohort [8]. The UMAP analysis successfully categorized the genetic profiles into three distinct subtypes: C1, C2, and C3, while also identifying various cell types, including immune cells, chromaffin, myeloid, and fibroblast cells across the samples (Fig. 4A). This diverse cellular representation highlights the complex biological composition of PPGLs and correlates with specific mutational profiles such as SDHB, VHL, and RET. Consistent with transcriptomic sequencing data from the TCGA cohort (Fig. 2A-D), subtype C1 predominantly harbors VHL and SDHB mutations, while subtype C2 is primarily enriched in mutations related to kinase pathways such as RET and HRAS. Subtype C3 is mainly characterized by mutations in SDHA and SDHD. After extracting the immune cells, we performed further cluster analysis. UMAP visualization revealed that the infiltrating immune cells were categorized into 12 distinct subsets, including CD8⁺ T, CD4⁺ T, NK, B cells, dendritic cells (DC), and others (Fig. 4B). Further examination of immune cell distribution using UMAP in subtypes C1, C2, and C3, compared to normal tissue, reveals distinct clustering patterns of immune cells (Fig. 4C-D). This visualization indicates a higher proportion of cytotoxic lymphocyte immune infiltration, specifically CD8⁺ T cells, in subtype C2 compared to subtypes C1 and C3, with no significant involvement of NK cells. Although no significant differences were observed in the overall macrophage content analysis across the different transcriptional subtypes (Fig. 2C-D), further clustering of macrophages using single-cell sequencing revealed an increased proportion of tumor-suppressing M1 macrophages and a significantly reduced proportion of tumor-promoting tumor-associated macrophages (TAMs) in the C2 subtype compared to the other two subtypes.

Figure 4.

To determine whether the increased CD8⁺ T cell infiltration in subtype C2 is associated with enhanced immunoregulatory function, we analyzed ligand-receptor interactions using CellChat. The results revealed that the interactions between CD8⁺ T cells and themselves, as well as with CD4⁺ T cells and NK cells, were significantly enhanced through chemokine pathways such as CCL-CCR. In contrast, the TGFB-TGFBR ligand-receptor interaction between CD8⁺ T cells and immunosuppressive M2 macrophages and TAM cells was significantly weakened (Fig. 4E). In conclusion, the increased infiltration of CD8⁺ T cells in subtype C2 may contribute to immune activation by interacting with other immune cells in the microenvironment.

Pseudotime analysis was used to elucidate the developmental pathways of tumor cells in PPGLs, highlighting dynamic shifts between transcriptional states, particularly the transition from a naive to an activated state within the tumor microenvironment [18]. Pseudotime analysis reveals that compared with normal adrenal cells, subtypes C2 and C1 follow distinct developmental trajectories, while subtype C3 occupies an intermediate position between the C1 and C2 pathways. This suggests that C3 may represent a transitional state between the two, which helps explain the significant differences observed between C1 and C2 in terms of neuroendocrine characteristics, tumor proliferation, HIF-1 pathway activity, and the tumor immune microenvironment, with C3 consistently lying between them (Fig. 4F).

In summary, our multi-omics analysis illuminates the transcriptional diversity and intricate immune interactions present in PPGLs, providing a comprehensive overview of the genetic and immune profiles that characterize each subtype. These results not only deepen our understanding of PPGL biology but also open avenues for more targeted and effective therapeutic strategies tailored to the distinct features of each subtype.

Subtype-Specific Gene Expression Profiles and Clinical Outcomes in PPGLs

To convert this transcriptional subtype into a clinically applicable classification tool, we initially employed a linear regression model to compare the effect values (β value) of various candidate marker genes (Supplementary Fig. 1) on the transcriptional subtypes. We then identified the genes with the most significant effect values and statistical differences for each subtype as the marker genes for that subtype. Ultimately, we identified the genes ANGPT2, PCSK1N, and GPX3, which are significantly overexpressed and have the most significant β values in subtypes C1, C2, and C3, respectively, as the marker genes for these three subtypes (Fig. 5A and Supplementary Fig. 1). ANGPT2, involvement in angiogenesis and vascular remodeling [12], had significantly higher expression in subtype C1 compared to C2 and C3, which may contribute to the aggressive characteristics of this subtype. PCSK1N, encoding the neuroendocrine protein ProSAAS [19], exhibited notably high expression levels in subtype C2. This is consistent with our finding of significantly elevated NE scores in subtype C2 (Fig. 1D), further suggesting that this subtype possesses the most prominent neuroendocrine characteristics. Additionally, GPX3, a significantly upregulated molecule in the C3 subtype, is involved in the regulation of redox homeostasis [20], which may be closely associated with the active metabolic state in this subtype. Immunohistochemical staining confirmed the differential expression patterns of ANGPT2, PCSK1N, and GPX3 in C1, C2, and C3 subtypes, respectively (Fig. 6B). These findings not only validate the transcriptomic data but also emphasize the potential biological roles these genes may play in the pathophysiology of PPGLs.

Figure 5.

Figure 6.

A Kaplan-Meier survival curve analysis revealed significant differences in patient outcomes based on subtype classification (Fig. 6C). Patients in subtype C1 exhibited markedly worse progression-free survival (PFS) compared to those in subtypes C2 and C3, underscoring the critical role of subtype stratification in predicting clinical outcomes. Additionally, Cox regression analysis quantified the impact of these subtypes on PPGL patient prognosis (Fig. 6D). The models, adjusted for various factors, confirmed that subtype C1 was associated with the highest risk of recurrence and metastasis. In contrast, subtype C2 demonstrated the most favorable PFS, significantly outperforming both subtype C1 and subtype C3.

In conclusion, our analysis highlights the potential of ANGPT2, PCSK1N, and GPX3 as subtype-specific biomarkers, offering valuable insights into the molecular heterogeneity and clinical prognosis of PPGLs.

ANGPT2 as a Key Regulator and Diagnostic Marker in PPGL Subtype C1

ANGPT2 had significantly higher expression in subtype C1 compared to C2 and C3, indicating its potential involvement in angiogenesis and vascular remodeling, which may contribute to the aggressive characteristics of this subtype. The Receiver Operating Characteristic (ROC) curve highlights the diagnostic potential of ANGPT2 in subtype C1, with an area under the curve (AUC) of 0.880 (Fig. 6A), demonstrating its high sensitivity and specificity in distinguishing between the subtypes based on their molecular profiles. The samples were categorized into ANGPT2^high and ANGPT2^low groups based on the cutoff of ANGPT2’s ROC in the PPGLs cohort. It was revealed that the ANGPT2^high group was associated with higher scores in the cell cycle, EMT, and HIF-1 signaling pathways (Fig. 6B-D). This suggests that ANGPT2 may not only serve as a specific biomarker for subtype C1 but also act as a regulator of tumor progression and the tumor’s ability to adapt to hypoxic conditions.

To further validate our findings, we employed the rat PPGL cell line PC12 and used CRISPR-Cas9 technology to knock out (KO) the ANGPT2 gene, which was confirmed by Western blot analysis (Fig. 6E). We then compared the tumor growth curves of ANGPT2 KO (ANGPT2^−/−) and ANGPT2^WT PC12 cells in subcutaneous xenograft models in BALB/c nude mice. The results showed that the in vivo growth of ANGPT2^−/− cells was significantly slower (Fig. 6F-G). RNA-seq analysis confirmed that, after ANGPT2 knockout, the expression scores of key genes involved in EMT, tumor metastasis, proliferation, and hypoxia-related pathways were markedly reduced in the PC12 cell line (Fig. 6H-L). This suggests that ANGPT2, a marker gene for subtype C1, functions as a regulatory molecule in the malignant biological behavior of PPGLs. Targeting this molecule may provide a new therapeutic approach for subtype C1 of PPGLs, which is associated with the highest malignancy and poorest prognosis.

Discussion

The significant biological heterogeneity and wide distribution of PPGLs contribute to substantial differences in patient outcomes, making their classification and risk assessment a critical focus in clinical practice. Traditional classification methods based on anatomical origin and pan-genomic approaches centered on pathogenic mutations have offered insights but fail to fully capture the complexity of these tumors, given the relatively low prevalence of detectable mutations.

Our study addresses these gaps by introducing a novel transcriptional subtype classification framework that effectively distinguishes PPGL patients based on transcriptomic expression. This classification demonstrates significantly improved differentiation capability compared to pan-genomic classification and has been validated across independent cohorts. Immunohistochemical analysis further confirmed the robustness of this classification, revealing subtype-specific clinicopathological features that closely align with their molecular characteristics, underscoring its clinical applicability.

Analysis of genomic characteristics revealed that transcriptional subtypes C1 and C3 were enriched in pseudohypoxic genomic subtypes but exhibited distinct mutation profiles: C1 was primarily associated with SDHB and VHL mutations, while C3 was predominantly linked to SDHA and SDHD mutations. Notably, more than half of the patients across all subtypes lacked detectable pathogenic mutations, yet their transcriptome expression profiles mirrored those of patients with such mutations. This finding highlights a previously overlooked patient population in genomic studies.

Further investigation into the tumor microenvironment revealed substantial differences among subtypes. Subtype C1 exhibited the most immunosuppressive microenvironment, characterized by low infiltration of CD8⁺ T cells and CD4⁺ Th1 cells, high infiltration of HSCs, and an elevated proportion of TAMs. Immune activation pathways were weaker in C1, with enhanced inhibitory interactions involving M2 macrophages and TAMs. In contrast, the C2 subtype demonstrated the most activated and inflammatory tumor microenvironment, with the highest inflammatory scores and robust immune cell interactions. These findings provide a foundation for stratified immunotherapy strategies.

The assessment of postoperative recurrence and metastatic risk in PPGLs remains one of the most critical challenges in current clinical practice [21]. Identifying high-risk patients with precision has become a major focus of research. While several approaches, such as pathology-based scoring systems (e.g., PASS and PAGG scores) and genomic mutation-based evaluation systems [22], have been explored, none have achieved satisfactory outcomes. In this study, subtype C1 was also identified as the poorest prognostic group, exhibiting heightened tumor proliferation, enhanced EMT, and the highest recurrence and metastatic risk. Importantly, we identified ANGPT2 was identified as a key marker gene and regulator driving the aggressive behavior of C1, providing both a novel tool for identifying high-risk patients and a potential therapeutic target.

In summary, this transcriptional classification framework advances our understanding of PPGL heterogeneity and provides a theoretical foundation for designing targeted therapeutic strategies. By integrating molecular and immune features, it enables precision medicine approaches that have the potential to improve outcomes for patients with this challenging cancer type.

Methods

Patient cohorts

We conducted a retrospective analysis of patients diagnosed with PPGLs who underwent surgery at the First Hospital of Jilin University between October 2019 and March 2022. Formalin-fixed paraffin-embedded (FFPE) tumor samples from these patients were retrieved and reevaluated by two pathologists using hematoxylin and eosin (H&E) staining and the original pathological reports. Patients with combined PPGLs or evidence of other tumors were excluded from the study. Clinical, pathological, and demographic data were obtained by reviewing electronic medical records. Follow-up information was gathered through regular visits or telephone interviews. RNA bulk sequencing was performed on FFPE samples from 87 patients with complete clinical data and confirmed PPGL diagnoses, along with 5 FFPE samples of normal adrenal tissue as controls. The study was conducted in accordance with the ethical standards of the Declaration of Helsinki and approved by the Ethics Committee of the First Hospital of Jilin University (23K101-001). Given its retrospective design, the need for informed consent was waived. Patient data were processed and analyzed anonymously to ensure confidentiality.

Cell lines

Rat PPGL cell lines PC-12 was purchased from American Type Culture Collection (ATCC) and was cultured in RPMI-1640 medium. Penicillin/streptomycin (#15140-122; Gibco) and 10% fetal bovine serum (#35-081-CV; Corning) were added to all media. All cell lines were cultured at 37°C with 5% CO2.

Mouse models

For CDX models, 1.0 × 10⁶ PC-12 cell lines were subcutaneously injected into the flank of 6-week-old female BALB/c nude mice (Hafukang Co., Ltd., Beijing, China). Tumor size was measured every 2 days, and animal survival rate was recorded every day. Tumor size was calculated as V = (L × W²)/2, where V is tumor volume, L is the length of the tumor (longer diameter), and W is the width. Mice with tumor size larger than 2.0 cm³ at the longest axis were euthanized for ethical consideration. The Animal Care and Use Committee of First Hospital of Jilin University approved all procedures concerning mouse studies.

RNA extraction and gene expression profiling

Total RNA was successfully extracted from 87 FFPE samples using TRIzol and RNeasy MinElute Cleanup Kit (Invitrogen). RNA purity was assessed using the NanoDrop Spectrophotometer (Thermo Fisher Scientific™, Waltham, USA). RNA integrity and concentration were measured with the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). Subsequently, mRNA libraries were created by using the NEB Next Ultra RNA Library Prep Kit (NEB, Beverly, MA, USA), following the manufacturer’s protocol. Geneplus-2000 sequencing platform (Geneplus, Beijing, China) was utilized to sequence the constructed RNA-seq libraries. The sequencing reads containing adapter sequences and low-quality reads were removed to obtain high-quality reads. Reads passing quality control were aligned to the human genome hs37d5 using STAR software²⁹. Transcript assembly was conducted by using StringTie2.

Immunohistochemical staining

Immunohistochemical (IHC) staining was performed to evaluate the expression of ANGPT2, PCSK1N and GPX3 in formalin-fixed, paraffin-embedded (FFPE) tissue sections. Tissue samples were first baked at 65 °C for 15 minutes. Then deparaffinize the sections in xylene (twice,5 minutes each) and rehydrate through graded ethanol solutions (100%,95%,70%, 5 minutes each). Antigen retrieval was conducted by heating the slides in Tris-EDTA buffer pH 9.0(MVS-0098, MXB) at 95 °C for 2 minutes using a pressure cooker. The slides were then allowed to cool to room temperature and washed with phosphate-buffered saline (PBS).

To block endogenous peroxidase activity, the sections were incubated with 3% hydrogen peroxide solution (BCCK5591, Sigma-Aldrich) for 20 minutes, and then washed with PBS. Non-specific antibody binding was minimized by pre-incubating the slides with 10% goat serum (C0265, Beyotime) for 1 hour at room temperature. The tissue sections were then incubated overnight at 4 °C with primary antibodies targeting ANGPT2, PCSK1N, and GPX3 (proteintech-ANGPT2-24613-1-AP, abnova-PCSK1N-H00027344-M02 and Proteintech-GPX3-13947-1-AP).

After washing with PBS, the sections were incubated with HRP signal enhancement solution (PV-9001, ORIGENE) for 20 minutes and washed, then incubated with enhanced enzyme-labeled goat anti-rabbit IgG polymer (PV-9001, ORIGENE) for 20 minutes at room temperature. For chromogenic detection, the signal was developed using DAB (3,3 ‘-diaminobenzidine) (DAB-4032, MXB) and counterstained with hematoxylin (G1004, Servicebio), differentiation solution (G1039, Servicebio) and bluing reagent (G1040, Servicebio), washed with distilled water. Slides were dehydrated through graded ethanol solutions, cleared in xylene, and mounted using neutral baisam (10004160, SCR).

Stained sections were examined under a light microscope (eclipse e100, Nikon), and images were captured. Evaluate the staining pattern and intensity. Use Fiji image J to measure the percentage of the stained area. Use the statistical software GraphPad Prism to perform T-test.

Weighted gene co-expression network analysis (WGCNA)

The bulk-seq data of 87 patients were analyzed by WGCNA using the R package WGCNA, which was divided into three subtypes: C1, C2, and C3. Then, the clustering dendrogram and gene module heat map of the genes were constructed to observe the significant gene modules of each phenotype. At the same time, the correlation analysis heat map between the modules was constructed. The parameters were set as follows: minModuleSize = 15, softPower=7.

Analysis of Immune Infiltration Using xCell

The immune infiltration in tumor samples was assessed using the R package xCell to predict the relative proportions of 19 immune and stromal cell types within the tumor microenvironment. The final results were visualized using a heatmap generated by the R package pheatmap, displaying significant differences in cell enrichment across various samples. Violin plots were utilized to illustrate the differences in immune and stromal cell distributions among the three tumor subtypes in the PPGLs Cohort and TCGA datasets.

Gene Mutation Analysis

For Gene Mutation Landscape, mutation profiles of each sample of TCGA data for genes related to PCPG, were analyzed using the R package mafools. A waterfall plot of the gene mutation landscape was generated using the oncoplot function. For tumor mutation burden (TMB), microsatellite instability (MSI) and mutant-allele tumor heterogeneity (MATH) analysis, TCGA data was utilized, and the R package maftools was employed to read and calculate the MSI, MATH, and TMB scores. The MSI, MATH, and TMB scores across the three subtypes were compared to observe the mutation patterns in the three subtypes.

Single-cell nuclear transcriptome sequencing(snRNA-seq)

For analysis of snRNA-seq data, single cell nucleus RNA sequencing datasets were downloaded from the European Genome-Phenome Archive (EGA) under accession code EGAS00001005861. The snRNA-seq data was analyzed using the R package Seurat v4.4.1. The scRNA-seq data underwent filtration to incorporate genes expressed in no less than 3 cells, as well as cells that exhibited gene expression between 200 and 5000 genes, with mitochondrial transcripts comprising less than 5%. Given patient-driven clustering observed in the analysis of immune populations, we applied the SCTransform v2 method to correct the batch effects between different patients. To ascertain that SCTransform v2 solely addressed technical discrepancies and not biological variations, we manually compared the resulting clusters and their associated markers to those obtained from a single library cluster analysis. After normalizing the data, dimensionality reduction and clustering were performed, with unsupervised clustering analysis used to identify different cell subpopulations. The predominant cellular subpopulations within the and immune compartments were delineated based on cluster-averaged gene expression levels of the subsequent gene markers: 1) Chromaffin cells, PCPG and ATPG: TH and CHGA, 2) Adrenocortical cells: STAR, 3) ECs: FLT1, 4) SCLCs: SOX10, 5) Fibroblasts: PDGFRB, 6) Myeloid cells: MSR1, 7) T/NK cells: CD2, 8) Mast cells: HDC, 9) B cells: CD79A. For in-depth investigation of immune cell populations, the dataset was narrowed down to exclusively include these specific cell types. Subsequently, the data underwent comprehensive processing from its raw form, following the aforementioned methodology. The DimPlot function was used to visualize UMAP plots for all cell types, infiltrated immune cells, and different subtypes. The ggplot function was utilized to create stacked bar plots of cell type proportions.

Cellchat analysis

CellChat is an R package designed for inference, analysis, and visualization of cell-cell communication from single-cell and spatially resolved transcriptomics. CellChat aims to enable users to identify and interpret cell-cell communication within an easily interpretable framework, with the emphasis of clear, attractive, and interpretable visualizations.

Pseudotime and enrichment analysis (monocle2)

The R package of monocle2 was used to perform the diffusion pseudotime analysis. Normalized data computed in Seurat was directly transferred to monocle2. The differential gene analysis function in monocle2 was used to identify genes with significant changes between C1, C2, C3, and Normal subtypes (q-value < 0.05), which were then utilized as input for temporal ordering of cells along the differentiation trajectory. The final results were subsequently transferred back to Seurat. Then, we performed KEGG and GO enrichment analysis by using clusterProfiler package (v3.18.1).

Survival Analysis

The data for survival analysis was derived from PPGLs Cohort. The R packages survival and survminer were utilized to plot Kaplan-Meier survival curves. The survfit function was applied to fit the survival model. Ultimately, the ggsurvplot was employed to plot the survival curves, depicting the survival status of patients in different subtypes.

Additional information

Author contributions

Yang Liu performed experiments, analyzed data, performed statistical analysis and wrote the manuscript; Xu Yan, Yibo Zhang, Zhenfu Gao, Siyu Shi and Jingyun Chen acquired data; Fengrui Nan performed statistical analysis; Lingyu Li conceived the experiments, performed statistical analysis and wrote the manuscript.

Funding

National Natural Science Foundation of China [grant number 82172690 and 82472799] to Lingyu Li

Data availability

All data generated or analyzed in this study are included in the manuscript and Supplementary Information. Source data are provided with this paper.

Additional files

Supplementary information.