Abstract

   Despite improvements in microscopically neurosurgical techniques made
   in recent years, the prognosis of adamantinomatous craniopharyngioma
   (ACP) is still unsatisfactory. Little is known about cellular atlas and
   biological features of ACP. Here, we carried out integrative analysis
   of 44,038 single-cell transcriptome profiles to characterize the
   landscape of intratumoral heterogeneity and tumor microenvironment
   (TME) in ACP. Four major neoplastic cell states with distinctive
   expression signatures were defined, which further revealed the
   histopathological features and elucidated unknown cellular atlas of
   ACP. Pseudotime analyses suggested potential evolutionary trajectories
   between specific neoplastic cell states. Notably, a distinct
   oligodendrocyte lineage was identified in ACP, which was associated
   with immunological infiltration and neural damage. In addition, we
   described a tumor-centric regulatory network based on intercellular
   communication in TME. Together, our findings represent a unique
   resource for deciphering tumor heterogeneity of ACP, which will improve
   clinical diagnosis and treatment strategies.
     __________________________________________________________________

   The landscape of cell atlas, tumor heterogeneity and intercellular
   communication network revealed by scRNA-seq in ACP.

INTRODUCTION

   Adamantinomatous craniopharyngiomas (ACPs) are rare intracranial tumor
   that tend to occur in the sellar region and originates from embryonic
   malformation of the craniopharyngeal epithelium (a remnant of the
   ectoderm, also known as Rathke’s pouch) ([76]1, [77]2). Epidemiological
   statistics support a bimodal age distribution of patients with ACP,
   mainly in children aged 5 to 14 years and adults aged 50 to 74 years
   ([78]2). Even extremely rare ACP cases occur in the prenatal and
   neonatal periods ([79]3–[80]5). ACPs are typically onset in childhood
   and adolescence and account for approximately 6 to 9% of pediatric
   brain tumors ([81]6, [82]7). As the most common lesion involving the
   hypothalamic-pituitary axis in childhood, ACP causes catastrophic
   damage to endocrine homeostasis, growth, and development in patients
   ([83]8, [84]9). Histopathologically, ACP is characterized by cysts,
   typically prominent calcifications and cholesterol crystals ([85]2,
   [86]10). Finger-like protrusions of tumor epithelium form at the
   junction of tumor and normal tissue and invade surrounding important
   neural structures ([87]6). This aggressive biological behavior makes
   treatment challenging and contributes to the poor prognosis ([88]11).
   Surgery is still the first-line treatment for ACP. However, considering
   the aggressiveness of the tumor and the close relationship between the
   tumor and the position of the pituitary stalk, optic chiasm, and
   especially hypothalamus, gross total resection is extremely challenging
   and often leads to serious neurological impairment and poor prognosis
   ([89]12–[90]14). Moreover, partial resection often causes frequent
   tumor recurrence ([91]15). Therefore, surgery alone cannot solve the
   clinical problems caused by ACP, and it is necessary to develop a
   precision treatment that can effectively eliminate the tumor and reduce
   the damage to the surrounding normal structures. There is an urgent
   need for further in-depth understanding of the biological
   characteristics of ACP and the molecular mechanisms behind the origin
   and development of ACP.

   Because of the lack of mature cell lines and organoid models, basic
   research on the tumor biology of ACP and its tumor microenvironment is
   progressing slowly. Previous studies supported the complex histological
   structure of ACP, while its cellular composition and functional
   characteristics are still unclear, which hinders our in-depth
   exploration of its tumorigenesis mechanism. Pathological staining can
   reveal the special internal organization of ACP, composed mainly of
   palisade-like epithelial cell layers, helical cell clusters, keratin
   nodules, and stellate reticulum ([92]16, [93]17). This diversity of
   histological structures may lead to marked differences in cellular
   function. ACP is generally considered to be driven by mutations in exon
   3 of the Catenin Beta 1 (CTNNB1) gene ([94]18–[95]21). Activating
   mutations of CTNNB1 prevent its encoded protein β-catenin from being
   degraded normally, leading to its accumulation in the nucleus and
   inducing persistent activation of the WNT pathway. Previous studies
   have suggested that the nuclear accumulation of β-catenin caused by
   constitutive mutations, such as S37F, may be a key driving factor for
   tumorigenesis ([96]22–[97]24). In a study using a mouse model of
   spontaneous ACP, SOX2^+ adult pituitary stem cells harboring a
   degradation-resistant form of β-catenin successfully induced
   tumorigenesis with features of ACP ([98]25). However, nuclear
   accumulation of β-catenin was seen only in sporadic cells, while the
   majority of cells exhibited normal membrane localization of β-catenin,
   suggesting the potential tumor heterogeneity of ACP ([99]26, [100]27).
   Furthermore, therapeutic strategies targeting β-catenin have not been
   efficacious, suggesting that there may be multiple underlying factors
   involved in the process of tumor progression ([101]7).

   Chronic inflammation and functional impairment of the central nervous
   system, which are closely related to patient prognosis and quality of
   life and survival, have accompanied most of the course of ACP
   development. However, the molecular mechanisms remain unclear. Reactive
   glial scars, which arise during ACP invasion into brain tissue,
   surround the tumor and endow ACP with a unique tumor microenvironment
   (TME) containing diverse immune cells and glial cells
   ([102]28–[103]30). Insight into the function and clinical significance
   of immune cells in ACP tissues has been provided by our previous study
   ([104]31). In contrast to most brain tumors, the lack of restriction of
   the blood-brain barrier allows ACP to more frequently come into contact
   with infiltrating immune cells ([105]32–[106]34). Therefore, further
   comprehensive characterization of the immune landscape of ACP is
   needed. Furthermore, as the buffering zone for ACP to invade brain
   tissue, the glial scar tissue is enriched with abundant astrocytes and
   oligodendrocyte lineage cells (OLs), which play a potentially vital
   role in remodeling local tissue microenvironment that facilitates ACP
   progression. Considering the severe impairment of ACP in the central
   nervous system that is unavoidable, such as endocrine deficits (52 to
   87%) and visual impairments (62 to 84%), dissecting the interaction
   between the tumor and the glial scar is crucial ([107]6).

   It is difficult to identify specific cell types and their cross-talk by
   analyzing traditional bulk transcriptomics data when the tumor is
   interrogated on a whole level. The latest development of single-cell
   RNA sequencing (scRNA-seq) technology provides us with the ability to
   obtain a more detailed perspective of intratumor heterogeneity and the
   composition of the TME at single-cell resolution ([108]35). Therefore,
   we performed scRNA-seq and analyzed the transcriptomes of 27,696 single
   cells from three ACP tumor specimens and adjacent glial scars in this
   study. We further integrated scRNA-seq data of 8453 peripheral blood
   mononuclear cells and 7889 normal brain tissue cells from public
   databases to reveal the functional characteristics of nonneoplastic
   cells in the ACP TME. Our work innovatively describes the unique
   classification of neoplastic cells and the complex TME of ACP, which
   will not only promote our understanding of the composition and
   structure of ACP cells but also provide in-depth mechanistic insights
   into the diagnosis, treatment, and progression of ACP.

RESULTS

Single-cell expression atlas and cell type identification in ACP TME

   We collected tumor tissues derived from three patients with ACP to
   perform scRNA-seq ([109]Fig. 1A and table S1). First, our pathological
   examination confirmed that their tumor radiological imaging and
   immunostaining of anti–β-catenin and anti-LEF1 antibodies were in line
   with the classic features of ACP (fig. S1, A and B). In addition, we
   also performed bulk RNA and exome sequencing on the aforementioned ACP
   tumors to gain insight into their basic properties. The transcriptional
   signatures and activating mutations of CTNNB1 exhibited by these tumor
   tissues were consistent with the typical genomic pattern of ACP (fig.
   S1, C and D, and table S2) ([110]12, [111]36). Eventually, we obtained
   a total of 27,696 cells for scRNA-seq, and the average number of genes
   detected per cell was 2701 (fig. S1, E to I). Following rigorous
   quality control, we performed the clustering analysis by using an
   unsupervised method and obtained 18 cell subpopulations (fig. S2, A and
   B). According to the transcriptomic characteristics of the
   subpopulations, 27,696 single cells were further grouped into eight
   main clusters ([112]Fig. 1B). Cluster-specific marker genes were used
   to annotate the classic cell types based on previous studies, including
   CDH1 and EPCAM for epithelial cells, LYZ and CD68 for myeloid cells,
   CD3E for T cells, CD79A for B cells, AQP4 and GFAP for astrocytes,
   OLIG2 and SOX10 for oligodendrocyte precursor progenitor cells (OPCs),
   MOG for oligodendrocytes, and PECAM1 for endothelial cells
   ([113]Fig. 1, B and C, and fig. S2B) ([114]12, [115]25). It is worth
   noting that the rare endothelial cells may reflect biases in tissue
   dissociation and microfluidic-based cell capture. Compared with other
   samples, the ACP-2, which was obtained from the junction of tumor and
   glial scar, contained more immune cells and glial cell components.
   Multiplex immunohistochemistry (IHC) staining identified the
   distribution of the main cell types, including CDH1 for epithelial
   cells, CD3 for T cells, CD68 for myeloid cells, GFAP for astrocytes,
   OLIG2 for the OLs, and the tumor-specific mutated gene CTNNB1 in the
   TME ([116]Fig. 1D). Furthermore, we identified neoplastic cells by two
   complementary approaches. First, previous studies have proposed that
   ACP is a type of tumor composed of squamous epithelial cells ([117]37).
   Multi-IHC staining indicated that the epithelial cells located in
   CDH1-positive staining area conformed to the classic ACP pathological
   characteristics, permitting us to spatially locate and further resolve
   neoplastic cells in ACP ([118]Fig. 1D and fig. S2C). Second, by
   performing inferCNV analysis, we found higher levels of copy number
   variation (CNV) for epithelial cells as compared with immune cells and
   glial cells ([119]Fig. 1E). As a result, 23,535 cells were defined as
   neoplastic cells. To verify the accuracy of the identification of
   neoplastic cells, we curated gene signatures of ACP from previous
   studies based on the bulk transcriptome data and examined their
   expression levels in various subpopulations ([120]19, [121]38).
   Subsequently, we showed substantially higher expression levels of
   canonical ACP marker genes in all the neoplastic cell clusters,
   although we noted a biased distribution (fig. S2D). Immunostaining
   results showed that most immune cells and glial cells in the TME of ACP
   were located outside of tumor structure, which were composed mainly of
   reactive glial scar tissue (fig. S2E).

Fig. 1. Single-cell transcriptomic analysis identifies cell types in ACP.

   [122]Fig. 1.
   [123]Open in a new tab

   (A) Workflow demonstrates the process of sample collection, single-cell
   suspension preparation, single-cell sequencing, and
   bioinformatics-based analysis to reveal tumor heterogeneity. (B) UMAP
   plots of 27,696 cells showing the major cell types and sample origins
   in ACP tissues based on analysis of 10X GENOMICS scRNA-seq data. (C)
   Violin plots show the expression levels of specific marker genes for
   eight major cell types. (D) Multi-IHC staining with anti-CDH1
   (epithelial cells), anti-CTNNB1 (ACP classic mutant gene), anti-CD68
   (macrophages), anti-CD3 (T cells), anti-OLIG2 (oligodendrocyte
   lineages), and anti-GFAP (astrocytes) antibodies to localize the
   histological distribution of major cell types in ACP. Scale bars, 50
   μm. (E) Heatmap showing the chromosomal map of single-cell large-scale
   CNVs inferred by scRNA-seq. Use immune cells as a reference, and
   amplifications (red) or deletions (blue) were inferred by averaging
   expression over 100-gene stretches on every chromosome.

Characterization of transcriptional heterogeneity of neoplastic cells in ACP

   The initial aggregation of neoplastic cells revealed patient-specific
   clusters (fig. S3A). Tumors derived from individual patients showed
   preferential expression of more than 150 genes, manifested by different
   levels of enrichment of classic tumor pathways and metabolic pathways
   (fig. S3, B to D). To characterize the intratumor heterogeneity of
   neoplastic cells, we analyzed clusters of coexpressed genes in each
   tumor based on nonnegative matrix factorization (NMF) and obtained a
   total of 25 intratumor signatures (table S3). Taking ACP-3 as an
   example, the NMF results demonstrated complex intratumoral expression
   features in ACP ([124]Fig. 2A). By assuming that some intratumoral
   signatures could reflect common features of ACP, we subsequently
   identified four prominent meta-signatures across tumors by performing
   correlation cluster analysis of 25 signatures ([125]Fig. 2B). These
   four meta-signatures contained their own unique gene sets ([126]Fig. 2C
   and table S4). The first meta-signature was related to epithelial
   development and was composed mainly of keratin genes, such as KRT6C and
   KRT7, as well as cell adhesion molecules, such as LGALS3 and CLDN4. The
   second meta-signature was composed of WNT signaling pathway molecules,
   such as LEF1 and WNT5A, and downstream target genes, such as FGF4 and
   BMP4. Most of genes in the third meta-signature were related to cell
   cycle, which was composed of mitotic genes, DNA replication–, and
   transcription-related genes such as UBE2C, CENPK, and TOP2A. Genes in
   the fourth meta-signature were related to stress responses and cellular
   activation, such as FOS, JUN, and EGR2 ([127]Fig. 2C). Next, we
   regrouped all 23,535 neoplastic cells into 10 clusters through Seurat
   integrated clustering (fig. S4A). Subsequently, on the basis of the
   abovementioned marker gene sets of meta-signatures, we calculated
   functional scores for all neoplastic cells, which permitted us to
   lastly define four major cell states ([128]Fig. 2D and fig. S4B).
   According to previous studies, ACP has unique histopathological
   structures, such as palisading epithelium, whorl-like structures, and
   stellate reticulum ([129]19). We hypothesized that neoplastic cells
   with similar transcriptional characteristics and functions tended to be
   spatially located in neighboring areas. Immunostaining of specific
   markers revealed histologically unique locations of different cell
   states ([130]Fig. 2E). According to the histological distributions and
   molecular functions of the neoplastic cells, we defined the four cell
   states with specific meta-signatures as: (i) keratinized-like
   epithelium (KE) with epithelial development feature, (ii) whorl-like
   epithelium (WE) with WNT feature, (iii) cycling cell (CC) with cell
   cycle feature, and (iv) palisade-like epithelium (PE) with stress
   feature. PE was distributed mainly along the outermost side of the
   tumor tissue, arranged closely in a fence-like manner, and was in
   contact with the outer stromal tissue. WE was distributed mainly in
   epithelial whorl-like structures, showing the expression pattern of
   nuclear accumulation of β-catenin. KE expressing a unique spectrum of
   cytokeratins was mainly distributed along the tumor inner side in which
   the stellate reticulum was located. CC was scattered mainly in the PE
   but rarely in other structures. In stark contrast, neoplastic cells
   clustered distribute in ACP, multiple types of cells grow intertwined
   each other in other kinds of tumors. Furthermore, we validated this
   distribution pattern in additional 18 ACP tumor specimens by IHC
   staining (fig. S5). These cellular states appeared to reflect essential
   features of ACP during tumor development. Although ACP-2 contained only
   a few neoplastic cells, the proportions of the four cell states were
   still similar to those of other samples. We also found that epithelial
   whorl-like structures, which were previously implicated in
   mutation-induced persistent activation of the WNT pathway, were
   composed of two cell states, WE and KE. Only WE shown nuclear
   accumulation of β-catenin (fig. S4C).

Fig. 2. Deciphering expression programs revealed distinct clusters of
neoplastic cells in ACP.

   [131]Fig. 2.
   [132]Open in a new tab

   (A) Heatmap showing gene expression signatures parsed from a
   representative tumor (ACP-3) using NMF. Each expression signature was
   numbered, and its corresponding genes were displayed (right). (B)
   Pearson correlation clustering of 25 intratumoral expression
   signatures. Four regions consisting of clusters of highly correlated
   signatures were identified as meta-signatures. (C) Heatmap showing the
   expression levels of meta-signatures marker genes in each intratumoral
   signature. (D) UMAP plots showing the distribution of cell states (a),
   clusters (b), and sample origins (c) in ACP neoplastic cells. Ten
   neoplastic clusters (b) were generated by Seurat analysis and
   summarized into four cell states (a) based on the AddModuleScore
   scores. (E) Multi-IHC showing the histopathological distribution of CC
   (anti-Ki-67, TK1, and UBE2C), WE (anti–β-catenin, MUC1, and ITGA1), KE
   (anti-KRT7, KRT6, and CALB1), and PE (anti-CD47, KRT7, and GFAP). Scale
   bars, 50 μm. (F) Heatmap showing differences in gene ontology (GO)
   pathway activity between neoplastic clusters of ACP scored by GSVA. (G)
   Heatmap showing differences in TF activity across neoplastic clusters
   of ACP scored by SCENIC. (H) Cartoon model showing the spatial
   distribution of major neoplastic cell states in ACP tissue.

   To illustrate the functional heterogeneity of neoplastic cells, we
   carried out clustering analysis via gene set variation analysis (GSVA)
   to infer the biological functions of neoplastic clusters ([133]Fig. 2F
   and fig. S4D). GSVA revealed that the WE was related to activation of
   the WNT signaling pathway, organ development, and dentinogenesis.
   Pathways such as epithelial development, cell adhesion, cell movement
   and migration, and cell polarity regulation were up-regulated in KE.
   Endoplasmic reticulum stress and inflammation were enriched in PE,
   which may be due to stronger inflammatory stimulation presented in the
   TME of ACP.

   Furthermore, we sought to correlate transcriptome differences with
   transcription factor (TF) motifs across all neoplastic clusters by
   inferring a single-cell regulatory network and conducting clustering
   (SCENIC) analysis. SCENIC analysis identified a set of TFs associated
   with biological signatures of distinct neoplastic clusters in ACP. The
   results showed that cell clusters within the same cell states exhibited
   similar patterns of TF regulation ([134]Fig. 2G). In addition, we found
   that the key transcription motifs of LEF1, TCF4, and CTNNB1 in the WNT
   signaling pathway were activated in WE_1 and WE_2 (fig. S4E). Previous
   studies demonstrated that the constitutive activating mutation of
   CTNNB1 is responsible for activating the WNT pathway in ACP, and this
   mutation covers almost all neoplastic cells ([135]22). In this study,
   WNT-related transcriptional motifs were substantially activated in WE
   relative to other cell states. Immunostaining confirmed that the
   β-catenin and LEF1 were specifically expressed in WE and showed nuclear
   accumulation (fig. S1B). In addition, we found that the osteogenesis-
   and dentinogenesis-related TFs RUNX2 and MSX1 were substantially
   activated in WE ([136]Fig. 2G). Motifs of MYC and MAF were activated in
   the KE ([137]Fig. 2G), thereby regulating downstream cell adhesion
   molecules such as CLDN3 and CLDN4 as well as the malignant tumor
   progression-related gene MGP ([138]39, [139]40). In short, we
   identified four distinct neoplastic cell states with characteristic
   transcriptional regulatory networks and regular histological
   distribution, which provided further insight into the intratumoral
   heterogeneity of ACP ([140]Fig. 2H).

Pseudotime trajectory analysis of neoplastic cells in ACP

   Although the WE and KE revealed notable transcriptomic characteristics,
   their proliferative activities were low ([141]Fig. 2E and fig. S4C). In
   addition, there were some transcriptional features and TF motifs shared
   by clusters belonging to different cell states, suggesting that
   transformation of cell states might occur among clusters. To test this
   hypothesis, we further performed RNA velocity analysis to explore the
   transcriptional dynamics exhibited by subsets of neoplastic cells.
   Transcriptional rate analysis suggested that ACP neoplastic cells might
   start from PE and follow two different evolutionary directions. One of
   which evolved toward the WE cell states that was highly activated by
   the WNT signaling pathway. Whereas the other evolved toward the
   direction of KE, in which epithelial differentiation was more
   pronounced ([142]Fig. 3A). Immunostaining of KE and WE showed mutually
   exclusive histological distribution, supporting that they belonged to
   cell groups with different evolutionary paths (fig. S4C).
   Differentiation marker, including LEF1, TCF4, and WNT5A for WE and
   KRT6C, KRT7, and KRT75 for KE, also showed a high transcriptional
   activity in the evolutionary paths ([143]Fig. 3B). Notably, whorl-like
   neoplastic cells with nuclear accumulation of β-catenin have been
   implicated as drivers of ACP tumorigenesis in previous studies ([144]6,
   [145]38). Experimentally, up-regulated β-catenin in nuclear has induced
   spontaneous tumor formation similar to human ACP in transgenic mice
   ([146]41, [147]42). Unlike classical tumor germinal centers, whorl-like
   structures exhibit biological features of senescence, including cell
   cycle arrest and senescence-associated secretory phenotype ([148]25,
   [149]42). In addition, inferCNV analysis revealed a remarkably complex
   pattern of CNVs in WE (fig. S6A). According to previous classifications
   of known cellular states, CNV patterns were associated with a
   remarkable heterogeneity and functional states of neoplastic cells
   ([150]43, [151]44). We reconstructed patient-specific CNV patterns
   based on principal components analysis (PCA) dimensionality reduction
   and cluster analysis, which revealed a unique distribution of WE (fig.
   S6B). Next, the clonality tree was generated using the UPhyloplot2
   algorithm ([152]45, [153]46). A specific set of CNV accumulations
   related to WE was observed, such as gain of chr11q, chr14q, and chr20p,
   which may be potentially linked to an increase in the genomic
   instability of the aging state (fig. S6C). To further delineate gene
   expression patterns depicting the evolution of neoplastic cells, we
   performed Monocle3 analysis of two evolutionary paths separately
   ([154]Fig. 3, C and F). WNT-related genes showed a trend of positive
   regulation along the pseudotime axis of PE-WE evolution ([155]Fig. 3, D
   and E). Similarly, epithelial differentiation marker genes such as
   KRT6C, KRT7, KRT75, KRT13, and cell adhesion–related markers TM4SF1
   showed a positive regulatory trend in the PE-KE evolution axis with
   pseudotime ([156]Fig. 3, G and H). In summary, our results revealed
   potential evolutionary trajectories among neoplastic cell clusters in
   ACP.

Fig. 3. The evolutionary trajectory of neoplastic cells in ACP.

   [157]Fig. 3.
   [158]Open in a new tab

   (A) Inferred evolutionary trajectories of all ACP neoplastic cells by
   RNA velocity. (B) UMAP plots showing RNA velocity and expression levels
   of WE (LEF1, TCF4, and WNT5A) and KE (KRT6C, KRT7, and KRT75) marker
   genes. (C) Pseudotime evolution trajectories inferred by Monocle3 after
   reclustering of neoplastic cell states PE, WE, and CC. (D) Traceplots
   of WE marker genes expression levels along pseudotime. (E) Traceplots
   of the WNT signature score along pseudotime. (F) Pseudotime evolution
   trajectories inferred by Monocle3 after reclustering of neoplastic cell
   states PE, KE, and CC. (G) Traceplots of KE marker genes expression
   levels along pseudotime. (H) Traceplots of the epithelial development
   signature score along pseudotime.

T cell clustering and subtype analysis in ACP

   After integrating immune cells with published scRNA-seq data of
   peripheral blood mononuclear cells (PBMCs) derived from one healthy
   volunteer, we regrouped the cells to obtain an atlas of immune cells
   containing more than 11,000 cells ([159]Fig. 4A and fig. S7, A and B).
   Regrouping T cells identified a total of 10 clusters, of which three
   clusters were derived mainly from tumor tissues ([160]Fig. 4B and fig.
   S7, C to E). We further interrogated the expression of typical T cell
   markers in 10 clusters and scored immune function based on the gene
   expression of classic markers ([161]Fig. 4C and fig. S7F). The analysis
   further revealed that the CD4-C4-FOXP3 and CD4-C5-CTLA4 clusters
   represented CD25^+/FOXP3^+ regulatory T cells (T[regs]). CD4-C4-FOXP3
   originated mainly from peripheral blood, while CD4-C5-CTLA4 originated
   mainly from tumor tissues. The CD4-C5-CTLA4 cluster had the highest
   scores of T[reg] and immune cells exhaustion, indicative of its
   immunosuppressive state. The CD4-C3-CD69 cluster from tumors had a
   strong immune exhaustion score but did not express any T[reg] markers,
   which may represent a population of exhausted CD4^+ T cells in the TME.
   The CD8-C3-GZMK, CD8-C4-CX3CR1, and CD8-C5-IFNG clusters showed higher
   expression levels of cytotoxicity markers, such as GZMK, NKG7, and
   GZMH. However, the CD8-C5-IFNG cluster did not have the highest
   cytotoxicity score but instead had the most substantial exhaustion
   score.

Fig. 4. Assessing the functional states of tumor-infiltrating immune cells in
ACP.

   [162]Fig. 4.
   [163]Open in a new tab

   (A) UMAP plot showing the reclustering of all immune cells in ACPs and
   PBMCs after the integration of public data. CD4^+ T cells, CD8^+ T
   cells, myeloid cells, B cells, and natural killer (NK) cells were
   identified as the main cell types. (B) UMAP plot showing 10 clusters
   identified in T cell reclustering. (C) Violin plots showing naïve,
   exhaustive, cytotoxic, T[reg], and resident characteristics in CD4^+ T
   cells and CD8^+ T cells scored by AddModuleScore function. Dependent
   gene sets refer to table S5. (D) Volcano plot showing differentially
   expressed genes between neoplastic and nonneoplastic tissue–derived
   T[regs]. (E) Volcano plot showing differential genes of neoplastic and
   nonneoplastic tissue–derived cytotoxic T cells. (F) UMAP plot showing
   the reclustering of myeloid cells. (G) UMAP plots showing the
   expression and distribution of canonical myeloid marker genes in
   myeloid cells. (H) Violin plots showing representative M1 and M2
   characteristics in monocytes/macrophages scored by AddModuleScore
   function. Dependent gene sets refer to table S6. (I) Dot plot showing
   the correlation between the M1 and M2 scores. (J) Heatmap showing
   differences in TF (row) activity in monocytes/macrophages (column). FC,
   fold change.

   Immune checkpoint blockade (ICB) has been proposed as a promising
   treatment for ACP, but its efficacy is interfered by the complex
   immunosuppressive TME in ACP ([164]47, [165]48). Previous studies in
   our laboratory and others have suggested complex expression patterns of
   B7 family ligands/receptors in craniopharyngiomas ([166]7, [167]31,
   [168]49, [169]50). To further elucidate the functional characteristics
   of tumor-infiltrating T cells, we performed differential expression
   analysis of T[reg] and cytotoxic T cells derived from neoplastic and
   nonneoplastic tissues ([170]Fig. 4D). For the two T[reg] clusters,
   CD4-C5-CTLA4 exhibited higher expression of chemokine receptors (such
   as CXCR6) and tissue-resident signatures (such as CD69 and ITGAE),
   which may represent a common feature of tumor-infiltrating T[regs].
   Furthermore, our analysis confirmed that tumor-derived T[regs]
   expressed stronger immune checkpoints, such as CTLA4, TIGIT, TNFRSF4,
   TNFRSF9, HAVCR2, LAYN, and PDCD1. These molecules are widely involved
   in tumor immunosuppression, suggesting a shift in the functional status
   of ACP-infiltrating T[regs] toward immunosuppression. We used a gene
   set of tumor-induced antigen-experienced T[reg] markers to assess the
   immune function status of T[regs] (table S5) ([171]51). ACP-derived
   T[regs] had substantially higher scores, which further supports that
   their activation led to functional tumor T[regs] (fig. S7G). The GSVA
   results suggested that ACP-derived T[regs] showed notable differences
   in multiple biological pathways associated with tumor progression, such
   as “blood vessel morphogenesis” and “transforming growth factor beta1
   production” (fig. S7H). Meanwhile, the activation of macrophage
   chemotaxis and stress-related pathways reflected the response and
   regulation of ACP-derived T[regs] in the tumor inflammatory
   microenvironment. In addition, tumor-derived cytotoxic T cells
   expressed IFNG and the chemokines CCL3, CCL4, and CCL5, suggesting that
   they might be activated ([172]Fig. 4E). In short, our study revealed
   the immunosuppressive status of tumor-infiltrating T cells in the ACP
   TME and their potential role in coordinating ACP progression.

Heterogeneity and state analysis of myeloid subpopulations in ACP

   We reclustered 4205 myeloid cells consisting of 2291 tumor-derived
   cells and 1914 nonneoplastic tissue–derived cells. Myeloid cells showed
   a notable degree of heterogeneity and were divided into 10 clusters,
   including monocytes/macrophages, dendritic cells, and neutrophils,
   based on the expression of known classic markers ([173]Fig. 4, F and G,
   and fig. S8, A and B). Most macrophages originated from tumor tissues,
   suggesting the influence of the ACP TME on the differentiation of
   macrophages. Accumulating evidence suggested that chronic inflammatory
   processes play a key role in ACP tumorigenesis and development, while
   the genesis of cystic components was closely associated with poor
   prognosis ([174]47, [175]52). To further define the characteristics of
   tumor-infiltrating monocytes/macrophages, we scored all
   monocytes/macrophages according to the classic characteristic genes of
   M1/M2 (table S6). We found that the tumor-derived monocyte/macrophage
   subpopulations had higher M1 and M2 scores ([176]Fig. 4H). There was a
   significant positive correlation of feature scores of M1 and M2 between
   all monocytes/macrophages (R = 0.67, P < 0.0001) ([177]Fig. 4I),
   indicating that activation of macrophages in the ACP TME may adopt the
   M1-M2 coupling activation mode. The results of differential gene
   expression analysis showed that tumor-derived monocytes/macrophages
   highly expressed polarization markers (such as CD163) and secreted
   abundant chemokines (such as CXCL8, CXCL16, CCL2, CCL3, and CCL4),
   which may partially explain the recruitment of tumor-infiltrating
   immune cells and the secretion of inflammatory cyst fluid in ACP (fig.
   S8C). The clustering analysis via GSVA suggested that the inflammatory
   response, inflammatory factor release, cell chemotaxis, and macrophage
   differentiation were substantially up-regulated in the
   tumor-infiltrating monocytes/macrophages (fig. S8D). SCENIC analysis
   identified a set of regulation patterns of tumor-specific TFs in
   monocytes/macrophages, such as MAF, NFIL3, CEBPD, CREM, etc.
   ([178]Fig. 4J and fig. S8E). Previous studies have shown that MAF and
   CEBPD are essential for the regulation of macrophage differentiation
   and inflammatory signaling ([179]53, [180]54). The results of the
   SCENIC analysis suggest that ACP-derived macrophages undergo functional
   activation and are involved in the regulation of local inflammatory
   responses in the TME, which may contribute to the up-regulation of
   pro-inflammatory mediator levels at the hypothalamus and pituitary
   during ACP invasion.

Deciphering the molecular phenotype and functional characteristics of glial
cells in the glial scar

   The process of ACP invading brain tissue induces gliosis outside the
   tumor and the formation of glial scar, which represents a unique TME
   for ACP ([181]10, [182]19). To identify the cellular heterogeneity and
   transcriptomic characteristics of glial scars, we extracted
   tumor-derived glial cells and integrated scRNA-seq data derived from
   eight cases of brain tissues for reclustering, which led to a
   transcriptome of 5802 glial cells. On the basis of the differential
   gene expression levels of each cluster and classical markers, we
   identified five cell types, including astrocytes, oligodendrocytes,
   OPCs, neural progenitor cells, and neurons ([183]Fig. 5A and fig. S9, A
   and B). We found that tumor tissue–derived astrocytes expressed higher
   levels of glial fibrillary acidic protein (GFAP), one of the major
   intermediate filament proteins of mature astrocytes, relative to brain
   tissue-derived astrocytes. The antennae of astrocytes derived from
   tumor tissue were thicker, and the GFAP^+ glial fibers of adjacent
   astrocytes overlapped, crossed, and entangled each other to form glial
   scars (fig. S9, C and D). Previous studies have suggested that this
   phenomenon is characteristic of activated astrocytes ([184]55,
   [185]56). Compared with nonneoplastic tissues, tumor-derived astrocytes
   showed a notable enrichment of pathways such as the inflammatory
   response, interferon-γ (IFN-γ) response, and astrocyte activation (fig.
   S9E). In addition, the abnormal activation of multiple metabolic
   pathways was observed in tumor-derived astrocytes, such as arachidonic
   acid metabolism and mucin-type O-glycan biosynthesis fatty acid
   biosynthesis, suggesting a metabolic adaptation of the TME regulated by
   ACP tumors (fig. S9F).

Fig. 5. Detailed characterization of oligodendrocyte lineage in ACP.

   [186]Fig. 5.
   [187]Open in a new tab

   (A) UMAP plot showing the distribution of glial cell clusters after
   reclustering. (B) Violin plots showing the expression levels of
   functional markers in oligodendrocyte lineages across different
   clusters. (C) Multi-IHC staining demonstrates the histological
   distribution of oligodendrocytes with immunological features.
   Arrowheads indicate iOLs. Scale bars, 50 μm. (D) Differential gene
   expression profiles of tumor- and brain-derived OPCs. Genes related to
   antigen presentation were highly expressed in tumor-derived OPCs (red).
   Representative genes are labeled. (E) Heatmap showing the activation
   levels of differential pathways for GO enrichment analysis in OPCs
   across different clusters. (F) Heatmap showing differential TF activity
   across different clusters in the oligodendrocyte lineage revealed by
   SCENIC.

   OLs are widely distributed along the outer side of tumor tissue and are
   involved in forming glial scars and remodeling the TME. To further
   explore its functional characteristics, we applied NMF to all
   oligodendrocytes and OPCs, which resulted in 10 expression signatures
   (fig. S10A and table S7). Among these signatures, OL-3 and OL-4 were
   specifically enriched in tumor tissue–derived OLs, which were composed
   mainly of genes related to oligodendrocyte development and central
   nervous system myelination, such as MAG, MOBP, and FGFR2 (fig. S10A),
   suggesting that tumor-derived OLs may be involved in neural repair
   during tumor invasion of brain tissue. Unexpectedly, we found that
   OL-9, a tumor-specific signature, is composed of immunoregulatory
   genes, suggesting that there may be a distinct immunophenotypic
   oligodendrocyte lineage in the ACP TME. Antigen presentation–related
   genes including major histocompatibility complex II (MHC-II)–encoding
   genes (such as HLA-DRA, HLA-DRB1, and CD74) and immune-type proteasome
   constituent molecules (such as PSMB8 and PSMB9) were enriched in OL-9.
   Previous studies have reported this immunophenotypic OLs (iOLs) in
   models of central system autoimmune diseases such as multiple sclerosis
   (MS) and experimental autoimmune encephalomyelitis (EAE)
   ([188]57–[189]60). Subsequently, we examined the gene expression levels
   of OLs in each cluster according to the classical features described in
   previous studies ([190]Fig. 5B). Compared to nonneoplastic tissue,
   neoplastic tissue–derived OLs had higher expression levels of MHC-I and
   MHC-II constituent molecules, including B2M, HLA-A, HLA-B, HLA-C,
   HLA-DRA, and HLA-DRB1, and genes related to antigen processing and
   presentation, including PSMB8, PSMB9, TAPBP, and TAP1, which suggested
   that they conformed to the classical features of immune-characterized
   OLs described in previous studies. Subsequently, we localized OLs with
   immunological properties in the TME by performing multi-IHC staining of
   MHC-II molecules (HLA-DR and CD74) and the OLs marker OLIG2
   ([191]Fig. 5C). CD68 was used as a negative marker to exclude myeloid
   cells. Furthermore, we compared the expression profiles between
   neoplastic and nonneoplastic tissue derived OPCs ([192]Fig. 5, D and
   E). The integrated analyses of differentially expressed genes and
   pathway enrichment revealed functional characteristics of tumor-derived
   OPCs. Except for MHC-II molecules and antigen-presenting genes,
   tumor-derived OPCs substantially expressed genes related to the IFN
   response (IFITM3 and IFI16), endoplasmic reticulum stress (DNAJB1 and
   HSPA1B), and DNA damage (CDKN2A and BTG2). Previous reports have
   suggested that the abnormal immune inflammatory response may be an
   inducing condition that promotes the functional phenotypic
   transformation of OLs and the expression of MHC-II molecules,
   especially the stimulation of IFN-γ ([193]58). The results of
   enzyme-linked immunosorbent assay (ELISA) and Western blot experiments
   indicated that ACP had a higher expression level of IFN-γ relative to
   brain tissue (fig. S10, B and C). In addition, the results of multi-IHC
   staining showed higher IFN-γ expression level in the ACP TME where OLs
   were located relative to the brain tissue environment (fig. S10D). We
   found that CD8^+ T cells were also distributed in a close proximity to
   OLs, which was not observed in brain tissue (fig. S10D). Previous
   studies have suggested that CD8^+ T cells are capable of recognizing
   iOLs and exerting cytotoxic effects through antigen presentation, which
   may contribute to neural damage upon ACP tumor progression ([194]59).
   Subsequently, SCENIC analysis revealed a potential regulatory network
   of TFs across subtypes in the OLs ([195]Fig. 5F). For example, IRF1 and
   STAT3 were identified as candidate TFs underlying the gene expression
   differences in iOLs. IRF1 and STAT3 are reported to be up-regulated
   after stimulation by immune and inflammatory factors (especially IFN-γ)
   and are involved in various downstream pathways, such as oxidative
   stress and apoptosis, as response regulators. To determine whether iOLs
   are also present in other brain tumors, we analyzed the glioblastoma
   (GBM) scRNA-seq data from public dataset ([196]61). After clustering
   and cell annotation, we calculated the immune score of each cluster
   based on the gene set of antigen presentation process (including
   HLA-DRA, HLA-DRB1, HLA-DRB5, HLA-DPA1, HLA-DPB1, CD74, PSMB8, PSMB9,
   TAP1, and TAPBP; fig. S11, A to C). However, we did not observe OL
   clusters with similar transcriptional characteristics to iOLs in GBM
   (fig. S11C). As an immunologically cold tumor, the immune
   microenvironment of GBM may not support the induction and
   transformation of iOLs.

   Then, we focused on transcriptomic differences in iOLs between tumor
   and autoimmune disease models. Differential gene expression analysis
   revealed that cell adhesion–related extracellular matrix components,
   such as proteoglycans (including VCAN, CSPG4, SDC3, LGALS1, and LGALS3)
   and VIM, were up-regulated in tumor-derived iOLs (fig. S10E). We found
   that OLs with high expression levels of VCAN were widely distributed
   outside the tumor tissue by immunostaining (fig. S10F). Previous
   studies support that VCAN is one of the key components of the
   extracellular matrix in various malignant tumors and contributes to
   tumor development and the formation of the preinvasive
   microenvironment, suggesting that OLs may be involved in ACP TME
   remodeling as potential regulators ([197]62, [198]63). We subsequently
   analyzed RNA sequencing data derived from MS samples in public datasets
   and did not find similar signatures, suggesting that the influence of
   the TME may induce the transformation of iOLs into a distinct
   tumor-promoting phenotype (fig. S10, E and F).

Establish an intercellular communication network in the ACP TME

   Next, we performed CellChat analysis to identify cell-cell interactions
   in the ACP TME. Potential ligand-receptor pair interactions between
   major cell types in the TME revealed a complex cellular communication
   network (fig. S12A). First, we focused on communication events that
   occurred within the tumor. WE_2 mediated the most efferent and afferent
   events from other subpopulations in the ACP, possibly representing a
   signaling center and secretory hub within the tumor ([199]Fig. 6A and
   fig. S12B). Considering that ACP is a WNT signaling–driven tumor, we
   focused on the distribution of WNT-related communication events across
   subpopulations. Our results indicated that WE mediated the most
   frequent WNT-related signaling among the four major neoplastic cell
   states, and the target cells receiving WNT signaling in the ACP TME
   consisted mainly of KE_3, WE_2, KE_1, and KE_2 ([200]Fig. 6B). Among
   all neoplastic cell states, WE expressed the highest levels of WNT
   pathway ligands, including WNT5A, WNT6, WNT10A, and WNT10B, as well as
   ligand genes for its downstream related pathways, such as fibroblast
   growth factor (FGF) (FGF4 and FGF20), bone morphogenetic protein (BMP)
   (BMP4 and BMP7), and Hedgehog (SHH) signaling pathways ([201]Fig. 6C).
   Their receptor distribution covered almost all neoplastic cells,
   suggesting that WE plays a crucial role in abnormal activation of WNT
   signaling and tumor development ([202]Fig. 6C). In addition, there were
   some typical ECM-receptor interactions between neoplastic cells, such
   as LAMB1-CD44, POSTN-ITGAV, and COL6A1-ITGB1, which might contribute to
   tumor proliferation and invasion during ACP progression.

Fig. 6. Complex cell-cell communication networks in the ACP TME.

   [203]Fig. 6.
   [204]Open in a new tab

   (A) Heatmap showing the strength of incoming and outgoing events in
   interactions between different clusters in the ACP. Histograms
   separately count the overall intensity of outgoing (y axis) and
   incoming (x axis) events for each cluster. (B) Interaction strength of
   WNT signaling pathway incoming and outgoing events between all clusters
   in ACP. (C) Dot plots show gene expression levels of receptor-ligand
   pairs involved in interactions between different clusters in ACP.

   We next focused on the interaction between neoplastic cells and stromal
   cells in the TME (fig. S12C). Neoplastic cells secreted multiple
   cytokines, including FGF2, FGF9, and PDGFA, the receptors of which were
   prominently expressed in the oligodendrocyte lineages. Previous studies
   have demonstrated that in models of brain injury or inflammation,
   abundantly secreted cytokines such as FGF2 and PDGFA recruit OLs to
   migrate to lesion sites to participate in neural repair
   ([205]64–[206]66). Furthermore, we observed that neurotrophic factors
   such as BDNF, NTF4, IGF1, and NRG2 were enriched in the TME of ACP.
   Their target receptors, such as IGF1R, ERBB3, ERBB4, NTRK2, and LIFR,
   were distributed in the oligodendrocyte lineage, which may contribute
   to the survival and proliferation of OL cells in an inflammatory
   environment.

   Infiltration of immune cells is a typical feature of ACP ([207]67).
   Abundant expression of immunochemokines (CXCL1, CXCL2, CXCL6, CXCL8,
   CXCL16, CCL3, and CCL4) was noted in the ACP TME, while the
   corresponding receptors were widely expressed in T cells and
   macrophages (fig. S12C). This chemokine-mediated cellular communication
   may play critical roles in maintaining inflammatory cell infiltration
   and immune microenvironment homeostasis in ACP. Furthermore, we
   speculated that some ligand-receptor pairs (such as MIF-CD74/CD44,
   TGFB1-TGFBR1, IL10-IL10RA/IL10RB, ADGRE5-CD55, and PROS1-AXL/MERTK) in
   the ACP TME may be involved in the transition of the
   monocyte/macrophage functional state to anti-inflammatory activation,
   and these secretory signals were secreted mainly by neoplastic cells
   (fig. S12C). In previous studies, secretory signals were involved in
   the recruitment and chemotaxis of tumor-associated macrophages as
   pleiotropic immune regulators. In addition, migration inhibition
   factor, transforming growth factor–β1 (TGF-β1), and interleukin-10
   (IL10) have been shown to be strong inducers of M2 polarization in
   tumor-associated macrophages ([208]68, [209]69). Pathway enrichment
   analysis of genes related to tumor monocyte/macrophage cell
   communication also suggested an induction of tumor paracrine signaling
   on macrophage chemotaxis and phenotypic differentiation ([210]Fig. 7A).
   Multi-IHC revealed that monocytes and macrophages were dominantly
   located within ACP tissues and exhibited M2 polarization, which would
   promote an immunosuppressive microenvironment to facilitate tumor
   progression ([211]Fig. 7B). To further explore the effect of tumors on
   the functional status of macrophages, we cocultured ACP neoplastic
   cells with PBMC-derived monocytes/macrophages and THP-1 (a human
   monocytic cell line) ([212]70). We dissociated ACP tissue and collected
   neoplastic cells for in vitro culture ([213]Fig. 7C). The results of
   Western blot showed that monocytes cocultured with neoplastic cell
   culture medium (CM) exhibited characteristics of macrophage
   differentiation, for example, the expression of macrophage
   differentiation markers (CD11c, CD11b, and LOX1) and scavenger receptor
   (CD204) ([214]Fig. 7D). We then focused on the effect of tumor
   coculture on macrophages polarization. PBMCs and THP-1 were stimulated
   with phorbol 12-myristate 13-acetate (PMA; 100 μM) for 24 hours to be
   induced into M0 macrophages. First, the quantitative polymerase chain
   reaction (qPCR) results indicated that cocultured macrophages expressed
   significantly higher levels of M2 polarization markers (IL10, TGFB1,
   CD163, CD206, and ARG1), relative to M1 polarization markers (TNF, IL6,
   and CD86) ([215]Fig. 7E). Similar results were obtained when detecting
   the expression of secreted proteins and cell membrane surface markers
   by ELISA and flow cytometry, respectively ([216]Fig. 7, F and G). Last,
   the transwell assay revealed that the neoplastic cells cultured in the
   lower chamber promoted the migration of macrophages in the upper
   chamber, suggesting that the paracrine effect of neoplastic cells may
   be involved in the regulatory mechanism of macrophage chemotaxis
   ([217]Fig. 7H).

Fig. 7. Coculture with ACP neoplastic cells induced M2 polarization of
macrophages.

   [218]Fig. 7.
   [219]Open in a new tab

   (A) GO enrichment analysis of genes involved in tumor
   monocyte/macrophage intercellular communication. The filled color
   indicates the number of genes involved in the pathways, and the length
   of the column indicates the P value. (B) Multi-IHC detection of the
   distribution of M2 phenotype macrophages in the TME. Scale bars, 50 μm.
   (C) Identification of dissociated primary ACP neoplastic cells by
   immunofluorescence. According to previous sequencing results, CDH1 was
   used as the positive marker, and GFAP was used as the negative marker.
   Scale bars, 20 μm. (D) After coculture of THP-1 and PBMC with ACP
   neoplastic cells, the protein expression of macrophage differentiation
   marker (CD11c, CD11b, CD204, and LOX1) was detected by Western blot.
   (E) Detection of mRNA levels of M1 (TNF, IL6, and CD86) and M2 (IL10,
   TGFB1, CD163, CD206, and ARG1) polarization marker genes in THP-1 and
   PBMC-derived macrophages by qPCR assay. (F) Detection of cytokine
   expression levels in cell supernatants in THP-1 and PBMC-derived
   macrophages by ELISA. (G) The expression levels of cell surface
   polarization marker proteins in THP-1– and PBMC-derived macrophages
   were detected by flow cytometry. (H) Transwell assay was used to detect
   the ability of macrophages in the upper chamber to migrate to the lower
   chamber. Scale bars, 50 μm. *P < 0.05, **P < 0.01, and ***P < 0.001.
   ns, not significant.

DISCUSSION

   Clinically, the outcome of existing treatment for ACP remains
   unsatisfactory. The tumor heterogeneity of ACP and the complex TME
   surrounding the outer side of finger-like protrusions have led to
   therapeutic resistance and tumor recurrence. Hence, a deep
   understanding of the cell composition and functional characteristics of
   ACP is urgently needed. In this study, we undertook single-cell
   transcriptomic analysis and established an atlas of cell maps tailored
   for ACP tumor tissues, which revealed the characteristics of neoplastic
   cells and the TME. Our in-depth analysis not only uncovered the
   phenotypic differences between the TME and normal tissue-derived
   stromal cells underlying tumor progression but also reconstructed
   complex cellular communication and molecular networks within the TME.

   To date, the molecular characteristics, cell function, dynamics, and
   plasticity of ACP neoplastic cells are still unclear. The solid part of
   ACP is composed mainly of PE, stellate cells, and whorl-like cell
   clusters, which reflects the complicated composition of ACP neoplastic
   cells. In this study, we identified 10 subpopulations of neoplastic
   cells, which were further redefined as four main cell states based on
   transcriptomic characteristics across tumors, including PE, WE, KE, and
   CC. Different from the growth mode of neoplastic cells infiltrating
   each other in most solid tumors, these main neoplastic cell states of
   ACP were characterized by independent histological localization and
   notable differences in cell morphology and polarity, which suggests
   that they had different molecular characteristics and cellular
   kinetics. PE was distributed mainly along the outermost side of
   finger-like protrusions and was immediately adjacent to reactive glial
   scars. Neoplastic cells in PE highly expressed genes involved in the
   early response to stress stimuli, which may be related to their
   response to the extensive inflammatory environment in glial scar
   tissues. In addition, two other neoplastic cell states, KE and WE, were
   distributed mainly medial to the PE in the tumor tissue area. KE had a
   specific expression profile of keratins, such as cytokeratin KRT6,
   KRT7, and KRT13 and hair follicle keratin KRT75, suggesting that they
   represented a more mature and unique neoplastic cell population at the
   stage of epithelial development. WE was distributed in the whorl-like
   structures inside the finger-like protrusions, which was characterized
   mainly by persistent activation of the WNT signaling pathway.
   Furthermore, we observed that WE secreted multiple cytokines, such as
   FGF4, FGF8, BMP2, BMP4, and SHH, which were regulated by the WNT
   signaling and were previously reported to contribute to the
   tumorigenesis and development of ACP. These indications suggest that WE
   tended to act as a functional core and were critical in ACP tumor
   progression, corroborating the views of previous studies ([220]19,
   [221]42). On the basis of tissue structure and gene expression
   profiling, early studies suggested that tumorigenesis and progression
   of ACP may be related to the development of ectopic ectoderm ([222]71,
   [223]72). It is worth mentioning that although WE and KE are widely
   distributed in ACP tissues and have distinct biological phenotypes,
   both of them exhibit lower proliferative activity compared to PE.
   Meanwhile, transcriptome splicing kinetics and pseudotime analysis
   revealed distinct evolutionary trajectories of WE and KE in ACP
   neoplastic cells. Consistent with previous study, only WE among the
   four cell states showed nuclear accumulation of β-catenin and
   activation of transcriptional regulatory networks downstream of the WNT
   pathway, although the vast majority of neoplastic cells harbored
   activating mutations of CTNNB1 that drive ACP tumorigenesis ([224]22).
   In general, CTNNB1 exon 3 mutations promote nuclear accumulation of
   β-catenin by preventing its phosphorylation-dependent degradation
   ([225]73). However, the widespread distribution of membrane-localized
   β-catenin implied that CTNNB1 mutation did not appear to be a
   sufficient condition for the nuclear accumulation of β-catenin and WNT
   signaling activation in ACP neoplastic cells. Relative to other
   neoplastic cells, WE present a cellular state that seems like more
   favorable for nuclear translocation and nuclear homeostasis of
   β-catenin. We observed that MUC1, a novel biomarker of WE, exhibited
   cell-wise colocalization with nuclear accumulation of β-catenin.
   Previous studies have shown that the cytoplasmic tail of MUC1 is
   capable of binding to β-catenin and promoting steady-state levels of
   nuclear β-catenin ([226]74, [227]75). Furthermore, MUC1 bound to
   β-catenin would inhibit GSK3β-mediated phosphorylation and downstream
   degradation of β-catenin by acting as competitive substrates ([228]76).
   Because of off-target effects, inhibitors targeting β-catenin mutations
   have not been incorporated into the current clinical treatment options
   for ACP ([229]7). However, WE is still the most promising clinical
   therapeutic target for ACP. Considering the regulatory role of WE as a
   paracrine hub for multiple cytokines on the development of other
   neoplastic cell states, targeting WE may interfere with the signaling
   cross-talk within ACP and inhibit tumor progression. In addition,
   because of the active proliferative capacity, PE may be another
   potential clinical therapeutic target for precision medicine in ACP.

   The progression of ACP is accompanied by the formation of glial scar,
   which causes local immune cell infiltration. In terms of T cell
   function in the TME, CD4^+ and CD8^+ T cells showed strong expression
   of tissue resident markers such as CD69 and ITGAE as well as varying
   degrees of immune depletion. Tumor-infiltrating T[regs] widely express
   immune checkpoints such as CTLA4, TIGIT, and PDCD1, suggesting that
   they may mediate an immunosuppressive microenvironment. Limited by
   inflammatory or postoperative adhesions, intracystic IFN-α therapy has
   been used in the management of patients with recurrent ACP. Notably,
   ICB targeting PD-1 and CTLA-4 has been proposed to produce synergistic
   effects with IFN-α therapy ([230]77). Combining ICB therapy plus IFN
   therapy is expected to be a promising nonsurgical treatment for the
   clinical management of patients with ACP. Macrophages represent the
   most abundant immune cell type in the TME of ACP and are involved in
   the regulation of inflammation levels. We found that unlike the classic
   model of M1-M2 polarization, ACP-infiltrated macrophages exhibited
   M1-M2 coupled activation, which was reported in previous studies of
   nasopharyngeal carcinoma and breast cancer ([231]78, [232]79). This
   activation mode of macrophages may lead to chronic inflammation in the
   glial scar during the growth and invasion of brain tissue by ACP.
   Furthermore, it may be the pathological basis of damage caused by ACP
   to the hypothalamus and pituitary function, which should be considered
   for tumor treatment and brain function maintenance.

   Reactive glial scars surround the outside of ACP neoplastic cells,
   constituting their unique tissue microenvironment. They may play a role
   in tumor progression, which was not elucidated in previous studies. One
   of the characteristics of glial scars is the activation of astrocytes.
   Affected by local chronic inflammation and tumor invasion, astrocytes
   exhibited larger GFAP-positive glial fibers than those in normal brain
   tissues. Their cell morphology and polarity are changed and arranged in
   a staggered manner. When the tumor invades the brain, the glial scars
   act as a buffer area in the organization process. Although we did not
   detect neurons, this area is enriched by a large number of OLs,
   possibly the result of migration of OLs induced by the paracrine
   effects of the tumor, suggesting that the infiltrating OLs in the TME
   may be involved in promoting tumor progression. It is worth noting that
   we discovered a cluster of OLs with immune characteristics (iOLs),
   which is manifested mainly by the expression of immune cell markers
   related to antigen presentation, such as HLA-DR and CD74. MHC-II
   components are usually expressed in professional antigen-presenting
   cells such as dendritic cells and macrophages. Previous studies
   reported this glial-immune dual characteristic OLs and its
   characteristics in autoimmune disease models such as MS and EAE, which
   has a molecular phenotype similar to iOLs in the ACP TME ([233]57,
   [234]58). In addition, ACP-derived iOLs have more unique
   characteristics, which are possibly influenced by tumor paracrine
   signaling. For example, abundant expression of proteoglycan-dominated
   extracellular matrix genes may contribute to tumor migration and
   invasion, indicative of a unique functional transformation of the
   oligodendrocyte lineage in tumor tissues.

   In summary, our research constructed a cellular atlas of ACP neoplastic
   cells and the tissue microenvironment and described their
   characteristics and interconnected molecular networks. This study will
   not only provide a solid foundation for further studies that aim to
   reveal the cellular and molecular basis of tumorigenesis and
   progression of ACP but also provide valuable resources for potential
   clinical treatment and drug development.

MATERIALS AND METHODS

ACP samples and normal glial tissue

   The study was approved by the Institutional Review Board of West China
   Hospital of Sichuan University (file no. 20211047A). All samples were
   collected at West China Hospital of Sichuan University, and informed
   consent was obtained from all patients. We recruited three patients who
   were pathologically identified as having ACP for single-cell sequencing
   (table S1). None of the patients had received radiotherapy,
   chemotherapy, or other forms of antitumor therapy before surgery. The
   aggressiveness of ACP causes the tumor and glial tissue to adhere
   seriously. Fresh ACP tumor tissue and adjacent glial tissue were
   obtained immediately after surgical resection. Before participating,
   all patients signed an informed consent form, and all research
   processes were conducted under the regulations of the Declaration of
   Helsinki.

Tissue dissociation and single-cell suspension preparation

   Fresh surgical specimens were cut into 1- to 2-mm^3 cubes in a petri
   dish and rinsed with prechilled Dulbecco’s phosphate-buffered saline
   (DPBS) three times to remove blood clots. The tissue block was
   incubated with a digestion solution made of 0.1% collagenase IV and
   0.1% papain, shaken, and digested in a 37°C constant temperature water
   bath for 30 min. The digestion juice was shaken gently every 3 to
   5 min. When no obvious large pieces of tissue remained, the digestion
   solution was filtered through a 70-μm cell sieve and then through a
   40-μm cell sieve again, and the filtrate was collected into a new
   centrifuge tube. The filtrate was centrifuged at 250g for 5 min, the
   supernatant was discarded, and the cells were washed three times with
   precooled DPBS. Then, the cell pellet was resuspended in 1 ml of
   precooled red blood cell lysate and incubated in a 4°C refrigerator
   with shaking for 15 min. Then, 10 ml of precooled DPBS was added and
   centrifuged at 250g for 5 min. The cell pellet was resuspended in 3 ml
   of DPBS without calcium and magnesium ions. The cells were counted with
   the Countes II Automated Cell Counter. Ten-microliter cell suspension
   was stained with trypan blue to evaluate cell viability.

Droplet-based scRNA-seq library preparation

   According to the manufacturer’s instructions, the Chromium Single Cell
   3′ Library and Gel Bead and Multiplex Kit (10X Genomics, V3) were used
   to encapsulate single cells into droplets to construct scRNA-seq
   libraries. In short, the single-cell suspension with viability higher
   than 90% was loaded into the Chromium Single-Cell Controller Instrument
   (10X Genomics) to generate single-cell gel beads in emulsions.
   Approximately 15,000 cells were added to each channel, aiming to
   capture 5000 to 10,000 cells per sample. Next, reverse transcription
   (RT) cleanup, cDNA amplification, and library quality assessment were
   performed according to the instructions of the kit. The scRNA-seq
   library was sequenced on the Illumina NovaSeq 6000 platform to generate
   150-bp paired-end reads.

Preprocessing of scRNA-seq data

   CellRanger (10X Genomics, 3.0.1) was used to map the sequencing data to
   the human reference genome (hg38). The scCancer package (2.2.1) based
   on R software (4.0.5) was used to perform quality control on sequencing
   data and filter low-quality cells based on the ratio of
   mitochondrial/ribosomal genes and dissociation-associated genes. The
   expression matrices of all samples were converted into Seurat objects
   through the Seurat package (4.0.0). To eliminate the potential impact
   of double cells on the results, we used the DoubletFinder package
   (2.0.3) and integrated the factors of gene number and unique molecular
   identifier (UMI) number to filter the cells. Last, we retained 27,696
   high-quality single cells for downstream analysis.

Integration, dimensionality reduction, clustering, and cell annotation

   To eliminate the batch effect, different samples were integrated
   through the default parameters of the “FindIntegrationAnchors” and
   “IntegrateData” functions in the Seurat package. Subsequently, to
   better reflect the biological functional characteristics of the cells,
   we used the “NormalizeData” function to normalize the expression matrix
   and the “FindVariableFeatures” function to determine the set of
   intercellular variant genes in the dataset as highly variable genes
   (HVGs). To perform dimensionality reduction and clustering on the
   expression matrix, we used “ScaleData” to scale the data, and “RunPCA”
   was used to analyze the first 50 main components of the expression
   matrix. Subsequently, “FindNeighbors” and “FindClusters” were used to
   cluster the first 25 main components and project them to
   two-dimensional Uniform Manifold Approximation and Projection (UMAP) or
   t-distributed stochastic neighbor embedding images for display. The
   characteristic genes of each cell subgroup were calculated through the
   built-in “FindAllMarkers” function of the Seurat package, and then the
   cell type of each subgroup was defined on the basis of the SingleR
   package and the recognized cell marker genes in previous studies.

CNV analysis and neoplastic cells identification

   To identify neoplastic and nonneoplastic cells, we used the inferCNV
   package (1.6.0) based on R software to calculate the CNV pattern in
   single cells. The default parameters and clusters composed of immune
   cells were used as a reference to infer real neoplastic cells. As a
   verification of the inference results, the expression levels of ACP
   tumor marker genes obtained based on bulk RNA-seq in previous studies
   were detected in each cell cluster. The clonality tree was generated
   using the UPhyloplot2 algorithm (v2.3) to analyze subclonal patterns in
   ACP ([235]https://github.com/harbourlab/UPhyloplot2).

NMF analysis of the sample heterogeneity expression signatures

   For neoplastic cells and oligodendrocyte lineages in ACP, NMF was
   performed to identify variable expression signatures based on the NMF
   package (0.23.0). For each sample, we reserved the first 2000 HVGs for
   the analysis of the expression program in the sample. All negative
   values in the expression matrix were replaced with 0. For neoplastic
   cells, we lastly extracted 25 potential expression signatures from
   three ACP tumors (table S3). For the oligodendrocyte lineages, 10
   expression signatures were extracted (table S7). To study the common
   characteristics of ACP across samples, we used the R AUCell package to
   calculate the score of each program in each tumor and then performed
   hierarchical clustering of all signatures and calculated the Pearson
   index. In the end, four meta-signatures across tumors were identified,
   representing the cell cycle, WNT signaling pathway, epithelial
   differentiation, and stress. For the determination of neoplastic cell
   states, we used the built-in “AddModuleScore” function of the Seurat
   package to score the function of single cells. The scoring gene set
   comes from the four meta-signature marker genes (table S4).

Gene set variation analysis

   The GSVA package (1.38.2) was used to perform pathway enrichment
   analysis based on scRNA-seq data. The pathways used in the enrichment
   analysis [including gene ontology (GO), HallMark, and Kyoto
   Encyclopedia of Genes and Genomes pathways] and the corresponding gene
   sets were derived mainly from the MSigDB database.

TF regulatory network analysis

   The SCENIC package (1.2.4) was used to construct a TF regulatory
   network in scRNA-seq data and identify the functional status of cells.
   To examine differences in the activity of TFs and their downstream
   target genes across cell clusters, we performed SCENIC analysis on all
   single cells. The quantification matrix of TF activity was imported
   into the Seurat object, and subgroup-specific TFs were found via the
   “FindAllMarkers” function.

RNA velocity and pseudotime analysis

   scVelo (0.2.4) was used to perform RNA velocity analysis to explore the
   potential evolutionary relationship between cell subpopulations. All
   parameters were set to default values, and the RNA velocity results
   were mapped to UMAP images. On the basis of the different trajectories
   displayed by the RNA velocity results, the corresponding cell
   subpopulations were extracted, and Monocle3 pseudotime sequence
   analysis was performed separately. The analysis results of Monocle3,
   including the pseudotime score and trajectory, were projected onto the
   UMAP graph of redimensionality reduction and clustering. The
   “plot_genes_in_pseudotime” function was used to display the genes with
   notable differences in the pseudotime trajectory.

Immune cell function status score

   We summarized the list of characteristic genes that define T cell
   functions (such as cytotoxicity, depletion, naïve cells, T[regs], and
   tissue retention; table S5) and macrophage functions (such as M1 and M2
   polarization; table S6) in published studies. The scoring of cell
   function based on a specific gene set was performed through the
   built-in “AddModuleScore” function of the Seurat package.

Cell-cell communication network

   We used CellChat (1.0.0), a database that can accurately analyze
   cell-cell signal connections and build communication networks, to
   explore potential cell-cell interactions in the ACP TME. By using the
   default parameters and the ligand-receptor pair database, if cell type
   A expresses the ligand that is paired with the receptor expressed by
   cell type B, then the event is regarded as a valid interaction.
   According to the molecular characteristics of the ligand-receptor pair,
   the interaction events are divided into secretory signals, cell-cell
   contact, and ECM receptors. We paid particular attention to the
   potential communication relationships between tumor cell subgroups,
   tumor immune cells, and tumor glial cells. Communication events between
   different cell types were displayed in the form of connecting lines.

IHC and immunofluorescence

   The tissue samples were immersed in formalin and fixed for 24 to 48
   hours immediately after surgical resection. The tissues were
   subsequently dehydrated, embedded in paraffin, and sectioned according
   to the usual procedures. The paraffin sections were baked in an oven at
   65°C for 2 hours and then dewaxed and hydrated with xylene and gradient
   alcohol. EDTA buffer (pH 9.0) or sodium citrate buffer (pH 6.0) was
   used for antigen retrieval. The sections were preincubated with
   antibody diluent to seal the nonspecific antibody binding sites of the
   tissue and then incubated with the primary antibody and secondary
   antibody sequentially according to the experimental protocol. A
   PerkinElmer Multiplex IHC staining kit was used for multi-IHC staining
   according to the manufacturer’s instructions. After incubating the
   primary antibody, the slide was incubated with the secondary antibody,
   and tyrosine signal amplification was performed. Subsequently, the
   sections were subjected to microwave heat repair, and the next round of
   staining was performed. After all the indicators were stained,
   4′,6-diamidino-2-phenylindole was used to mark the nucleus. Vectra
   Polaris (PerkinElmer) and Qupath (0.3.0) were used for multispectral
   slice scanning and image analysis. All antibody information is included
   in the Supplementary Materials (table S8).

Isolation of PBMCs

   Peripheral blood from patients with ACP was anticoagulated with heparin
   after collection. Peripheral blood was subjected to density gradient
   centrifugation with mononuclear cell separation solution (Haoyang,
   China) to obtain PBMCs. After incubating peripheral blood mononuclear
   cells with FITC-CD14 monoclonal antibody, flow cytometry was performed
   to obtain purified CD14^+ PBMCs.

Cell culture

   Primary ACP cells were cultured in EpiLife medium (Thermo Fisher
   Scientific), and human keratinocyte growth supplement (Thermo Fisher
   Scientific) was added. Half of the medium was changed every 48 hours.
   THP-1 cells and PBMCs were cultured in 1640 medium supplemented with
   10% fetal bovine serum. PMS (100 nM, MCE) was used to treat THP-1 cells
   and PBMCs for 24 hours, and then the medium was changed to induce them
   to differentiate into nonactivated (M0) macrophages.

Tumor monocyte/macrophage coculture in vitro

   To explore the effect of neoplastic cell paracrine signaling on the
   functional status of monocytes/macrophages, we collected the
   supernatant of neoplastic cell culture medium and added it to
   monocyte/macrophage culture dishes. Half of the medium was changed
   every 24 hours.

Flow cytometry

   Macrophages were blocked with Fc receptor blocking reagent (BioLegend,
   #422301) before immunostaining. The cells were incubated with
   monoclonal antibodies labeled with the following fluorescent dyes and
   detected by flow cytometry: CD86-FITC (BioLegend, #374203), CD163-Alexa
   647 (BioLegend, #326507), CD206-APC (BioLegend, #321109), and CD14-FITC
   (BioLegend, #325603).

RNA extraction and real-time quantitative PCR

   The extraction of total mRNA, RT, and qPCR were carried out in
   accordance with the published methods of previous studies. TRIzol
   (Invitrogen) was used to isolate and extract total mRNA from cells or
   tissues, and after determining the concentration, the PrimeScript RT
   reagent Kit (Takara) was used for RT to obtain cDNA. According to the
   manufacturer’s requirements, TB Green Premix Ex Taq II (Takara) was
   used for qPCR analysis. The amplification step in the PCR consisted of
   40 cycles of denaturation (95°C, 15 s) and annealing extension (60°C,
   60 s). After the last cycle, all samples were subjected to melting
   point analysis with a temperature range of 60° to 95°C. Using ACTB as a
   housekeeping gene, the relative expression value of the gene was
   calculated. Refer to table S9 for the primer information required for
   the reaction.

Western blot

   Total protein was collected from cells or tissues with
   radioimmunoprecipitation assay lysis solution (Thermo Fisher
   Scientific) and boiled to denature the total protein. The protein
   concentration was detected using a bicinchoninic acid kit. Protein
   samples (10 to 20 μg) were run on SDS–polyacrylamide gel
   electrophoresis gels for electrophoresis. After transferring the
   protein to the polyvinylidene difluoride membrane, according to the
   experimental requirements, the membrane was blocked and incubated with
   the corresponding primary antibody at 4°C for 12 to 16 hours.
   Subsequently, the membrane was washed three times, after which it was
   incubated with the horseradish peroxidase–conjugated secondary
   antibody, and the color reaction was performed with the enhanced
   chemiluminescence (ECL) kit to detect the Western blot. β-Tubulin was
   used as a housekeeping gene to evaluate the relative expression value
   of each index. Refer to table S8 for all antibody information.

Cytokine detection

   According to the manufacturer’s specifications, an ELISA kit (R&D
   Systems) was used to detect the levels of IL4, TGFβ1, IL12, and tumor
   necrosis factor–α in the supernatant of macrophages. IFN-γ is detected
   in tumors and brain tissues. Standards with different concentration
   gradients were tested to construct a standard curve by the
   four-parameter logistic fitting curve method.

Cell migration experiment

   Transwell chambers (Corning, catalog no. 3422) with 8-μm pores were
   used to detect the migration of macrophages. In short, 1 × 10^4 M0
   macrophages (from THP-1 cells and PBMCs pretreated with PMA for
   24 hours) were seeded in the upper chambers, and 1 × 10^4 ACP primary
   neoplastic cells were placed in the lower chambers. The culture plate
   was placed in a 37°C incubator for 48 hours, and then the transwell
   chambers were removed and fixed with 4% paraformaldehyde. After
   staining the chambers with crystal violet, a moistened cotton swab was
   used to carefully wipe away the remaining cells inside the chambers. A
   microscope was used to observe and count the cells that migrated on the
   bottom surface of the chambers.

Statistical analysis

   Data statistics are performed based on GraphPad Prism 8.0 software. We
   obtained the P value of the difference between different groups by
   performing Student’s t test and one-way analysis of variance (ANOVA).
   When P < 0.05, differences were considered statistically significant.

Acknowledgments