Cellular map of ACPs shows significant inter-tumoral heterogeneity

To generate a cellular map of ACP tumors, we used single-cell RNA (scRNA-seq) and T cell receptor (TCR) sequencing (scTCR-seq) to profile surgically resected fresh tumors from 12 ACP patients, among which 7 are primary tumors and the other 5 are relapsed tumors (Suppl. Table 1). The ages at the initial diagnosis range from 3 to 50 years old. To localize the spatial distribution of different types of the cells in the tumor microenvironment, we also used the Visium Spatial Transcriptomics technology to profile 3 formalin-fixed paraffin embedded (FFPE) tumor tissue sections out of the 12 ACPs.

After rigorous quality control, the transcriptomic profiles of a total of 70,682 cells were obtained. We were able to identity 15 major cell types based on their expression of known marker genes (Fig. 1 and Suppl. Table 2). These cell types can be grouped into six main categories (1) tumor cells; (2) myeloid immune cells (including various dendritic cells [DCs], conventional type 1 DC [cDC1], myeloid DCs [mDCs] and plasmacytoid DCs [pDCs]); (3) lymphoid immune cells (T cells, NK cells, mast, plasma and B cells); (4) pituitary cells; (5) neural cells; and (6) endothelial cells and other stromal cells.

Fig. 1

Single-cell transcriptomic profiling of adamantinomatous craniopharyngiomas revealing intra- and inter-tumor heterogeneity. a The UMAP of scRNA-seq data from a total of 12 ACPs. Cell clusters are denoted by color, and labeled with inferred cell types (top). The UMAP are also colored by the mutational status of the CTNNB1 exon 3 in each cell (bottom left, gray dots for the cells with no reads covering the mutation loci, red dots for the cells with at least one read carrying the mutation bases [ALT], i.e., possibly mutated cells, and cyan dots for the cells with no alternative base found in the aligned reads [REF], i.e., likely the wildtype cells), and by the inferred cell cycling phase (bottom right). The circled clusters are the cycling cells (in G2M/S phases) with strong proliferating capability and they are found in the tumor cells, immune myeloid cells and T cells, and stromal cells. b Heatmap showing the expression of marker genes in 500 randomly sampled cells per cluster and a few representative gene labels are highlighted (see Suppl. Table 2 for the gene list). c Proportions of tumor cells, B cells and stromal cells among 12 ACPs, which reveals that there exist three 3 subgroups of ACPs. UMAP, Uniform Manifold Approximation and Projection; scRNA-seq, single-cell RNA sequencing; ACPs, adamantinomatous craniopharyngiomas; NK, natural killer; DCs, dendritic cells; cDC1, conventional type 1 DC; mDCs myeloid DCs; pDCs, plasmacytoid DCs

Whole-exome sequencing (WES) of the tumors and the matched blood samples was applied to determine the somatic variants in each tumor. Bulk RNA sequencing profiles were also available for 8 out of the 12 tumors. All the ACPs carried a driver mutation in CTNNB1 exon 3, called from either WES or bulk RNA-seq data. The mutations in a majority of the tumors (7/12) occur at the late phosphorylation codon (4 at Ser33) or the flanking codons (2 at Asp32 and 1 at Gly34), 3 tumors at Thr41, and 2 tumors at Ser37. The variant allele fraction (VAF) ranges from 0 to 0.424, estimated from the WES reads, indicating a diverse range of tumor cell purity. These variants were used to annotate the cellular mutational status of CTNNB1 exon 3 in the scRNA-seq data. Using inhouse code, each cell was annotated either as wildtype (REF, all the reads in the cell carrying the reference bases), mutated (ALT, at least one read carrying the alternative bases) or unknown (no reads covering the variant loci). With this information, two clusters of tumor cells were unequivocally identified. The cells expressing the marker genes of whorl-like cluster (WC) cells obtained from RNA-sequencing of laser-capture microdissection specimens, such as inhibitors of WNT/β-catenin signaling pathway, DKK4 and NOTUM, forms an isolated cluster, showing distinct expression patterns from the other tumor cells.

As it has been observed in bulk RNA profiles and DNA methylome profiles, significant intertumoral heterogeneity was observed in terms of distribution of different cell populations. Four specimens (P455, P432, P452, P420) have a high proportion of tumor cells (> 8% in the corresponding specimens), while B cells are mainly found in P431, P419, and P457 tumors (> 14% in the corresponding specimens) in which the stromal cells are mostly missing (< 5% in the corresponding specimens; Fig. 1c).

The cell-cycle analysis results indicate that there are highly proliferative cycling cell subsets, in which cells are in S and G2/M phases only, in the tumor cell population, myeloid cell population, T cell population, and stromal cell population. The cycling tumor cells are likely the source of tumor proliferation while the cycling immune cells are reactive to killing the tumor cells.

Palisading epithelial cells consist of a proliferating subset and whorl-like clusters are not tumor stem cells

The cluster of the tumor cells other than the WC cells shows strong heterogeneity. In order to further characterize them, all the tumor cells were extracted and subject to re-clustering, which resulted in 6 subpopulations (Fig. 2 and Suppl. Table 3). The cluster 3 (T3) is the WC cells and the other tumor cells are divided into 5 clusters. The expression level of CTNNB1 is the highest in T3 and decreases in the clusters 1 and 6 (T1 and T6) while the expression is very low in the other three clusters (T2, T4 and T5). In spite of the significant differences in the CTNNB1 expression level among these subpopulations, the discrepancy is much smaller in the proportion of the mutated cells among the cells where CTNNB1 exon 3 was detected, ranging from 31% in T5 to 74% in T3, which are not very far from 50%, the expected proportion of a heterozygous mutation, indicating that these cells are indeed tumor cells. T6 is the subpopulation of cycling cells and they express cell proliferation and division marker genes, e.g., TYMS, UBE2C, CKS1B, CENPE, TOP2A and MKI67. The cell-cycling hallmark pathways, such as ‘G2M_CHECKPOINT’, ‘E2F_TARGETS’ and ‘MITOTIC_SPINDLE’, are activated in this cluster. Other than these cell-cycle marker genes, the expression profile of T6 is similar to that of the T1 subpopulation. SFRP1, FRZB (also called as SFRP3), and a few other genes that negatively regulate the WNT signaling are specifically expressed in T1 and T6.

Fig. 2

Subpopulations of the tumor cells. a The UMAP and re-clustering result of the tumor cells. Clusters are denoted by color. Cluster 3 is the whorl-like cluster tumor cells and cluster 1 is the PE cells. b The UMAP are colored by the mutational status of the CTNNB1 exon 3 in each cell (the color and legend have the same meaning as in Fig. 1). c The proportion of cells with sequenced reads covering the CTNNB1 exon 3 (top) and the proportion of mutated cells among the cells sequenced with reads covering the exon 3 (bottom). d The UMAP are colored by the inferred cycling phases. Cluster 6 is the cycling tumor cells (circled with the dotted line). (e) Dot plots showing the log-transformed expression of four representative marker genes. DKK4 and CXCL14 are mainly expressed in cluster 3, and FCSP and S100A8 are mainly expressed in cluster 5. f Violin plots showing the distribution of the percentages of mitochondria reads (top) and ribosomal reads (bottom) in the six tumor clusters. g Bubble plots showing the log-transformed expression of representative marker genes for clusters 1, 3, 5 and 6, where the size represents the percentage of cells expressing the genes and color hue represents the log-transformed mean expression levels. f Heatmap of the pathway enrichment values in each cell calculated with the scAUC algorithm for cancer hallmark pathways and tooth-related GO pathways

To further identify these tumor cells, we applied the ‘anchor’-based integration workflow in the Seurat package to map the tumor subpopulations to spatial transcriptomic spots (Fig. 3). The results reveal that T1, the most abundant subpopulation, corresponds to the palisading epithelial (PE) cells and the adjacent squamous cells (those not orderly aligned as PE) as observed in the hematoxylin & eosin (H&E) staining images of the Visium slide from P455. The T3 subpopulation is mapped to the whorl-like cluster spots, confirming that T3 are indeed WC cells. The WNT/β-catenin and non-canonical WNT signaling pathways and hedgehog signaling pathway are activated in T3. Most of the cells in this subpopulation are in the G1 phase, without proliferating cells. Multiplex immunofluorescence staining also shows that Ki67-positive cells are distributed only in the PE cells, not found in the WC cells, indicating that proliferating tumor cells originate from the former not the latter. The genes SOX2 and SOX9, which are commonly used to label stemness, are mainly expressed in the T1 subpopulation as well, with little expression in T3 (Suppl. Figure 1). Spatially, the two genes are also expressed in the PE spots only, not in the WC spots (Fig. 3). CD44 has been widely implicated as a cancer stem cell (CSC) marker, but it is widely expressed in the stromal and immune cells at a higher level than in the tumor cells. Although the expression level is higher in T3 than in T1, it is also expressed in T2 and T5 at a similar level (Suppl. Figure 1). This also indicates that in ACPs, the WC cells are unlikely the tumor stem cells.

Fig. 3

Spatial location of the subpopulations of the tumor cells. a Prediction scores for each spot in the Visium slide of P455 for cluster 1 of the tumor cells (T1) in the scRNA-seq data and the expression levels of SOX2 and SOX9 in each spot. b Prediction scores for cluster 5 of the tumor cells (T5) and the expression levels of FDCSP and S100A9 in each spot. c The zoomed image of the square in the leftmost hematoxylin and eosin (H&E) staining image in (b) where representative tumor cells of T2, T3 and T5 are labeled. Prediction scores for T2 and T3 are shown on the bottom left. T3 are mapped to the whorl-like clusters and T2 to the densely packed tumor germinal centers. d Multiplex immunofluorescence imaging of S100A9, KI67, DAPI in an ACP specimen

Senescent tumor cells are identified with specific expression of inflammatory factors and they have unique cell morphology

It has been hypothesized that the seeding cluster cells in ACPs drive the formation of tumors through a paracrine mechanism by secreting inflammatory and growth factors to form a senescence-associated secretory phenotype (SASP) [9, 29]. Although the T3 cells express inflammatory cytokines, e.g., CXCL14, and fibroblast growth factors, such as FGF4 and FGF19, they have the lowest SASP enrichment score calculated with the SenMayo gene set [31]. This indicates that the WC cells are not the senescent cells either. Instead, we identified a novel subset, the T5 subpopulation, in which the SenMayo gene set are enriched (Fig. 2h). Many inflammatory factors, e.g., FDCSP, S100A2/8/9 and IL1RN are expressed in this subset only (Fig. 2g). FDCSP is a secreted protein expressed in follicular dendritic cells which specifically binds to activated B cells [22]. S100A8 and A9 encode proteins to form calprotectin which is a calcium- and zinc-binding protein and plays a prominent role in the regulation of inflammatory processes and immune response, such as neutrophil chemotaxis and adhesion [36]. Inflammatory hallmark pathways are activated in this subset, including the ‘Complement’ system, the ‘IL6_JAK_STAT3’ signaling pathway, the ‘TNFα_via_NFκB’ signaling pathway, and the interferon α and γ response pathways, indicating that this subpopulation exhibits an inflammatory phenotype.

On the Visium spatial slide of P455, T5 cells were not mapped to the usual PE or WC spots, but instead to the loose squamous epithelial cell region (Fig. 4) where FDCSP and S100A8/9 are expressed. From the H&E staining image, it can be seen that T5 shows different morphological features from the usual PE or WC cells. Their nuclei are elongated to a narrow rectangle and they are not arranged in order. Immunofluorescence staining also found that S100A9 positive cells are at different regions from the PE cells (Fig. 3). In addition to the inflammatory pathways, estrogen signaling pathways are also activated. These results indicate that the cells in the T5 subpopulation are likely the SASP cells.

Fig. 4

Molecular characteristics of manually selected spots in the Visium slides. a The selected spots which are annotated as T1, T3, T4 and T5 in the slide of P455 (top) and the scaled expression of the marker genes in each spot (bottom). Three zoomed images show the regions for one whorl-like cluster (T3) and two tumor germinal centers (T2). b The selected spots of four histological components, PE, SR, WC and WK, in the slide of P457 (top) and the scaled expression of the marker genes in each spot (bottom). c The selected spots which are annotated as T1 and T5 in the slide of P452 (top) and the scaled expression of the marker genes in each spot (bottom). Two zoomed images show the regions for one WC region surrounding WK along a thin layer of PE (left) and the SASP cells (T5). The center of each circle is located at the center of the corresponding spot and the diameter is 100 µm (inter-spot distance)

In the scRNA-seq data, T5 cells were found in all the tumors except in P433, and in 7 specimens there are more than 20 cells of T5 (Suppl. Table 3). As a comparison, only in P432 and P455 there are more than 10 cells of T3. On the spatial slides of the other two samples, the unique regions of T5 cells were also found (Suppl. Figures 2 and 3). These results suggest that, like PE cells, the senescent cells of T5 are ubiquitously present in ACPs.

Single-cell and spatial sequencing also identifies tumor germinal cells

T2 is the second most abundant tumor cells. The scRNA-seq profiles of T2 and T4 subpopulations are similar as that of T1, and most of them were mapped to the adjacent region to T1 on the spatial slides. The proportion of mitochondrial genes in T2 is higher (median is 10%) than those in the other subpopulations, while the proportion of ribosomal genes in T4 is higher (median is 40%) than those in the other tumor cells (Fig. 2f). Therefore, like PE, they are the epithelial cells at different states.

In the Visium slide of P455, there are a significant number of tumor spots identified as T2 (Fig. 3). In particular, there are a few isolated spots with unique morphology in the stromal region with no other adjacent PE or squamous cells (Fig. 4). These cells are tightly packed together in a much higher density (approximately 1 cell per µm2) than other tumor cells. The size of each cell is small while the nucleus represents the majority of a cell, indicating that these burgeoning cells are the tumor germinal centers. Furthermore, the packing density is looser in the middle than that in the borderline suggesting that these germinal centers are expanding. Although the number of reads sequenced in these T2 spots are much smaller than those in T1 spots, the top-expressed genes in these T2 spots are similar with those in T1 spots.

We applied the RNA velocity algorithm, scVelo, and the monocle pseudo-time analysis to infer the development trajectory of tumor cells in P455, the specimen in which there are abundant cells of all 6 subpopulations. Both algorithms revealed that there are two different development trajectories starting from the T2 cells (Fig. 5). Figure 5a shows the velocity streamlines, which describes the speed and the direction of the state transitions. One root state is found at the intersection between T1 and T2 clusters. Besides the major development trajectory in T1, there is another trajectory which ends in T3 (the white arrow). The monocle’s pseudo-time analysis also reveals a similar two-branch development trajectory. The trajectory starts from T2 and develops into T1 initially, and then branches into two fates, in which T3 is one terminal state, while the other terminal state consists mainly T1 and T5 cells. Similar results were also obtained in P432 (Suppl. Figure 4a). This result is consistent with that the T2 cells are the germinal tumor cells.

Fig. 5

Developmental trajectories and transcriptional factor analysis results of tumor cells. a The embedded velocity streamlines in the RNA velocity analysis of the tumor cells in P455, showing major directions of cell progression in the transcriptional space. The zoomed regions show the source (left) and the drain (right) in the velocity field. The white arrow shows the trajectory from T1 to T3. b The trajectory of the tumor cells in P455 based on the pseudo-time analysis in Monocle2. The cells are colored by the cell cluster (top, the same color scheme as in a) and the pseudo-time (bottom). The arrows indicate the two trajectory branches. c The heatmap of the activity (scaled AUC scores) of the cluster specific regulons. d The violin plots of the AUC scores for two regulons, PITX1 and SOX15, in 5 tumor clusters. The plots for the stromal cells are also shown for comparison in which the two regulons are inactive

During the development process, the expression level of CTNNB1 increases from the root to both the terminal states. Along the two branches of the trajectory, different sets of inhibitory genes of the WNT/β-catenin signaling are upregulated, DKK4 and NOTUM in the branch to WC, and SFPR1 and FRZB in the other branch (Suppl. Figure 4 and Suppl. Table 4). Different cytokines, e.g., CXCL14 in one branch and S100A10 in the other, are also associated with the development trajectories.

Regulons of inflammation factors are activated in the senescent tumor cells

The SCENIC algorithm was used to determine the activity of regulons in the scRNA-seq data. The cluster-specific regulons are shown in Fig. 5c. Several inflammation-associated regulons are activated in T5, including interferon regulatory factors 5 and 7 (IRF5 and IRF7), NFYC which regulates the expression of major histocompatibility complex (MHC), and MAFB which is a critical regulator of complement component C1q. In WC (T3), the regulon of TCF7 is activated, and this regulon also acts as feedback transcriptional repressor of CTNNB1 besides forming a transcription complex with β-catenin. In T2, the activated regulons include LEF1, SOX17 and others. LEF1 is a key nuclear mediator of Wnt/β-catenin signaling and SOX17 inhibits Wnt signaling. The PITX1 regulon is activated in all tumor cells except that its activity is stronger in T1 (Fig. 5d) than in other subpopulations. PITX1 is a transcriptional regulator involved in basal and hormone-regulated activity of prolactin. The SOX15 regulon is activated in all tumor subpopulations except T3. This result suggests that targeting the transcription factors specific to the senescent tumor cells has the potential to inhibit their growth.

Odontogenesis associated genes are expressed in the senescent tumor cells

Although in the scRNA-seq data the gene ontology (GO) biological process (BP) pathways of the regulation of odontogenesis, the regulation of chondrocyte proliferation and chondroblast differentiation, are activated in T3, the pathways of tooth eruption, tooth mineralization, regulation of vascular endothelial growth factor production and vascular endothelial growth factor production are activated in T5, not in T3 (Fig. 2h). Among odontogenesis related genes, AMBN is specifically expressed in T3, other genes including ENAM, AMELX and AMELY are not expressed in all the tumor cells.

Representative spots of T1, T3, T4 and T5, annotated with the scRNA-seq data were manually picked from the Visium slides to identify marker genes. As a comparison, representative spots of PE, SR, WC and WK were manually selected from the slide of P457, which are rich in various histological components, to identify marker genes in each of them (Fig. 4). From the expression heatmaps of the marker genes in Fig. 4, we can see that T1 are indeed PE cells with specific expression of EPHA7, TLE2, WIF2, PITX2, FGFR2, EGR1, etc. The T3 spots overlap with the WC regions, with characteristic expression of NOTUM, DKK4, SHH, FGF4, KRT23, etc. Besides a number of cytokines, e.g., FDCSP, IL1RN, LIF, other highly expressed genes in the T5 spots include some epidermal growth factors, e.g., AREG, proteinase MMP13, odontogenesis genes, e.g., ODAM and AMTN, and laminin genes, e.g., LAMA3, LAMB3, LAMC2. The odontogenesis associated genes, e.g., ENAM, AMELX, AMTN, AMBN, and the laminin gene LAMC2 are expressed in the WK spots not in WC spots. These results suggest that senescent cells might play a critical role in odontogenesis while the WC cells regulates the process.

Infiltration of immune cells and clonal expansion of cytotoxic T cells in the tumor microenvironment

In the scRNA-seq data, a total of 21,866 lymphoid immune cell expression profiles were obtained. Cluster analysis identified 12 subclusters, including T cells, natural killer T (NKT) cells, natural killer (NK) cells, B cells, mast and plasma cells (Fig. 6). A subcluster of cycling T cells was found, indicating the presence of rapidly proliferating T cells. Other T cell subpopulations include CD8 T effector memory cells (TEM), regulatory T cells (Treg) and CD4 naïve T cells. GZMK is highly expressed in CD8 TEM (clusters 1 and 9) and the expression levels of GZMB and GZMH are higher in NKT (cluster 7) than in CD8 TEM while GZMA and CCL5 are expressed in both the CD8 TEM and NKT subpopulations. The results show that they are the cytotoxic effector cells. Tregs (cluster 5) show high expression of TIGIT, CTLA4, FOXP3, ICOS, IL2RA and TNFRSF4.

Fig. 6

Subpopulations and clonal expansion of lymphocytes. a The UMAP of the re-clustering result of lymphoid cells, colored by the cell types. b Violin plots showing the distribution of the percentages of mitochondria reads (top) and ribosomal reads (bottom) in the clusters. c Bubble plots showing the log-transformed expression of representative marker genes, where the size represents the percentage of cells expressing the genes and color hue represents the log-transformed mean expression levels. d Clonal proportion of top 10 clonotypes of T cells of which the TCR was retrieved in each sample (see Suppl. Table 5). e The UMAP are colored by the number of TCR contigs in each cell. (f) The UMAP are colored by the top 4 clonotypes in P416 (clonotype 1, 168, TRA:CATGNSGNTPLVF, TRB:CASSLMGQAMGELFF; clonotype 2, 83, TRA:CAMRGIRSGGSNYKLTF, TRB:CSASPPGGVGANVLTF; clonotype 3, 80, TRB:CASSLMGQAMGELFF; clonotype 4, 57, TRB:CSASPPGGVGANVLTF). TCR: T cell receptor; CDR3, complementarity-determining region 3

We further integrated the scTCR-seq data with the scRNA-seq results. Although most of the clonotypes are singletons, i.e., detected in only one cell, there are a few clonotypes in each specimen showing a high level of clonal expansion, particularly in P416 and P455 (Fig. 6d, and Suppl. Table 5). Most of the TCR contigs are found in the T cells, supporting the correct assignments of the T cells. We colored the cells carrying with the four largest TCR clonotypes of P416 in Fig. 6f. A majority of them are located in the CD8 TEM subpopulation, suggesting that these clonal-expanded cells are indeed cytotoxic T cells.

The scRNA-seq profiles of a total of 25,532 myeloid immune cells were obtained, in which the microglial cells are the predominating subpopulation, with specific expression of P2RY12 and a higher-level expression of GPR34 and chemokines, e.g., CCL3, CCL4, CCL3L1 and CCL4L2, than other cells (Fig. 7). The other myeloid subpopulations include monocytes, macrophages, DCs and a few neutrophil cells. In monocytes there are a small cluster of cycling cells. The major histocompatibility complex (MHC) class II molecules are expressed at a higher level in DCs than in microglial cells or monocyte/macrophages, but the CD1 family genes (C1QA, C1QB and C1QC) are mainly expressed in microglial cells, monocytes and macrophages, and the expression levels are higher in microglial cells and macrophages than in monocytes. The CD1 molecules are known presenting a broad range of lipid, glycolipid and lipopeptide antigens.

Fig. 7

Subpopulations of myeloid cells. a The UMAP of the re-clustering result of myeloid cells, colored by the cell types. b The UMAP colored by the inferred cell cycle phases. The circled subpopulation is the cycling cells in the G2M/S phases. c Scatter plots of proportions of B cells versus those of the microglial cells in 12 ACPs. The 3 samples in the bottom right square are microglia dominant tumors and the 3 samples in the top left square are the B-cells dominant tumors while both are absent in the other specimens. d Bubble plots showing the log-transformed expression of representative marker genes, where the size represents the percentage of cells expressing the genes and color hue represents the log-transformed mean expression levels. e Score distribution of B cells, T cells, plasma cells and monocytes/macrophages on the spatial slide of P455 (see Fig. 4a)

We notice that there is a mutual exclusion between B cells and microglial cells in ACPs (Fig. 7c). Most of the microglial cells come from P416, P432, and P418 tumors, in which the proportion of B cells is very low. In P419, P421, and P457, there is high abundance of B cells, but the proportion of microglial cells is very low. In the other six tumors, both cell types are relatively lacking.

The spatial distributions of several major types of immune cells in P455 are illustrated in Fig. 7e. Besides the accumulation of monocytes/macrophages in the inflammatory stromal regions, we observe that B, T and other immune cells cluster around the T5 cells. There is also clustering of these immune cells adjacent to the T1 region, whereas the scores for these immune cells are notably low in the spots adjacent to the T3 cells. A few spots with high T-cell scores are scattered among the T1 spots. These results suggest that the SASP cells possess the ability to recruit immune cells, whereas the WC cells lack such capability.

Comparison between the primary and relapsed tumors

Among the 12 ACPs studied, 7 are the primary tumors and the remaining 5 are the relapsed tumors. Although an increasing trend is observed in the proportion of T1 tumor cells in the relapsed tumors, this difference is not significant (Wilcoxon’s P = 0.6, Fig. 8a). In contrast, the proportion of T3 cells is significantly higher in the relapsed tumors compared to the primary ones (Wilcoxon’s P = 0.05). The β-catenin staining results reveal that there is also an increasing proportion of nuclear positive cells, indicator of WC cells, in the relapsed tumors, but the difference is not significant (Wilcoxon’s P = 0.6, Fig. 8b). Additionally, we noted that the proportion of nuclear β-catenin positive tumor cells in the IHC results is significantly higher than that of T3 cells in the scRNA-seq dataset, reflecting the bias of single-cell sequencing technique which requires live tumor cells.

Fig. 8

Comparison of tumor microenvironment between the primary (P) and relapsed (R) ACPs. a Box-whisker plots of the proportions of three tumor subpopulations in the scRNA-seq dataset grouped by the tumor status (P vs. R). b Box-whisker plot of the proportions of nuclear β-catenin positive tumor cells. c Box-whisker plots of the proportions of monocytes, B cells and NKT cells in the scRNA-seq dataset. d Box-whisker plot of the numbers of the CD20 + cells. e Typical IHC images of CD20 staining in primary (top) and relapsed (bottom) tumors. Scale bar, 50 μm. The P values were calculated with the Wilcoxon rank sum test

Contrary to T1 and T3, the proportion of T5 cells is significantly reduced in the relapsed tumors (Wilcoxon’s P = 0.05). This observation may stem from the fact that T5 cells are at a later developmental stage, rendering their proliferation rate less as rapid compared to T1 cells.

Among the immune cell populations, the proportion of B cells is significantly lower in the relapsed tumors compared to the primary tumors (Wilcoxon’s P = 0.02, Fig. 8c) while the proportion of NKT cells is significantly higher (Wilcoxon’s P = 0.002). IHC staining against CD20 supported the trend of the B cells (Wilcoxon’s P = 0.08, Fig. 8d, e). Additionally, the infiltration of T cells is also lower in the relapsed tumors although the difference is not significant (Suppl. Figure 5). These findings are in line with the observation of a lower proportion of T5 tumor cells in the relapsed tumors, suggesting potential alterations in the immune landscape during tumor relapse.