Molecular Signatures of Cancer Stemness Characterize the Correlations with Prognosis and Immune Landscape and Predict Risk Stratification in Pheochromocytomas and Paragangliomas

TL;DR

This study identifies cancer stemness signatures in pheochromocytomas and paragangliomas, linking them to prognosis, immune features, and treatment response.

Contribution

A novel stemness scoring system is developed to predict prognosis and guide treatment in PPGLs.

Findings

High mRNAsi scores correlate with metastasis and poor prognosis in PPGLs.

A stemness-based risk model achieved an AUC of 0.908 for prognosis prediction.

Low-risk PPGLs show higher immune infiltration and better immunotherapy response.

Abstract

Background: Pheochromocytoma and paragangliomas (PPGLs) caused refractory hypertension in clinics. The sustained risk of local or metastatic recurrences or new tumor development prompted more research on diagnosis, prognosis prediction, and immunotherapy. Method: The tumor stemness is closely related to the heterogeneous growth of tumor, metastasis, and drug-resistance, and mRNA expression-based stemness indices (mRNAsi) could reflect tumor stemness. This was calculated based on OCLR machine learning algorithm and PPGLs patients’ TCGA RNAseq data. The relationship between clinical, molecular, and tumor microenvironment (TME) features and tumor stemness was analyzed through the hub genes that best captured the stem cell characteristics of PPGLs using weighted gene co-expression network analysis (WGCNA), Cox, and LASSO regression analysis. Results: Our study found that metastatic PPGLs had higher mRNAsi scores, suggesting the degree of tumor stemness could affect metastasis and progression. HRAS, CSDE1, NF1, RET, and VHL-mutant subtypes displayed significant difference in stemness expression. Patients were divided into stemness high-score and low-score subtypes. High-score PPGLs displayed the more unfavorable prognosis compared with low-score, associated with their immune-suppressive features, manifested as low macrophages M1 infiltration and downregulated expression of immune checkpoints. Furthermore, from the viewpoint of stemness features, we established a reliable prognostic for PPGLs, which has the highest AUC value (0.908) in the field so far. And this could stratify PPGLs patients into high-risk and low-risk subtypes, showing the significant differences in prognosis, underlying mechanisms correlated with specific molecular alterations, biological processes activation, and TME. Notably, high immune infiltration and tumor neoantigen in low-risk patients and further resulted in more responsive to immunotherapy. Conclusion: We indicated that tumor stemness could act as the potential biomarker for metastasis or prognosis of PPGLs, and integrated multi-data sources, analyzed valuable stemness-related genes, developed and verified a novel stemness scoring system to predict prognosis and guide the choice of treatment strategies.

Figures

Full text

1. Introduction

Pheochromocytomas and paragangliomas (PPGLs) are rare clinical tumors with a prevalence of 500 to 1600 cases per year in the general population. They have a wide range of symptoms and indications, ranging from asymptomatic disease to non-specific symptoms, leading to refractory hypertension and catastrophic hypertensive crisis [1]. Pheochromocytomas are tumors that develop from catecholamine-secreting chromaffin cells in the adrenal medulla, whereas paragangliomas are tumors that arise from extra-adrenal sympathetic or parasympathetic paraganglia [2]. Although the two are located in different anatomical locations, their histological manifestations are highly similar due to their common origin in neural crest cells. Surgical excision could help control the large proportion of tumors. However, patients display a wide range of different clinical courses; nearly 15–20% of patients metastasize to the liver, bones, and lungs, which are distant from the primary tumor. There remains a sustained risk of local or metastatic recurrences or the development of a new tumor [3,4]. In the past decade, immunotherapy has obtained promising results in several types of solid cancers by applying novel strategies that strengthen the natural immune antitumor military. However, immunotherapy for PPGLs is in the early stages of development; we must better understand the mechanisms related to immune system activation.

The cancer stem cells (CSCs) hypothesis postulates that a subpopulation of tumor cells is capable of self-renewal and is responsible for long-term tumor maintenance. Furthermore, numerous studies have demonstrated a connection between the biological properties of tumor stem cells and resistance to chemotherapy and radiation therapy, tumorigenesis, metastasis, tumor self-replication ability, tumor metabolic reprogramming, and tumor immune microenvironment remodeling, suggesting that it could be a promising therapeutic target. Recently, a study focusing on non-Hodgkin’s lymphoma demonstrated that CD47 (a tumor stemness-associated gene) is highly upregulated in the tumor cells of some patients, and is resistant to monoclonal antibody (rituximab) therapy. These tumor cells provide a “do not eat me” signal to defend against being phagocytosed by macrophages. Additionally, tumor stemness is associated with the poor prognosis of some cancers, such as soft tissue sarcomas [5], breast [6], and ovarian cancers [7]. However, few studies focus on the role and impact of CSCs in PPGLs.

In this study, we aim to analyze the CSCs’ characteristics in PPGLs patients from The Cancer Genome Atlas (TCGA) and explore the relationship between cell stemness and tumor immunity. We may develop a stemness scoring system that can be applied to predict prognosis and patient responsiveness to immunotherapy after identifying factors highly correlated with tumor stemness. Additionally, the newly identified predictors of response may help determine the subgroup of patients most likely to respond favorably to this therapy. Finally, new trials are underway, promising to reveal further information on a therapy that marks a major advancement in cancer.

2. Materials and Methods

2.1. Data Collection and Processing

A total of 178 PPGLs patients were analyzed in this study. The mRNA expressions and follow-up information were obtained from the TCGA website (https://portal.gdc.cancer.gov/) (accessed on 20 June 2023). Gene mutation information (somatic or germline mutation) and molecular subtypes (TCGA subtype and miRNA cluster) were collected from the Supplementary Materials in Fishbein et al.’s study [8]. Malta et al. [9] developed a machine learning algorithm named “OCLR” to define stem cell features based on the stem cell data (stemness associate gene expression, methylation, and mutation). The two dependent stemness indices (mDNAsi and mRNAsi) for TCGA-PPGLs patients were obtained from their literature. The data used in this article was obtained from publicly available databases.

2.2. The Clinical and Molecular Features Associated with the Stemness Index (mRNAsi)

Since gene expression data are more accessible than epigenetic data, mRNAsi was chosen to explore their association with phenotypic characteristics. The correlation between clinical features and mRNAsi was constructed using the Beeswarm R package (https://www.rdocumentation.org/packages/beeswarm/versions/0.2.3/, accessed on 25 June 2023), involving gender, age at beginning pathologic symptoms, vital Status, pheochromocytoma or paraganglioma, tumor location, and metastatic disease. The correlation between mutation status or molecular subtypes and mRNAsi was also constructed using the Beeswarm R package 0.2.3.

2.3. Associations Between the Stemness Index and PPGLs Immune Signatures

Based on the gene expression profiles of PPGLs samples, the single-sample Gene Set Enrichment Analysis (ssGSEA) method was used to quantify the enrichment levels of tumor-infiltrating immune cells (TIICs) using R package GSEABase. This provides a set of gene signatures for a 28 human tumor microenvironment (TME) immune cell fractions, such as activated T cells, regulatory T cells, and natural killer (NK) T cells. The PPGLs samples were further divided into two stable clusters (high- and low-immunity clusters) according to the scores of ssGSEA using the unsupervised clustering methods of the ConsensuClusterPlus R package. To better comprehend the connection between the amount of immune invading cells and the mRNAsi, we evaluated the correlation between the mRNAsi and immune score, stromal score, and ESTIMATE score, respectively.

The immune score and stromal score representing immune and invasion of stromal cells at the tumor site, respectively, and the ESTIMATE score representing tumor purity, were calculated using the ESTIMATE R package (v3.6.3, https://www.r-project.org/, accessed on 25 June 2023). The relative proportion of 22 immune cells in PPGLs patients was quantified using the CIBERSORT algorithm. The variations in these three immunological correlation scores, the mRNAsi, and immune-related cell types between the high- and low-immunity clusters were compared. The association between the three immunological correlation scores, immune-related cell types, and mRNAsi was estimated using Pearson correlation analysis.

The median value of mRNAsi was chosen as the boundary to divide TCGA-PPGLs patients into the high- and the low-mRNAsi groups. The expressive heterogeneity of immunological checkpoints between the high- and low-mRNAsi groups was determined using unpaired Student t-tests.

2.4. Functional Enrichment Analysis

The DEGs were determined with fold changes > 1.5 and a statistically significant corresponding p < 0.05 using the limma R program. The Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted on these DEGs to reveal their potential function. The GOplot R package was used for the z-score calculation (z-score > 0, referring to positive regulation; z-score < 0, referring to negative regulation; a larger absolute value reflects a relatively greater difference in the number of highly and lowly expressed elements, indicating a potentially higher degree of regulation). The representative pathways enriched in molecular function (MF), cellular component (CC), and biological processes (BP) were displayed in chord diagrams in order of absolute Z value, with an adjusted p-value < 0.05.

2.5. Identification of mRNAsi-Related DEGs and Candidate Target Drugs

The co-expressed gene modules that group similar gene expression patterns were identified via WGCNA co-expression network analysis. The following are the primary stages that comprise the analysis process: (1) obtained the optimal soft threshold power β value through the pick-Soft-Threshold algorithm; (2) the β value was selected to be used in the formation of the proximity matrix, (3) which was subsequently converted into the topological overlap matrix (TOM); and (4) a hierarchical clustering structure was created to incorporate the nearby gene modules.

Correlation analysis was employed to screen the most significant modules associated with stemness levels (mRNAsi and EREG-mRNAsi indices). Following the determination of the crucial stemness-related module, all of the genes contained in the module were utilized to search for possible target drugs for PPGLs using the signature matching methods (XSum algorithm) based on the CMap dataset (CMap Build 2: 3 cell lines, 1309 compounds) [10]. Subsequently, the XSum score ranging from −1 to +1, representing the pairwise similarity metrics between the expression patterns of cells treated with drugs and disease signature, was calculated. Briefly, a drug with XSum scores between −1 and 0 suggested its stronger reversal potency and implementation potential.

2.6. Construction of Risk Score Model Based on mRNAsi-Related DEGs

The Kaplan–Meier (K-M) survival analysis was utilized to eliminate the mRNAsi-associated prognostic genes. Furthermore, LASSO regression analysis was conducted to construct the best gene signature derived from the mRNAsi-related prognostic genes using the glmnet R package. The mRNAsi-associated prognostic risk model was created with the following formula: risk score = Expression mRNAi × Coefficient mRNAi, where i represented the sequence of mRNAsi-relevant genes.

The TCGA-PPGLs cohort’s risk scores for each specimen were estimated in line with the model, and the median risk score was chosen as the split point to classify patients into high- and low-risk subgroups. The impact of the risk characteristics on the prognosis of PPGLs in the high- and low-risk subgroups were assessed using K-M curves. Next, receiver operating characteristics (ROC) curves and principal component (PCA) evaluations were used to assess prediction accuracy. A multivariate logistic regression analysis was utilized to ascertain if the risk characteristic could independently reflect the clinical outcome in PPGLs patients.

The relative quantity of the immune cells was calculated in TCGA-PPGLs samples using CIBERSORT [11], EPIC [12], xCell [13], quantTIseq [14], MCP-counter [15], and TIME algorithms [16]. The ssGSEA algorithm was employed to estimate the tumor infiltration fraction and the activities of all steps in the anti-cancer immunity cycle for TCGA-PPGLs samples based on the TIP database (http://biocc.hrbmu.edu.cn/TIP/, accessed on 30 June 2023) [17]. The TIDE algorithm [18] was used to determine the response of the TCGA-PPGLs samples to immunotherapy, allowing the patients to be categorized as non-responders or responders. The diverse immune signatures were compared in high- and low-risk populations, including tumor-infiltrating lymphocyte infiltration, anti-cancer immune activities, several immune checkpoints expression, immune scores, stromal scores, ESTIMATE scores, and immunotherapy response.

The association between immune cells and risk score was assessed using the R package Corrplot based on the data from the immune infiltration obtained through the CIBERSORT algorithm.

It is difficult to discover more PPGLs datasets with prognostic or follow-up data for immunotherapy. The immunotherapy dataset IMvigor210 [19] was chosen as the independent external validation cohort. The patients were divided into high- vs. low-risk groups based on the median value of all risk scores. The prognostic value of the mRNAsi-associated signature was validated using the K-M method. The response for immunotherapy based on the IMvigor210 cohort included complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). Immunotherapy responses were evaluated between high- and low-risk subgroups. Meanwhile, we analyzed the differences in genomic mutation landscape and heterogeneity between high- and low-risk subgroups, including tumor mutational burden (TMB), mutant-allele tumor heterogeneity (MATH), microsatellite instability (MSI), tumor neoantigen, homologous recombination deficiency (HRD), and loss of heterozygosity (LOH). Tumor heterogeneity played important roles in the progression, invasion, and drug sensitivity of tumor cells.

2.7. Statistical Analysis

Student’s t- or Wilcoxon Rank-Sum tests were employed to compare continuous parameters between groups. Spearman’s correlation analyses were used to estimate the relevance between the two parameters of interest. A p value less than 0.05 was deemed statistically significant. All the aforementioned were achieved by various R packages implanted in the R Foundation (v3.6.3, https://www.r-project.org/ accessed on 30 June 2023).

3. Results

3.1. Correlation Between mRNAsi and PPGLs Molecular Mutation

We examined mRNAsi differences in distinct groups after classifying patients into different groups according to their molecular and clinical state. Figure 1A–C illustrates that the patients with HRAS, CSDE1, or NF1 somatic mutation have higher mRNAsi than wild-type patients. We observed that high mRNAsi was associated with RET germline mutation when combined with germline mutation data (Figure 1E). However, the stemness index was significantly lower in VHL-somatic or germline mutant samples than in VHL-wild-type samples (Figure 1D,F). Consistent with earlier studies, the TCGA expression subtypes displayed distinctive mutational profiles; germline and somatic mutations in VHL were particular to the pseudohypoxia subtype.

3.2. Correlation Between mRNAsi and PPGLs Clinical Features

Regarding the clinical features, numerous retrospective studies suggest that patients with SDHB mutations are more susceptible to experiencing metastatic disease [20,21]. So, we investigated the connection between mRNAsi and the likelihood of metastasis in patients lacking SDHB mutations. Figure 2A shows that the metastatic patients had higher mRNAsi, suggesting that the mRNAsi index may influence the PPGLs malignant process. We also found that the pseudohypoxia subtype showed lower mRNAsi scores than the kinase signaling and the Wnt-altered subtypes (Figure 2B).

Additionally, we discovered that the mRNAsi values were significantly related to clinical outcomes. Figure 2C–F demonstrate that the dead patients have higher mRNAsi scores than alive ones, while patients with different ages, gender, and tumor location of PPGLs did not indicate significant differences in mRNAsi levels. This indicates that miRNAsi may only influence the prognosis of PPGLs, but this cannot be considered as clinically significant evidence due to limited data on death cases. Therefore, the survival analysis was performed, and the results further indicated the K-M curve (Figure 2G–I) shows significant difference in overall survival (OS), disease specific survival (DSS), and disease free interval (DFI) between the two groups with high and low mRNAsi expression in the cases we studied. The high-mRNAsi cluster tended to show a relatively poorer survival outcome, with the shorter OS/DSS or shorter disease-free survival period than low-mRNAsi subtypes(OS, HR = 6.52, 95% CI 1.30–32.66; DSS, HR = 11.55, 95% CI = 1.33–100.16).

3.3. Correlation Between mRNAsi and PPGLs Immune Signatures

The ssGSEA analysis was employed to evaluate the TME of the TCGA-PPGLs cohort. Unsupervised clustering was utilized to divide TCGA-PPGLs patients into immunity-H and immunity-L groups (Figure 3A). The high-immunity subgroup was characterized by relatively higher average immune, stromal, ESTIMATE, and mRNAsi scores than the low-immunity subgroup (Figure 3B–F). Meanwhile, we found that the mRNAsi score had a negative association with immune, stromal, and ESTIMATE scores (Figure 3G–I), suggesting that immunosuppression in PPGLs has a significant correlation with the higher degree of stemness.

Additionally, we compared the TME infiltrating fractions between the two subgroups using CIBERSORT analysis. Figure 4A illustrates that immunity-H has higher anti-tumor TME components like activated memory CD4+ T cells, macrophages M1 cells, and dendritic cells. The correlation analysis revealed a relationship between weak CD4+ T and mast cells and mRNAsi abundance (Figure 4C). We further observed elevated expressions of several immune checkpoints in the high immune infiltration group, such as CD86, CD274, CD80, CD276, CTLA4, PDCD1, and PDCD1LG2. However, this trend was reversed in the high mRNAsi group (Figure 4B–D). These findings imply that the tumor stemness index is inversely correlated with TME infiltration. The high mRNAsi group may be linked to poor immune checkpoint treatment efficacy.

3.4. Identification of Differentially Expressed Genes and Functional Enrichment Analyses

A total of 397 DEGs were identified between the high- and the low-mRNAsi groups (Figure 5A,B). We performed the GO/KEGG analysis to investigate the latent biological behavior and the related pathways involved in these DEGs. Figure 5C displays that 397 DEGs mainly enriched several classic tumor signaling pathways, including ERK1 and ERK2 cascade regulation and extracellular structure organization. Protein digestion and absorption, inorganic anion transport, and mononuclear cell migration are active pathways during tumorigenesis and pathways related to immune response regulation. This could explain why the high mRNAsi group had lower immune activity and a higher malignant potential than the low mRNAsi group.

3.5. Construction of the mRNAsi-Related Prognostic Model for PPGLs

As previously described, we employed WGCNA to identify the key gene modules of highly correlated DEGs that correlated to mRNAsi and EREG-mRNAsi traits. When soft threshold β = 4, R^2^ > 0.9, this allowed the construction of a scale-free co-expression network (Figure 6A,B). Moreover, DEGs with comparable expression patterns were merged into the same module (Figure 6C). Ultimately, three cohesive gene modules were yielded (Figure 6D). The blue module was the module most relevant to both mRNAsi indices (cor = −0.79; p < 0.001) and EREG-mRNAsi indices (cor = 0.54; p < 0.001). Then, 141 mRNAsi-related DEGs from the blue module were employed for querying CMap. Five agents (X4.5. dianilinophthalimide, TTBPB, imatinib, X5224221, calmidazolium) were then considered to have anti-stemness potency (XSum score < 0) (Figure 6E).

The K-M survival analysis was utilized to estimate the prognostic value of 141 mRNAsi-related DEGs in PPGLs. Thirty-three genes connected to PPGLs patient’s OS were identified using p < 0.05 as the cutoff (Supplementary Table S1).

The LASSO algorithm was applied to the 33 prognostic-related genes, and six genes were identified to develop the mRNAsi-risk signature (Figure 7A,B). The six genes were CCL26, GPR148, FAM3B, GATA4, C1orf94, and CRYBB3. Each PPGLs patient risk score was determined as follows:

CCL26 expression × (0.351137846036265) + GPR148 expression × (0.601409843037207) + FAM3B expression × (−0.746029637361062) + GATA4 expression × (−1.04632515461906) + C1orf94 expression × (1.0041025435075) + CRYBB3 expression × (−0.446009476652402).

The risk score for each TCGA PPGLs patient can be generated by applying this risk model. We observed that the mRNAsi score was positively associated with the risk score (Figure 7F). All TCGA-PPGLs patients were divided into two clusters using the median risk score. Figure 7C displays that the group with high levels of risk score had a greater number of fatalities than the low mRNAsi risk group. The PCA analysis (Figure 7D) showed that the high- and low-risk patients had clinical outcomes about OS. The model’s projected accuracy was assessed using the site under the ROC curve, and it showed outstanding prognostic efficacy (AUC = 0.9082) (Figure 7E). Additionally, univariate Cox analysis demonstrated that clinical aggressive status and risk scores were closely related to PPGLs prognosis (Figure 7G). The prognostic performance of this cancer stemness-based prediction model suggested that it can achieve satisfactory predictive accuracy in survival analysis for PPGLs patients.

3.6. The mRNAsi Risk Group Displayed Distinctive TIME and Immunogenomic Profiles

We performed the CIBERSORT, EPIC, xCell, quantTIseq, MCP-counter, and TIME algorithms to assess the relative abundance of immune-associated cell types based on the expression patterns of PPGLs patients (Figure S1A–G). We detected that at least five algorithms elucidate the macrophages, especially the M1 macrophages, were more abundant in low-risk patients (Figure S1B–G). In addition, low mRNAsi risk was associated with high immune checkpoint expressions, such as PDCD1, CD276, CD80, PDCD1LG2, and CD86 (Figure 8A). Meanwhile, we observed more immune-related functions like tumor killing, NK cells, and neutrophils in low-risk patients. Patients with low tumor risk exhibited higher anti-tumor immunity than high tumor risk (Figure 8B). Correspondingly, the immune, stromal, and ESTIMATE scores were negatively related to the risk score (Figure 8C–E). Furthermore, the analysis revealed that PPGLs patients with a low mRNAsi risk were more responsive to immunotherapy (Figure 8F). To further verify the prognostic value of the mRNAsi risk model for clinical outcomes and immunotherapy sensitivity, an independent external validation cohort verified the risk score model based on six mRNAsi-related DEGs. In the IMvigor210 cohort, we used the same method to assign 348 patients to low-risk groups (N = 174) or high-risk groups (N = 174). As expected, patients in the low-risk group had a higher CR/PR rate (Figure 8G).

Meanwhile, we also could observe that CCL26, GPR148, FAM3B, GATA4, and CRYBB3 expression were significantly associated with the immune and stromal scores (Figure 9A). Genomic heterogeneity characteristics, for instance, TMB, MSI, and tumor neoantigen, were important predictors for tumor immunotherapy response within the tumor microenvironment. In our study, although there was no significant difference in the genomic mutation landscape between the high- and low-risk subtypes, we noticed that the low-risk PPGLs patients showed higher tumor neoantigens (Figure 9B,C). Therefore, the favorable prognosis of low-risk patients was found to correlate with the immune microenvironment (immune cell infiltration), high immune checkpoint expressions and tumor neoantigen.

3.7. The mRNAsi Risk Group Displayed Distinctive Therapeutic Drug Sensitivities

Temozolomide (TMZ) was the recommended chemotherapy drug for PPGLs postoperatively, however, we found that there was no significant difference in TMZ response between the two groups. Subsequently, we determine whether distinct mRNAsi risk groups exhibited different susceptibility to other therapeutic drugs. Patients with a high mRNAsi risk had higher estimated IC50 values for Erlotinib, I-BET-762, AZD5153, and WEHI-539 (Figure 9D).

4. Discussion

Numerous types of research demonstrated that CSCs were related to the self-renewal and long-term tumor maintenance of the subpopulation of tumors [22,23]. Additionally, it has been demonstrated that CSCs are crucial for carcinogenesis, metastasis, tumor self-replication capacity, tumor metabolic reprogramming, and the remodeling of the tumor immunological microenvironment [24]. Factors highly correlated with tumor stemness were identified using a comprehensive bioinformatics analysis to obtain insights into the stemness features of PPGLs. A stemness scoring system was developed to predict patient prognosis and immunotherapy responsiveness.

The mRNAi is an index that describes the similarity between tumor and stem cells to explore their association with phenotypic characteristics. Now, PPGLs cannot be distinguished between benign and malignant. Instead, all PPGLs are considered potentially malignant [25]. Malignant metastasis occurs in 15% of pheochromocytomas and half of abdominal paragangliomas. Many studies also explore histological and molecular predictive markers to predict tumor metastasis [26]. Although current reports show that the probability of metastases in patients without SDHB mutations is low, some patients still develop metastasis [20,21]. Our study explored mRNAsi in patients without SDHB mutations and found that PPGLs patients with metastasis had higher mRNAsi scores. This suggests that the degree of tumor stemness of PPGLs affects the patient metastasis and progression, and it may predict the metastasis probability in patients without SDHB and cancer progression through changes in its stemness characteristics in the future.

The TME plays an important role in tumor cell renewal, proliferation, and differentiation, especially immune-related signals in the microenvironment that can affect patients’ prognosis and survival rate. According to previous studies, CSCs have a regulatory effect on tumor immunity, which is consistent with our findings in PCPG [27]. Our study revealed that mRNAsi scores were negatively correlated with immune, stromal, and ESTIMATE scores, suggesting that the stronger the stemness characteristics of PPGLs, the weaker their immune function, and the tumor may escape immune clearance under the condition of low immune function. Patients with a low mRNAsi risk exhibited high immune activation characteristics, and ultimately with a good clinical prognosis. As shown in Figure 4A, our study analyzed the infiltration characteristics of immune cells, between immune-activated (Immunity_H) and immune-suppressive (Immunity_L) subtypes, suggesting the TME landscapes in low mRNAsi and high mRNAsi risk patients. Notably, various immune cells that exert anti-tumor roles, such as M1-like macrophages, NK cells, and neutrophils, were enriched in patients with low mRNAsi risk. In addition, the infiltration of M2 macrophages with pro-tumor roles also significantly decreased than the high mRNAsi subtypes. Previous studies have reported that M1-like macrophages have an antitumorigenic role when the innate and acquired immune systems are activated. The dramatic infiltration of M1-like macrophages impaired the growth, migration, and angiogenesis of various solid tumors [28], which was related to a better survival rate. According to several studies, tumor-associated macrophages, particularly M2-like macrophages, directly communicate with CSCs to support their ability to migrate, self-renew, and survive. However, the effects of classical (M1)-polarized macrophages on CSCs are unclear. Luo et al. reported that unlike the alternative (M2)-polarized macrophages, M1-like macrophages inhibit the self-renewal and expansion of the hematopoietic stem cell (HSC) by inducing nitric oxide (NO) production from catalytic intracellular L-arginine, because the exposure to NO further induces HSC and hematopoietic progenitor cell apoptosis [29]. The above result also explains why PPGLs patients with low tumor stemness risk have a high proportion of M1-like macrophage infiltration. Many studies indicated that TAMs in the TME often exhibited M2 macrophages, characterized by the production and secretion of immunosuppressive cytokines. M2-like TAMs secreted cytokines, inhibited the function of cytotoxic T lymphocytes, blocked Th1-type immune activity, promoted anti-immunity Th2-type activity, thus allowing tumors to evade immune surveillance and thereby promoting tumor cell growth, drug resistance, and angiogenesis [30,31,32]. Therefore, compared with the high mRNAsi, the lower M2 macrophages in the low mRNAsi subtypes, also caused better clinical prognosis. Meanwhile, several immune checkpoint genes were suppressed in the high mRNAsi risk population, attenuating the immune response. Thus, we proposed that the CSC marker expression in tumor tissues of PPGLs has a close association with the immunosuppressed immune microenvironment and revealed poor survival outcomes.

Recently, immunotherapy has been widely investigated and applied to various tumors. Patient-tailored immunotherapy has provided tremendous advantages to patients with neural crest-derived cancer, especially for melanoma and neuroblastoma. However, the clinical scientific experience of these phase II clinical trials with checkpoint inhibitors and other types of immunotherapy for PPGLs was not satisfactory [33,34]. In the current study, PPGLs patients in the low-mRNAsi risk group had higher antitumor immune activity and immune checkpoint expression. We had a hunch that the mRNAsi risk value could reveal how much patients would benefit from immunotherapy. As expected, immunotherapy had a substantially better effect on patients classified as a low-risk group. IMvigor210, an authentic immune therapy cohort that investigated the efficacy of PD-L1 inhibitor therapy on patients with metastatic endothelial cancer, was used for validation. Our study demonstrated that the risk score model might be appropriate for predicting the efficacy of immunotherapy in PPGLs and beneficial for patients with other types of cancer. Medical staff can choose a different treatment strategy depending on the two different mRNAsi signature subtypes above. We believed that PPGLs patients with low miRNAsi risk, characterized by the immune-activated TME, could benefit more from immunotherapy, and in addition, their estimated IC50 for drugs such as erlotinib and I-BET-762 is lower. Some studies also indicated that the immunosuppressive tumor microenvironment (TME) has lower sensitivity to immunotherapy, for example, in prostate cancer [35]. Zhou et al. demonstrated the low-score stemness subgroups were the significant responders to anti-CTLA4/PD-1 treatment in clear cell renal cell carcinoma [36]. These findings were consistent with our results. Therefore, the response of tumor to treatment was also closely related to the interaction of TME. Many studies have found that neoadjuvant chemotherapy combined with immunotherapy could increase CD8+ T cell infiltration more significantly than neoadjuvant chemotherapy alone, prolonging overall patient survival [37]. Tumors with better treatment effects (higher pathological complete response) exhibited a high M1 polarization level before treatment, and appear to have a higher M1/M2 ratio after treatment [38,39].

Targeting tumor stem cells is the top priority in treating and preventing tumor recurrence. Bruno et al. proposed that disrupting CSC-TAM crosstalk or combining CSCs and TAMs targeting may be a new therapeutic strategy [40]. In the future, we can thoroughly grasp the mechanism of tumor stem cells and completely eliminate tumors at the stem level. Although we performed a prognostic analysis of our risk model in the TCGA dataset, we did not find other external sequencing datasets to validate it due to the rarity of PPGLs. We may discover relevant and complete follow-up data to supplement the experimental validation. Moreover, the mechanism of interaction between CSCs and various immune components needs further experimental exploration in the TME of PPGLs.

5. Conclusions

In conclusion, we explored the relationship between tumor stemness characteristics and the clinical properties and immune microenvironment of PPGLs in TCGA to develop a stemness scoring system that can be applied to predict prognosis and patient immunotherapy responsiveness. The occurrence, metastasis, and prognosis of malignancies are strongly correlated with the stemness characteristics of tumors. Therefore, evaluating tumor stemness may help determine the extent of the tumor and prognostic indicators and serve as a guide for future pharmacological target therapy.