Cell signaling communication within papillary craniopharyngioma revealed by an integrated analysis of single-cell RNA sequencing and bulk RNA sequencing

Abstract

Objective

This study aims to elucidate the primary signaling communication among papillary craniopharyngioma (PCP) tumor cells.

Methods

Five samples of PCP were utilized for single-cell RNA sequencing. The most relevant ligand and receptor interactions among different cells were calculated using the CellChat package in R software. Bulk RNA sequencing of 11 tumor samples and five normal controls was used to investigate the pair interactions detected by single-cell RNA sequencing.

Results

Fibroblasts were not found in ACP, whereas they were detected in PCP. InferCNV revealed high CNV scores for the clusters of epithelial cells and fibroblasts using immune cells as a reference. Epithelial Mesenchymal Transition, Interferon Gamma Response, p53 Pathway, and Estrogen Response Early are pathways commonly shared by fibroblasts and epithelial cells, ranking high in priority. The Wnt signaling pathway and PI3K-Akt signaling pathway play a crucial role in facilitating communication between epithelial cells and fibroblasts. Neutrophils were recognized as the main receivers of incoming signals, with ANXA1-FPR1 and MIF-(CD74 + CXCR2) being identified as the primary signals transmitted from fibroblasts to neutrophils.

Conclusion

Through analyzing the communication of essential signaling pathways, ligands, and receptors among epithelial cells, fibroblasts, and neutrophils in PCP tumor tissues, we have identified certain molecules with promising prognostic and therapeutic potential.

Introduction

Craniopharyngiomas (CPs) primarily originate from the epithelial cells of the craniopharyngeal duct, making them a complex type of brain tumor due to their deep location and close proximity to key nerve and vascular structures [1]. CPs are typically located in the sella turcica or above it, and while they can develop anywhere along the pituitary-hypothalamic axis, approximately 50% of cases originate at the base of the third ventricle, extending into its cavity [2, 3]. According to the World Health Organization Classification of Tumors of the Central Nervous System, CPs are classified into Adamantinomatous craniopharyngioma (ACP) and Papillary craniopharyngioma (PCP) [4], with ACP being more common in children and PCP having a higher incidence rate in adults [5].

PCP is a rare histological subtype with squamous-papillary characteristics, typically adhering strongly to the tumor and hypothalamus [6], making surgical resection likely to damage the hypothalamus and lead to psychiatric disorders [7]. Although PCP is a benign tumor, complete surgical resection is not always feasible due to its proximity to the nervous system, pituitary gland, and third ventricle [8]. Tumors that are completely or partially removed through surgery also have a high recurrence rate. Therefore, finding better drug targets is crucial for the treatment of PCP [9].

Single-cell RNA sequencing (scRNA-seq) involves collecting disease samples and/or normal samples, followed by sequencing cell populations in tissues using tools like 10 × Genomics to quantitatively measure molecular changes at a single-cell resolution [10, 11]. This approach provides a more sensitive understanding of each cell's role in tissue function and disease development [12]. In tumor research, scRNA-seq can help identify cell interactions and transformations [13], aiding in the exploration of physiological and pathological tumor mechanisms [14]. Since its inception, scRNA-seq has proven to be a powerful tool for studying disease onset and progression, offering new possibilities for clinical treatment approaches by analyzing cell-to-cell interactions in more detail. While there is a substantial amount of scRNA-seq research on tumor tissues or normal tissues, to our knowledge, there is a lack of scRNA-seq research on PCP. In this study, we performed scRNA-seq on samples from five PCP patients and bulk RNA sequencing on 11 tumor samples and five normal controls. Our aim was to gain insight into the tissue characteristics of PCP at the single-cell level and identify a series of key genes, ultimately providing data reference for the precise treatment of PCP.

Materials and methods

Patient samples

From 2019 to 2022, five patients diagnosed with PCP at the Beijing Tiantan Hospital were selected for scRNA-seq. The tumor tissue extracted from the patient was sequenced using the 10X Genomics platform, with Shanghai OE Biotech Co., Ltd. conducting the sequencing and initial sequencing data analysis. This study was approved by the ethics committee of Beijing Tiantan Hospital. Informed consents were received from the parents prior to the study. Detailed information of the patients is included in this study was included in Supplementary Table 1. Single-cell sequencing data, along with bulk RNA sequencing data (GSE243682), pertaining to various neuropathologies including High-grade glioma (GBM), Pituitary neuroendocrine tumors (PitNET), Lewy Body Dementia (LBD), Arteriovenous malformations (AVM), Low-grade glioma (LGG), Lung brain metastasis (MET), and Meningioma (MEN), were retrieved from the Gene Expression Omnibus (GEO) database for confirmatory analysis. The specific data sets for single-cell sequencing include the following GEO accession numbers: GSE242044, GSE256493, GSE235914, GSE278450, and GSE244101. Access to the database can be facilitated through this link: https://www.ncbi.nlm.nih.gov/gds/?term = .

scRNA-seq data processing

Initially, the samples underwent quality control using the official software Cell Ranger from 10 × Genomics, which internally integrates the STAR software. After aligning the reads to the reference genome, quality control results such as the number of high-quality cells, genes, and genome mapping rates were obtained from the raw data to evaluate the quality of each sample. Following the initial quality control by Cell Ranger, further quality control was conducted on the cells: each cell had to have a gene expression count of at least 200, UMI count greater than 1000, mitochondrial UMI ratio less than 20%, and red blood cell gene ratio less than 5% to be considered high-quality cells.

Processing of scRNA-seq data involved subjecting the cells obtained from the previous step to downstream analysis using the CCA method from the R package "Seurat" to integrate the cells, and remove batch effects. This process resulted in 16 cell clusters through dimension reduction and clustering. Twelve cell types were identified using the SingleR function: NK_cell, T_cell, B_cell, Endothelial, Epithelial, Fibroblast, Macrophage, Monocyte, DC, Mast_cell, Neutrophil, and Plasma_cell. Further analysis was conducted on Fibroblast and Epithelial cells using PCA and UMAP dimension reduction and clustering. Subsequently, the FindAllMarkers function was utilized to identify marker genes for each cluster.

Gene set functional analysis

The gene sets MSigDB_Hallmark_2020 and WikiPathway_2021_Human were downloaded through the Enrichr website (https://maayanlab.cloud/Enrichr/). Gene set enrichment analysis on marker genes was performed using the R packages of "enrichr" and "clusterProfiler". The "AUCell" package in R was used to determine the activity levels of enriched pathways in individual cells. Initially, the AUCell_buildRankings function was used to rank all genes within each cell. Subsequently, the AUCell_calcAUC function was applied to assign scores to the ranked gene sets. Finally, a UMAP graph was generated based on these scores to visually represent the activity levels of pathways in each cell.

InferCNV analysis

The copy number variations (CNVs) of epithelial cells and fibroblasts were examined in R using the "infercnv" package. The analysis followed the official instructions to obtain the raw count matrix, cell type annotation file, and gene/chromosome position annotation file. NK cells were selected as the standard normal cells, with a denoising value set at 2.5, while all other parameters were kept at their default values.

Cell–cell communication analysis

We referenced the CellChatDB database to conduct an analysis of cell–cell interactions in secreted signaling, utilizing the “CellChat” R package. We randomly selected a subset of 20,000 cells to analyze their interactions. The cell types were categorized by examining both incoming and outgoing signals with the aid of the “NMF” package. Based on the Cophenetic and Silhouette indices, we classified the signals into four distinct patterns for both incoming and outgoing communications. We identified ligand-receptor gene pairs that were closely associated with each of these patterns.

Initially, we used the “computeCommunProb” function to ascertain the probabilities of cell communication, thereby inferring the communication networks. Subsequently, we applied the “filterCommunication” function to eliminate signaling networks that did not meet the minimum threshold of min.cell = 10 cells in the pathway. After this filtering, we employed the “computeCommunProbPathway” function on the remaining signaling pathways to gauge the communication intensity between cells, with a significance threshold set at thresh = 0.05 to discount minor cell communication events. We then compiled an overview of the entire communication network across all cells using the “aggregateNet” function. To evaluate the potential for overfitting or biased results from the “CellChat” analysis, we also scrutinized cell–cell communication networks within specific cell subtypes. In this study, we carried out “CellChat” analyses on all cell types, as well as on interactions where epithelial cells acted as ligands and fibroblasts as receptors, and vice versa, with fibroblasts as ligands and neutrophils as receptors.

Pseudotime trajectory analysis

Cell subtype developmental and evolutionary paths were inferred through pseudotime analysis conducted with Monocle (version 2.26.0). The software focused on genes with greater variability in expression levels across different cells. Utilizing the DDRTree algorithm within the reduce Dimension function, a minimum spanning tree was then built to outline cell differentiation paths. This method guarantees connectivity among cells while reducing the cumulative edge weight. Finally, the derived trajectory was depicted using the plot_cell_trajectory function.

Single-cell regulatory network inference and clustering (SCENIC)

To delve deeper into the mechanisms underlying the heterogeneity among various cell clusters and to examine the expression patterns of differentially expressed genes, we utilized single-cell regulatory network inference and clustering to construct a transcription factor regulatory network (SCENIC, version 1.1.0.1) [24]. In essence, SCENIC identifies co-expression modules between transcription factors (TFs) and potential target genes, as well as calculates regulatory activity scores (RAS) for individual cells. The regulatory specificity score (RSS) is calculated to forecast the precise relationship between regulatory factors and specific cell types. The connectivity specificity index (CSI) is employed to indicate the interplay between different regulatory rules. Regulatory genes with a higher CSI are likely to co-regulate downstream genes and collectively contribute to cellular functions.

Analysis of PCP bulk RNA sequencing data

Detailed information on bulk RNA sequencing of 11 tumor samples and five normal controls was previously reported [15]. Screening criteria of Fold Change > 2 and P-value < 0.05 were applied, and heat maps of differential genes were created based on the findings. Subsequently, pathways obtained from CellChat were used to perform correlation analysis of ligand-receptor pairs. The screening criteria for correlation were set at Cor > 0.5 and P-value < 0.05. Potential therapeutic targets were identified by analyzing the results of this correlation analysis.

Results

Cell subtypes in PCP tissues revealed by scRNA-seq

To investigate the internal composition of PCP tumor tissue, scRNA-seq was utilized to sequence samples from five PCP patients. Following standard procedures such as quality control and removal of empty cells, 29,058 high-quality cells were successfully obtained. Through standardization, normalization, and dimensionality reduction clustering, 12 distinct cell types were identified (Supplementary Fig. 1A). The widely recognized crucial role of epithelial-mesenchymal transition (EMT) in tumor development led us to focus primarily on analyzing the status of epithelial cells, resulting in a total of 8957 cells being collected. The top 1000 highly variable genes were selected for the MNN batch correction algorithm and dimensionality reduction clustering. Subsequently, four distinct cell clusters were identified (Supplementary Fig. 1B). These clusters were characterized by marker genes CTNNAL1 and PDPN, MME and CYSRT1, LRRC32 and GZMB, and TPSB2 and CD69, respectively. Violin plots were used to visually represent the marker genes for each cluster, followed by a UMAP visualization illustrating the expression of these markers in the cells (Supplementary Fig. 1C, D).

Growing evidence suggests that communication and interactions between cancer cells and cancer-associated fibroblasts (CAFs) play a crucial role in tumor metastasis and progression [16]. Given that PCP and ACP are both subtypes of craniopharyngioma, our initial analysis focused on the fibroblasts reported in ACP. The results revealed that fibroblasts were not identified in ACP (Supplementary Fig. 1E), consistent with previous findings [17]. Our study confirmed the presence of fibroblasts in PCP tumor tissues, indicating distinctions between the two subtypes of craniopharyngioma and suggesting a potential role of fibroblasts in the pathogenesis of PCP tumors. A total of 518 fibroblasts were identified in the five PCP samples. The top 1000 highly variable genes were chosen for principal component analysis (PCA). Using the FindNeighbors and FindClusters functions, all cells were assigned to different clusters. Following dimensionality reduction clustering, five distinct cell clusters were identified (Fig. 1A). Heatmaps generated with the top 10 marker genes were displayed in Fig. 1B. A violin plot was created using marker genes for each of the five clusters: RERGL and MYOCD for the first cluster; APOD and SFRP4 for the second cluster; MTHFD2 and MX2 for the third cluster; VWF and RAMP2 for the fourth cluster; and PTGER1 and REM1 for the fifth cluster. As depicted in Fig. 1C, the majority of the marker genes effectively differentiated between the cell clusters.

Fig. 1

Cell types of fibroblasts revealed by scRNA-seq. A UMAP plots of fibroblasts in PCP tissues based on scRNA-seq data. B Heatmap showing the top 10 marker gene expression signatures. C Violin plots displaying the expression levels of distinct marker genes for five primary fibroblast subclusters

InferCNV analysis on epithelial cells and fibroblasts

To investigate the CNVs, InferCNV was conducted using NK cells as reference normal cells. The graph revealed evident chromosome deletions and amplifications in epithelial cells (Fig. 2A). Notably, the CNV score of the fourth cluster of epithelial cells was the highest, suggesting a higher level of malignancy within that cluster (Fig. 2B). The five clusters of fibroblasts also exhibited chromosomal deletions and amplifications, with cluster 4 displaying the highest CNV score (Fig. 2CD). In summary, the presence of a significant number of CNVs suggests that fibroblasts may contribute to the development and progression of PCP tumors.

Fig. 2

InferCNV analysis results of epithelial cells and fibroblasts. A The heatmap displays the chromosomal mapping of epithelial cells’ single-cell large-scale copy number variations inferred through scRNA-seq. B CNV score of epithelial cells. C The heatmap shows the chromosomal mapping of single-cell large-scale CNVs in fibroblasts inferred through scRNA-seq. D CNV score of fibroblasts

Functional analysis of epithelial cells and fibroblasts

To elucidate the functional roles of epithelial and fibroblasts, unique marker genes for each cell type were selected, and the corresponding pathways were enriched using data from the HALLMARK and WikiPathways databases. The analysis revealed that epithelial cells were predominantly enriched in pathways such as Cytoplasmic Ribosomal Proteins, Allograft Rejection, Vitamin D Receptor Pathway, Epithelial Mesenchymal Transition, TNF-alpha Signaling via NF-kB, and Estrogen Response Late (Fig. 3A and Supplementary Table 2). Fibroblasts were primarily enriched in pathways including Focal Adhesion-PI3K-Akt-mTOR-signaling pathway, PI3K-Akt signaling pathway, Focal Adhesion, Epithelial Mesenchymal Transition, Apical Junction, and Interferon Gamma Response (Fig. 3B and Supplementary Table 3). The top10 pathways of each datasets shared by both epithelial cells and fibroblasts included the Epithelial Mesenchymal Transition, Interferon Gamma Response, p53 Pathway, and Estrogen Response Early.

Fig. 3

Enrichment analysis result of epithelial cells and fibroblasts. Functional enrichment of epithelial cells (A) and fibroblasts (B) using HALLMARK and Wiki Pathways databases. C AUCell scoring of the activity of the p53 pathway in different cell clusters of epithelial cells showed significantly lower activity in cell cluster 4. D The activity of the p53 pathway was relatively low in the fibroblast cell clusters

The protein p53 plays a crucial role in regulating the cell cycle, DNA replication, and controlling cell division to prevent excessive cell proliferation [18]. Subsequently, AUCell scoring was conducted on the pathway to assess the activity of the p53 pathway in different cell clusters of epithelial cells and fibroblasts. In epithelial cells, the activity of the p53 pathway was notably lower in cell cluster 4 compared with the other clusters (Fig. 3C). This suggests that the p53 protein may be less effective in regulating cell cluster 4 of epithelial cells, potentially contributing to increased malignancy. This observation aligns with the findings from the InferCNV analysis. In fibroblast cell clusters, the activity of the p53 pathway was relatively low (Fig. 3D). Mutations in normal epithelial cells are primarily responsible for tumor formation, while fibroblasts may interact with other cell types in the TME to influence tumor development and invasion.

Signaling communication between epithelial cells and fibroblasts through outgoing patterns

To explore cellular interactions in PCP samples and identify highly correlated ligand-receptor pairs between them, an analysis of output signals was conducted using CellChat. The findings indicated that epithelial cells significantly contributed to output signals across various cell subtypes (Fig. 4A). Subsequently, the “NMF” package was utilized to classify cell types and their associated signaling pathways (Fig. 4B). As depicted in Fig. 4C, cell patterns were classified into four categories based on Cophenetic and Silhouette metrics: Pattern 1 (Epithelial_1, Epithelial_2, and Epithelial_3), Pattern 2 (DC, Monocyte, Macrophage, and Neutrophil), Pattern 3 (Fibroblast_2, Fibroblast_5, Fibroblast_1, and Fibroblast_3), and Pattern 4 (NK_cell and T_cell). Communication patterns were also divided into four types, and the ones with the highest contributions were selected for subsequent analysis. The primary contributors for Pattern 1 included GDF, CSF3, ncWNT, BAG, SAA, PTN, MK, WNT, CX3C, PERIOSTIN, and SEMA3 (Fig. 4C and Supplementary Table 4, 5).

Fig. 4

Signaling communication between epithelial cells and fibroblasts through outgoing patterns. A Epithelial cells played a significant role in producing output signals among different cell subtypes. B, C Based on Cophenetic and Silhouette metrics, cell patterns were categorized into four groups using Nonnegative Matrix Factorization (NMF) method. D Communication interactions between epithelial cells and fibroblasts uncovered by analyzing Pattern 1 within the communication patterns. E Enrichment analysis result using the fibroblast subpopulation of Pattern 3 as the receptor and the epithelial cell subpopulation of Pattern 1 as the ligand

To investigate the communication interactions between epithelial cells and fibroblasts, an enrichment analysis was conducted using Pattern 1 from the Communication patterns. Initially, the analysis focused on epithelial cells as ligands and other cells as receptors. The pathways that showed significant enrichment included Wnt signaling in kidney disease, Osteoblast differentiation, Neovascularisation processes, ESC Pluripotency Pathways, and Epithelial to mesenchymal transition in colorectal cancer (Fig. 4D and Supplementary Table 6). Subsequently, an enrichment analysis was carried out using the fibroblast subpopulation of Pattern 3 as the receptor and the epithelial cell subpopulation of Pattern 1 as the ligand. This analysis highlighted enrichment in Wnt signaling in kidney disease, PDGF Pathway, ESC Pluripotency Pathways, and Cytokines and Inflammatory Response pathways (Fig. 4E and Supplementary Table 7). Based on the outcomes of the enrichment analyses, it can be inferred that the Wnt signaling pathway in kidney disease, ESC Pluripotency Pathways, and Epithelial to Mesenchymal Transition in colorectal cancer pathways play crucial roles in mediating communication between epithelial cells and fibroblasts.

Important ligand-receptor pairs between epithelial cells and fibroblasts

Through the aforementioned analysis, we have identified the significant signaling pathway between cells. Given the pivotal role of epithelial cells in the output signals of PCP tumor cells, we were intrigued to explore the signaling communication between epithelial cells and fibroblasts. In Fig. 5A and Supplementary Table 8, the ligand-receptor pairs involved in the signaling pathway were depicted along with the significance of these interactions. The most critical ligand-receptor pairs in Epithelial_1, Epithelial_2, and Epithelial_3 were identified as WNT5A-MCAM/FZD4. Additionally, the most crucial ligand-receptor pair in Epithelial_4 was found to be VEGFA-VEGFR2/VEGFR1R2. It is noteworthy that the interactions of cluster 4 epithelial cells exhibited a significant decrease compared with the other three clusters. Then, the ncWnt signaling pathway was selected for further analysis, with epithelial cells acting as ligands (Fig. 5B). Interactions among Epithelial_1, Epithelial_2, Epithelial_3, and all subgroups of fibroblasts, as well as interactions among different clusters of epithelial cells, are observable. Based on the ranking of ligand-receptor contributions, the WNT5A-MCAM pair emerged as the most significant regulatory factor in the pathway (Fig. 5C). It is worth noting that, apart from not being found in the fibroblasts cluster 2, the WNT5A-MCAM pair is widely present in other clusters of fibroblasts (Fig. 5D). To assess the potential for overfitting and biases due to the computational methods employed, we isolated the epithelial and fibroblast subpopulations for “CellChat” analysis. The findings revealed that the WNT5A-MACM ligand-receptor pair remained the most significant contributor across the analyses (Fig. S2). Since no comprehensive studies on PCP have been conducted based on single-cell sequencing technology, to verify our analysis results, single-cell sequencing data for seven neurological diseases was downloaded from the GEO database and corresponding analyses were conducted. Initially, using the “Seurat” algorithm, we found that the fibroblasts subpopulation exists only in GBM, LGG, AVM, MET, and MEN samples. Subsequently, through “Cell Chat” analysis, we discovered that in GBM samples, when fibroblasts act as receptor cells, the ligand-receptor pair WNT5A-MCAM ranks first in contribution (Fig S3A, B). Therefore, only GBM samples will be focused upon for subsequent validation and analysis.

Fig. 5

Important ligand-receptor pairs between epithelial cells and fibroblasts. A Ligand-receptor pairs between epithelial cells and fibroblasts that were involved in the signaling pathway. B The Wnt signaling pathway involves epithelial cells acting as ligands. C WNT5A-MCAM pair emerged as the most significant regulatory factor in the pathway. D The detailed communication pattern of the WNT5A-MCAM pair in subclusters of epithelial cells and fibroblasts

Next, to explore the evolutionary trajectories of the various fibroblast subpopulations in PCP disease and the changes in the ncWNT pathway as the disease progresses, we commenced our investigation with pseudotime analysis on each fibroblast subpopulation using “monocle”. This analysis uncovered three distinct differentiation states among the subpopulations (Fig. 6A), starting with Fibroblast_2 and transitioning through Fibroblast_1, Fibroblast_3, Fibroblast_5, before ultimately differentiating into Fibroblast_4. A surprising discovery was that Fibroblast_1 and Fibroblast_3 exhibited unique differentiation states during the mid-phase of the differentiation process. This led us to hypothesize that the differentiation states of each subpopulation might be influenced by the genes with high variability within those subpopulations.

Fig. 6

Pseudotime analysis of fibroblast subpopulations and highly variable genes. A Pseudotime ridge plot of fibroblast subpopulations.Pseudotime plot of highly variable genes in the Fibroblast_1 (B), Fibroblast_2 (C), Fibroblast_3 (D), Fibroblast_4 (E), and Fibroblast_5 (G) subpopulation. F Pseudotime plot of highly variable genes in the Pseudotime analysis of the WNT pathway in fibroblasts

Subsequent analysis of highly variable genes within the fibroblast subpopulations revealed that genes such as MCAM, WNT4A, ITGA5, and VMF in Fibroblast_1, TAGLN, TPM2, TPM1, and CPE in Fibroblast_3, as well as CCL3, CCR1, CXCL2, and RGS5 in Fibroblast_5 were upregulated during the mid-phase of cell differentiation (Fig. 6B–F). Conversely, genes including WNT5A, FZD1, DPT, and C1QTNF3 in Fibroblast_2 were active in the early stage of cell differentiation, whereas ANXA1, ANXA3, FPR1, and CCL14 in Fibroblast_4 showed upregulation at the terminal stage of cell differentiation. Pseudotime analysis indicated a significant increase in the expression of WNT5A, WNT4, and MCAM genes in the ncWNT pathway within the Fibroblast_1 and Fibroblast_2 subpopulations. However, pseudotime analysis of the ncWNT pathway indicated that it is primarily active during the mid to late stages of cellular differentiation (Fig. 6G). In line with the previously discussed analysis, we hypothesize that the Fibroblast_1 subtype likely assumes a significant role in the ncWNT pathway, and MCAM may represent a potential therapeutic target. In GBM samples, through pseudotime analysis of the fibroblasts subpopulation, we found that fibroblasts are predominantly present in the early stages of cell differentiation. Similarly, the ncWNT pathway, as well as the genes WNT5A and MCAM, exhibit the highest expression levels in the early stages of differentiation. This suggests that the binding of the WNT5A ligand factor to the MCAM receptor factor in the fibroblasts subpopulation affects the transmission of the ncWNT signal, which may play a significant role in the occurrence and development of tumors. This result is consistent with the analysis results of PCP (Fig S4A–C).

To delve deeper into the functions of the molecules that interact within the signaling pathway, an Enrichr analysis was initially performed on the ligand-receptor pairs found in epithelial and fibroblast cells. The results revealed significant enrichment in pathways such as the PI3K-Akt signaling pathway, Focal Adhesion, Focal Adhesion-PI3K-Akt-mTOR signaling pathway, TNF-alpha Signaling via NF-kB, Hypoxia, Glycolysis, Epithelial Mesenchymal Transition, and Angiogenesis in both epithelial and fibroblasts (Fig. 7A, B and Supplementary Table 9). Following this, a “SCENIC” analysis was performed on the fibroblast population. The results of the analysis indicatedthat the transcription factor activities of EGR1 and HMGB1 within the Fibroblast_1 subgroup were the most elevated (Fig. 7C). EGR1 transcription factor modulates the expression of numerous downstream genes that are implicated in cellular differentiation, proliferation, and inflammatory responses [19, 20]. It is well documented that in cases of acute kidney injury, EGR1 enhances the expression of SOX9 in renal tubular epithelial cells (TECs) by directly binding to the promoter region of the SOX9 gene [21]. This action promotes the proliferation of SOX9-positive renal tubular cells and boosts the kidney’s repair capabilities via the activation of the Wnt/β-catenin pathway [21]. Consequently, it is hypothesized that EGR1 may facilitate tumor progression by modulating downstream target genes, thereby activating the WNT/PI3K/AKT signaling pathway. To sum up, our results collectively demonstrate that pathways such as the PI3K-Akt signaling pathway and WNT signaling pathway play a crucial role in facilitating communication between epithelial cells and fibroblasts.

Fig. 7

Enrichment analysis of ligand-receptor pairs between epithelial cells and fibroblasts, and SCENIC analysis of fibroblast subpopulations. A, B Enrichr analysis of ligand-receptor pairs in epithelial cells and fibroblasts. C SCENIC analysis of fibroblast subpopulations

Incoming patterns of fibroblasts acting as ligands in the signaling communication with neutrophils

Subsequently, we applied a similar strategy as with the outgoing signals and utilized CellChat to analyze the incoming patterns in PCP samples. It was observed that Neutrophils were the primary recipients of incoming signals, followed by monocytes and macrophages (Fig. 8A). By analyzing the Cophenetic and Silhouette indices of incoming signaling, cell types and signaling pathways were categorized into five distinct patterns (Fig. 8BC). Pattern 1 included Epithelial_1, Epithelial_2, Epithelial_3, and Epithelial_4. Pattern 2 comprised Monocytes, DCs, and Macrophages. Pattern 3 consisted of Fibroblast_2, Fibroblast_1, Fibroblast_3, and Fibroblast_5. Pattern 4 encompassed Endothelial cells and Fibroblast_4. Pattern 5 primarily involved Neutrophils. In communication Pattern 5, the most influential factors were IL1, CXCL, SAA, ANNEXIN, and CSF3 (Fig. 8C and Supplementary Table 10, 11).

Fig. 8

Incoming patterns of fibroblasts acting as ligands in the signaling communication with neutrophils. A Incoming patterns in PCP samples revealed by scRNA-seq. B, C Cell types and signaling pathways were classified into five distinct patterns through the analysis of the Cophenetic and Silhouette indices of incoming signaling. D Enrichr analysis result using Neutrophils as the receptor and other cells as ligands. E Enrichment analysis result using fibroblasts as the ligand and Neutrophils as the receptor

Given that Neutrophils were identified as the primary recipients of incoming signaling, an Enrichr analysis was conducted with Neutrophils as the receptor and other cells as ligands. The analysis revealed significant enrichment in Epithelial Mesenchymal Transition and immune-related pathways such as the Cytokines and Inflammatory Response, Photodynamic therapy-induced NF-kB survival signaling, TNF-alpha Signaling via NF-κB Inflammatory Response, and IL-6/JAK/STAT3 Signaling (Fig. 8D and Supplementary Table 12). Subsequently, another Enrichr analysis was carried out with fibroblasts as the ligand and Neutrophils as the receptor. The results of this analysis revealed enrichment in pathways similar to those observed when using all other cells as ligands (Fig. 8E and Supplementary Table 13). In summary, our results suggest that there is extensive communication of immune signal pathways between fibroblasts and neutrophils, indicating the potential for immunosuppressive therapy targeting fibroblasts.

Important ligand-receptor pairs between fibroblasts and neutrophils

Our study revealed that neutrophils are the receptors with the most significant contribution to the communication between fibroblasts and immune cells. To elucidate the signaling molecules between fibroblasts and neutrophils, the transmission probability between these cell types was evaluated. The top outgoing signals from fibroblasts to neutrophils were identified as ANXA1-FPR1 and MIF-(CD74 + CXCR2) (Fig. 9A and Supplementary Table 14). Subsequently, an analysis of the ANNEXIN signaling network with neutrophils as receptors was conducted. The data presented in Fig. 9B clearly demonstrated the strong communication between the subgroups of epithelial cells and fibroblasts acting as ligands and neutrophils within the ANNEXIN signaling network. Moreover, a contribution analysis of the regulatory factors within the pathway highlighted that ANXA1-FPR1 makes the most significant contribution (Fig. 9C). To address the potential for overfitting and biases associated with the computational methods applied, we conducted independent analyses of cell–cell communication for the fibroblast subpopulations and neutrophils via the “CellChat” method. Moreover, the ligand-receptor pair ANXA1-FPR1 was pinpointed as the most prominent in these interactions (Fig. S5). Therefore, our results suggest that ANXA1 and FPR1 could potentially serve as targets for immunotherapy in the treatment of PCP. To validate this analysis, we also conducted a “CellChat” analysis in GBM samples. The analysis revealed that ANXA1-FPR1 is the only ligand-receptor pair when fibroblasts act as ligands and neutrophils act as receptor cells, interacting with each other. This result is consistent with the analysis results of PCP (Fig S6A, B).

Fig. 9

Important ligand-receptor pairs between fibroblasts and neutrophils. A The transmission probability between fibroblasts and neutrophils. B The ANNEXIN signaling network with neutrophils as receptors. C ANXA1-FPR1 pair emerged as the most significant regulatory factor in the pathway. Enrichr analysis on all ligand-receptor interactions between fibroblasts (D) and neutrophils (E).

Subsequently, an Enrichr analysis was conducted on all ligand-receptor interactions between fibroblasts and neutrophils, revealing significant enrichment in pathways such as Cytokines and Inflammatory Response, KRAS Signaling Up, Interferon Gamma Response, IL − 6/JAK/STAT3 Signaling, Glycolysis, Allograft Rejection (Fig. 9DE and Supplementary Table 15). Chronic inflammation, the infiltration of immune cells, and the ability of cancer cells to evade the immune response are regarded as significant characteristics in the progression of cancer [22]. CAFs, as the main components of the TME, have been reported to primarily function in immunosuppression to assist in cancer immune evasion [23]. Therefore, we speculate that fibroblasts in PCP should also have an immunosuppressive effect, which can influence the function and role of immune cells in the TME.

The correlation of ligand-receptor genes in the Wnt pathway in PCP samples

The communication and interactions between cancer cells and fibroblasts are crucial in tumor metastasis and progression [16]. In our study, scRNA-seq revealed that the WNT5A ligand and MCAM receptor had the highest contributions between epithelial cells and fibroblasts. To investigate the expression of ligand and receptor genes associated with the Wnt signaling pathway, bulk RNA sequencing was conducted on 11 PCP samples and five normal controls. A differential analysis was then carried out on the original count matrices of the normal and tumor samples. As depicted in Fig. 10A, the ligand gene (WNT5A) and receptor genes (MCAM, FZD1, FZD6) identified by both scRNA-seq and bulk RNA sequencing in the Wnt pathway exhibited significant differences between normal and tumor samples, with all showing markedly increased expression in the tumor group. Subsequently, correlation analysis was performed on the ligand-receptor pairs, revealing strong correlations in all pairs (Cor > 0.5). WNT5A and MCAM showed the strongest correlation (r = 0.951), with other receptor genes also showing a positive correlation with the ligand (Fig. 10B–D). We downloaded bulk RNA data for GBM samples from the GEO database and found through correlation analysis that WNT5A and MCAM have the highest correlation (Fig S7A–D).These results of the bulk RNA sequencing analysis were consistent with those of scRNA-seq, indicating that the Wnt pathway plays a crucial role in the cellular communication between epithelial cells and fibroblasts.

Fig. 10

Correlation of ligand-receptor genes in the Wnt pathway in PCP samples using bulk RNA sequencing. A Heatmap of Wnt pathway gene expressing pattern of genes identified by both scRNA-seq and bulk RNA sequencing. B–D Correlation analysis showed strong correlation between the ligand-receptor pairs of Wnt pathway

Discussion

PCP originates from embryonic remnants of the Rathke's pouch epithelium, leading to the formation of a differentiated squamous cell network [24]. In this study, scRNA-seq showed a significant number of epithelial cells in all PCP samples, aligning with previous research findings [25]. This suggests that epithelial cells play an important role in the pathogenesis of PCP. The TME typically comprises extracellular matrix (ECM) and various cell types. Previous studies have indicated that tumor development is not solely attributed to the malignant transformation of a specific cell type, but rather to the synergistic interaction between cancer cells and their stroma [26]. Normal cells within the TME are vital in tumor development, infiltration, and invasion [26]. Although CPs are benign tumors, inflammatory cell infiltration is often present around them due to their invasiveness and attachment to areas such as the hypothalamus [24]. Correspondingly, we observed a substantial infiltration of immune cells in PCP samples in this study. CAFs play a crucial role in tumor metastasis and progression [16].Interestingly, fibroblasts were not detected in ACP; however, they were observed in PCP in this study. The InferCNV analysis revealed a significant presence of CNVs in fibroblasts, indicating that fibroblasts may play a role in the development and progression of PCP tumors. Therefore, our focus shifted to investigating the role of fibroblasts in the progression of PCP.

Studies suggest that ACPs arise due to somatic mutations in CTNNB1. This gene encodes β-catenin and enhances its stability, leading to the dysregulation of the WNT pathway [27]. We hypothesize that the WNT pathway may serve a similar function in the PCP subtype as it does in ACP. As a result, we subsequently performed an “Enrichr” analysis between epithelial cells and fibroblasts. Our findings suggest that, like in ACP, the Wnt pathway is also significantly activated in PCP tumor tissues. The Wnt signaling pathway includes both canonical and non-canonical routes [28]. The canonical Wnt pathway predominantly controls cell proliferation, whereas the non-canonical Wnt pathway governs cell polarity and migration. The Wnt signaling pathway is instrumental in the self-renewal of specific mammalian tissues [29]. For example, it is involved in the development and regeneration of the small intestinal epithelium and facilitates the differentiation of Paneth cells located at the crypt bases [30]. Additionally, the Wnt signaling pathway is closely linked to processes such as liver metabolism and regeneration, lung tissue repair and metabolic functions, hair follicle cycling, hematopoietic system development, osteoblast maturation, and activity [31,32,33]. Consequently, we hypothesize that in PCP, the worsening of the disease is mainly attributable to the tumor cells’ proliferative responses mediated by the Wnt pathway. CellChat is widely employed to accurately represent cell–cell signaling interactions and perform effective systems-level analyses of these interactions [34]. To explore the intercellular communication within PCP, we performed a “CellChat” analysis. It is revealed that epithelial cells played a substantial role in transmitting output signals among different cell subtypes,the majority of epithelial cell clusters collectively contribute to the formation of pattern 1, with prominent contributors including GDF, CSF3, ncWNT, BAG, SAA, PTN, MK, WNT, CX3C, PERIOSTIN, and SEMA3. The most critical ligand-receptor pairs identified between Pattern 1 and Epithelial_4 with fibroblasts were WNT5A-MCAM/FZD4 and VEGFA-VEGFR2/VEGFR1R2. Upon further analysis of the ncWnt signaling patterns, it became apparent that the WNT5A-MCAM pair emerged as the predominant regulatory factor in the pathway. Bulk-RNA sequencing data from 11 PCP samples and 5 normal controls revealed significant correlations between the ligand gene (WNT5A) and receptor genes (MCAM, FZD1, FZD6). Tsai and colleagues have demonstrated that active β-Catenin is contained within exosomes secreted by cells that express the atypical WNT ligand WNT5A, thereby eliciting β-Catenin-dependent responses in recipient cells [35]. In our analysis, we discovered that the epithelial cell subpopulation within the ncWNT pathway acts via an autocrine signaling mechanism, whereas the fibroblast subpopulation engages in paracrine signaling. Consequently, we hypothesize that the epithelial cell subpopulation releases the WNT5A ligand, which then activates the β-Catenin protein in both the epithelial cells themselves and in fibroblasts acting as receptor cells. This leads to the abnormal activation of the canonical WNT pathway and downstream target genes, such as MMPs, thereby promoting the abnormal proliferation of epithelial cells and ultimately contributing to tumor formation. Notably, a correlation analysis was conducted across seven neurological disorders, which revealed that key pathways and protein molecules akin to PCP were exclusively identified in GBM. The distinct nature of PCP is underscored by this finding, as is the critical need for more targeted analyses dedicated to the specific understanding of PCP.

To validate our hypothesis, we next performed an “Enrichr” analysis on the key ligand-receptor pairs in both epithelial and fibroblast cells. The results revealed a significant enrichment of the PI3K-Akt signaling pathway in PCP samples, confirming the close relationship between the canonical WNT pathway and the PI3K-Akt signaling pathway. The aberrant activation of these pathways aids tumor cells in evading apoptosis and promotes their proliferation, survival, and migration. Interestingly, a strong connection between the PI3K/AKT/mTOR and WNT/β-Catenin signaling pathways has been documented in various cancers, including colorectal cancer (CRC) [36], pancreatic cancer [37], lung cancer, breast cancer, and hepatocellular carcinoma (HCC), as well as in glioblastoma multiforme (GBM) [38, 39]. Prossomariti and colleagues reported that the inhibition of the WNT/β-Catenin signaling pathway in CRC is associated with the activation of the PI3K/AKT/mTOR signaling cascade. Conversely, the inhibition of the PI3K/AKT/mTOR pathway results in the abnormal activation of the canonical WNT pathway [40]. Our analysis, therefore, indicates that the crosstalk between the PI3K/AKT/mTOR and WNT/β-Catenin pathways in PCP is a key factor in tumor development and progression. Consequently, targeting the ligand-receptor pairs of the WNT/β-Catenin and PI3K/AKT/mTOR pathways could offer novel therapeutic strategies for PCP. Aditionally, the signaling of WNT5A leads to the depalmitoylation of pro-metastatic cell adhesion molecules CD44 and MCAM by acyl protein thioesterase 1 (APT1), ultimately enhancing melanoma invasion [41]. In gastric cancer, the downregulation of WNT5A protein through antibody treatment inhibited cell migration in vitro and suppressed cell invasion and metastasis in vivo [42]. Similarly, silencing the WNT5A gene and reducing protein expression in PC3 prostate cancer cells resulted in a decreased migratory ability of the cells [43]. Based on this, we propose that the enhanced WNT5A-MCAM pair between epithelial cells and fibroblasts plays a key role in promoting the malignancy of PCP tumors. Moreover, our results indicate that investigating WNT5A as a target for research could aid in developing precise therapies against particular molecular targets in PCP.

In another one of our studies (to be published), we discovered that the epithelial cells Epithelial_4 exhibit a notable degree of malignancy. It is well recognized that activated VEGFA-VEGFR2 can trigger angiogenesis and enhance glioma growth [44]. Furthermore, ACE2 inhibiting breast cancer angiogenesis by suppressing the VEGFA/VEGFR2/ERK pathway, leading to the reduction of VEGFA-VEGFR2 expression and phosphorylation [45]. Building upon our findings and existing research, we suggest that increased signaling between Epithelial_4 and fibroblasts promotes angiogenesis and accelerates the progression of PCP. Therefore, we hypothesize that communication through these pathways between epithelial cells and fibroblasts may contribute to the development of various cancer types. However, confirming this hypothesis will necessitate extensive research in the future.

In our previous study of bulk-RNA sequencing, we observed a higher abundance of neutrophils in the PCP tumor tissues compared with the normal brain tissues [15]. This study revealed that incoming signals were predominantly received by neutrophils, indicating a significant role for these cells in PCP tumor tissues. Additionally, our earlier findings showed a strong correlation between IL1A and the invasiveness of PCPs [15]. Consistent with this, IL1 was identified as the most influential factor in this study. When considering neutrophils as receptors, along with other cell types and specifically fibroblasts as ligands for functional enrichment, the analysis highlighted enriched pathways related to Epithelial Mesenchymal Transition and immune pathways such as TNF-alpha Signaling via NF-κB Inflammatory Response. As vital components of innate immunity and inflammation, neutrophils not only secrete proteases that degrade tissue but also fuel inflammation by releasing cytokines and chemokines [46]. Tumor necrosis factor alpha (TNFα) is a pivotal mediator of this inflammatory response. Helen et al.'s research revealed that in rheumatoid arthritis, neutrophils exhibit heightened TNFα expression, and the activation of the transcription factor NF-κB (a TNFα target) may trigger a self-perpetuating inflammatory cycle [47]. Therefore, the chronic inflammation observed in PCP may be linked to the TNFα-activated NF-κB pathway, suggesting that immunotherapy targeting neutrophils could be a beneficial strategy for PCP treatment.

ANXA1-FPR1 and MIF-(CD74 + CXCR2) were identified as the primary outgoing signals from fibroblasts to neutrophils. It is noteworthy that a robust ANNEXIN signaling network was also observed between epithelial cells and neutrophils. Annexin A1 (AnxA1) is known for its role in promoting metastasis and angiogenesis, with its expression positively correlating with the progression of various cancer types [48]. The ANXA1/FPR1 signaling axis is recognized for its ability to suppress the innate immune response in glioma patients by creating an anti-inflammatory environment and supporting a Treg-driven TME [49]. Additionally, it is reported that blocking the autocrine axis of AnxA1/FPR1 leads to decreased growth and aggressiveness of MDA-MB-231 breast cancer cells both in vitro and in vivo. Studies have shown that blocking the autocrine axis of AnxA1/FPR1 leads to reduced growth and aggressiveness of MDA-MB-231 breast cancer cells both in vitro and in vivo [48]. Therefore, we propose that the AnxA1/FPR1 pathway, connecting epithelial cells, fibroblasts, and neutrophils, plays a role in promoting the increased growth and aggressiveness of PCP cancer cells. Our previous findings demonstrated a positive correlation between IL6 levels and the invasion of PCP tumors in the hypothalamus [15]. Consistently, Enrichr analysis of all ligands and receptors interactions between fibroblasts and neutrophils revealed a significant enrichment of immune-related pathways, including the IL-6/JAK/STAT3 signaling pathway. IL6 is a multifunctional cytokine whose levels in the central nervous system (CNS) are reported to rise in conditions such as cerebral ischemia, brain injury, and neurodegenerative diseases [50]. The overexpression of IL6 in the brains of transgenic mice has been linked to dysfunctions of the blood–brain barrier, neuronal damage, and neuroinflammation [51, 52]. Brain pericytes are a TNF-α-sensitive cell type that acts as a specific effector of TNF-α during brain inflammation by releasing IL6 [53, 54]. TNF-α can activate multiple intracellular signaling pathways, including those involving the Janus family of tyrosine kinases (JAK) signaling transducers, the JAK-STAT3 pathway, and the inhibitor kappa B (IκB)-nuclear factor kappa B (NFκB) pathway [53, 54]. The crosstalk between the JAK-STAT3 and IκB-NFκB pathways promotes the production of IL6 and provides positive feedback for IL6 signaling [55, 56]. This interplay between IL6-mediated signaling pathways is likely a significant contributing factor to the progression of PCP.

In conclusion, this study utilized scRNA-Seq data to explore the fundamental interactions among cells in PCP, offering a comprehensive understanding of the roles of ligand-receptor interactions in signaling pathways. We have identified that the ligand and receptor genes in the WNT/PI3K/AKT signaling pathway, AnxA1/FPR1, and TNF-α/IL-6/JAK/STAT3 signaling pathway may have potential prognostic significance. The expression of these regulatory factors at the single-cell level reveals the molecular mechanisms driving biological processes, which could assist in identifying drug targets for PCP treatment. It is essential to acknowledge the limitations of our study. The shortage of subsequent validation studies has significantly hampered our ability to validate the analytical results, thereby diminishing the generalizability of our research conclusions. This is mainly attributed to the present scarcity of commercial PCP cell lines and the challenges involved in the isolation and characterization of primary cells. To ascertain the robustness of our findings, we conducted a thorough analysis of single-cell and bulk RNA sequencing data for PCP. Additionally, we retrieved seven datasets pertaining to various neurological disorders from the GEO database to perform comprehensive validation analyses. Moving forward, we stand prepared to collect tissue samples with the intent of identifying key molecular markers. Our objective will be achieved by leveraging a comprehensive suite of analytical strategies, including quantitative polymerase chain reaction (qPCR), western blotting, and immunohistochemical techniques. These selected approaches are designed to construct a thorough molecular profile, thereby refining our understanding of the underlying pathobiological mechanisms of PCP. Additioanlly,, the small sample size used in this study may limit the broader applicability of our therapeutic insights for PCP. In our future work, we aim to continue collecting samples to corroborate the findings of this analysis, which will enhance our understanding of PCP and provide more comprehensive treatment guidelines for clinical use.