Gliomagenesis mimics an injury response orchestrated by neural crest-like cells

Mouse models

The transgenic mice used in this study were obtained from Jackson Laboratories, with the following exceptions: Sox2^creER (B6;129S-Sox2^{tm1(cre/ERT2)Hoch}/J) from Konrad Hochedlinger⁵⁶; Trp53^f/f from Chi-chung Hui⁵⁷; Sox2eGFP (Sox2^tm1Lpev) from Freda Miller⁵⁸. The glioma and control mouse models were generated and housed in 12-h dark–light cycle facilities maintained at appropriate temperature and humidity and in which mice had free access to water and chow. The mice were monitored daily and euthanized once they developed end-point symptoms of raised intracranial pressure or focal neurological abnormalities. The mouse experiments were all performed following the ethical and legal regulations. The experiments and animal use protocols were approved by the Animal Care Committees in the different institutions at the University of Toronto, including the Hospital for Sick Children and University Health Network.

MRI and collection of mouse tissue samples for single-cell profiling

MRI was performed as reported previously⁸, and fresh brain tissue was collected from mutant mice at four MRI-defined stages: the ‘preneoplastic’ stage, at which the brain imaging shows no signs of neoplastic lesion development; the ‘early-lesion’ stage, characterized by small abnormalities seen on T2-FLAIR MRI sequences; the ‘mid-lesion’ stage, when the lesion has reached a larger size, as indicated by a T2-FLAIR-bright mass, in asymptomatic animals and occupies a substantial fraction of the brain hemisphere; and, finally, the ‘end-point’ stage, when mice develop symptoms of raised intracranial pressure or focal neurological abnormalities, with the tumour extending over a large portion of the brain hemisphere(s), typically with midline shift. Each brain tissue sample was dissociated into single cells as described previously⁸. This was followed by FACS to separate the live tdTomato⁺ cells and tdTomato⁻ cells. Following sorting, both populations were characterized by scRNA-seq using the 10x Genomics platform as reported previously⁸. Tissue processing was performed in the same way for the Sox2eGFP mouse brain injury samples and was followed by cell sorting to separate the live GFP⁺ and GFP⁻ cells. For scATAC-seq, the nuclei were isolated from sorted cells and processed using the 10x Genomics Single Cell ATAC-seq workflow and following the manufacturer’s protocol (user guide CG000168 Rev A). In brief, 5 μl of nucleic suspension was added to the transposition reaction, which was used for the gel beads-in-emulsion (GEMs) generation and barcoding on the 10x Chromium Chip E for PCR amplification. DNA was recovered using Dynabeads MyOne Silane beads, which proceeded to library preparation to add Chromium i7 Sample Index N Set A library barcodes. Libraries were validated on the Agilent 2100 Bioanalyzer to check for size and quantified by quantitative PCR using Kapa Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems). Validated libraries were pooled in equimolar quantities and paired-end-sequenced on the Illumina Novaseq platform following Illumina’s recommended protocol to generate paired-end reads of 50 bases in length. Spatial transcriptomics was performed on formalin-fixed and paraffin-embedded mouse tissue sections using the 10x Genomics Visium platform and following the manufacturer’s protocols (user guides CG000407 Rev D and CG000409 Rev C).

Implantation of intracranial cannula and osmotic pumps

Mice were prepared for surgery and anaesthetized using isofluorane. The intracranial cannula (Brain infusion kit 3) was implanted at 1.5 mm lateral and 0 mm posterior to bregma, and the saline osmotic pump (Alzet Model 1007D) was implanted subcutaneously following the manufacturer’s instructions. Each osmotic pump was filled with normal saline to infuse the brain at a flow rate of 0.5 µl h⁻¹ for 5 consecutive days.

Analysis of the mouse scRNA-seq tumorigenesis atlas

Raw sequencing data from the 30 scRNA-seq samples (Extended Data Table 1) were aligned to the mm10 mouse reference using Cell Ranger (v3.1.0). Doublets and multiplets were identified by scDblFinder (v1.4.0), and low-quality cells (percentage of mitochondrial >12%; number of genes detected per cell <800; number of unique molecular identifier (UMI) per cell <500) were removed as part of the quality control process. Seurat (v4.5) was used to integrate and cluster data from all samples together. We used Seurat’s merge function to combine data from all samples together. For clustering, we used the following: Seurat’s LogNormalize method with a default value of 2,000 variable features; scaling was performed with regression of mitochondrial percentages; rann method for finding nearest neighbouring cells; and principal component analysis (PCA) reduction with the first 30 dimensions. Cell type assignment was based on the differentially expressed gene analysis and using the known cell type markers from the literature. We computed the top differentially expressed markers for each of the clusters compared to all other clusters using Seurat’s Wilcoxon rank sum method with a minimum cell fraction of 0.25 and a minimum fold difference of 0.25, and then ranked them by their fold changes and adjusted P values.

Analysis of the malignant cells

We used Seurat’s subset function to isolate the malignant cells, which belong to the Sox2CEPPT and nestinCPPT mouse model samples, from the main atlas. The following parameters were used to recluster these malignant cells: data normalization was performed with the LogNormalize method and a scale factor of 10,000; a default value of 2,000 variable features with the vst selection method was used to determine the top variable features; data scaling was performed with regression of mitochondrial percentages. We used Harmony (v1.0)⁵⁹ to generate a batch-corrected embedding space for downstream analysis. The rann method was used to identify k-nearest neighbouring cells. We then annotated the cells on the basis of the differentially expressed gene analysis and using the known cell type markers from the literature.

To guarantee that downstream analysis is not biased by the integration method, we used FastMNN⁶⁰, implemented in SeuratWrappers, and BBKNN⁶¹, implemented in scanpy⁶², to redo the batch correction. These methods all gave similar results (Supplementary Fig. 6).

Analysis of the cycling PC-like cells

To explore the heterogeneity within the cycling PC-like cells, we first isolated these cells from the malignant cells and then removed all cell cycle mouse orthologues associated with both S and G2/M phases⁶³ as well as all of the genes from the KEGG (Kyoto Encyclopedia of Genes and Genomes) cell cycle pathway reference (https://www.genome.jp/entry/pathway+mmu04110). Removing these 205 cell cycle genes was necessary to mask out cell cycle genes expressed at high levels to identify the cycling subtypes. Reclustering revealed six cycling cell types (cycling NSC-, cycling OPC-, cycling NPC-, cycling MSC-, cycling AC- and cycling NCC-like cells). To compute the fractions of cycling and non-cycling cells for each of the precursor-like malignant states (NSC, OPC, NPC, MSC and NCC), we calculated the number of cycling and non-cycling cells for each precursor state in each of the samples. The fractions of cycling and non-cycling cells for each PC-like state were calculated as the number of cells over the sum of both cycling and non-cycling cells. We plotted the averages along with their standard errors.

Relative abundance of the malignant cell states

To determine the relative abundance of the individual malignant cell states (NSC-, OPC-, NPC-, MSC-, AC-, NCC-, immature OL-, cycling NSC-, cycling OPC-, cycling NPC-, cycling MSC-, cycling AC- and cycling NCC-like cells) across the four time points in an area plot, we computed the number of cells from all samples for each of the cell states in each of the four time points. We then computed the total number of cells of all states for each of the four time points. For each time point, the fraction of each cell state was calculated as the number of cells for that cell type over the total number of cells in that time point (Fig. 2e). Further, to compare the relative abundances of the cell population categories across the four stages of tumour development, we first combined all of the non-cycling PC-like states (NSC-like, NPC-like, OPC-like and MSC-like) and designated them as non-cycling lineage PC-like. Similarly, we combined the differentiated states (AC-like, cycling AC-like and immature OL-like) and designated them as differentiated-like. The third and fourth groups were the cycling PC-like and NCC-like cells. We then computed the fraction of each of the four populations in each of the tumour samples. We grouped these fractions on the basis of the tumour stages of the samples and computed their averages and standard errors. We plotted the averages along with their standard errors for each of the categories with respect to the four tumour stages (Extended Data Fig. 2b).

Inferring CNAs

We used the R package inferCNV (v1.7.1), inferCNV of the Trinity CTAT Project (https://github.com/broadinstitute/inferCNV), to estimate the CNAs from scRNA-seq data. We adopted the recommended parameters for 10x Genomics data from the inferCNV tutorial (cutoff = 0.1). As a reference, we selected three samples (from P39, P111 and P365 mice) from our normal brain atlas⁸. Running inferCNV resulted in a continuous, gene-level relative CNA profile for each cell. Visualizing the distribution of the relative CNAs, we observed three main peaks, with one bigger peak centring around 1 representing a gene without clear CNAs, and two peaks roughly symmetric about 1, which we interpreted as a gain of copy and a loss of copy, respectively. Theoretically, one copy of gain or loss should cause a shift of 0.5 in the relative copy number, meaning that the centre of the peak should be around 1.5 for a gain of copy and 0.5 for a loss of copy. Owing to the high level of noise in scRNA-seq data, the inferred CNAs are far from perfect, presenting a smaller one-copy shift than 0.5. To infer the absolute copy numbers, we rescaled the inferred CNAs to make the two CNA peaks located around 1.5 and 0.5 before rounding them to integers. Specifically, we first removed genes without clear CNAs (abs(CNA − 1) < 0.01). We then identified the two CNA peaks by clustering the CNA values into two clusters using k-means. We then used the centroids of the two clusters to calculate the one-copy shift in the data (mean(abs(centroid − 1)). Then a rescaling factor can be calculated as the fold change between the theoretical one-copy shift (0.5) and the estimated one from the data. After rescaling the CNA profile, we rounded it to integers to calculate the absolute number of copies (round(CNA × 2)).

CNA neighbourhood analysis

For each of the tumour samples, we stored the relative CNA as a new ‘modality’ (cna) in Seurat and processed with the log-normalization pipeline with the default parameters. We therefore obtained a new set of PCA embeddings, clusters and UMAP visualizations. We observed that most CNA clusters (clones) contain different cell states. To further elaborate on this, we fed the CNA PCA embeddings to the R package miloR v1.5.0⁶⁴ and clustered cells into (overlapped) neighbourhoods. Nodes are neighbourhoods of CNA, with colour indicating cell state of the neighbourhood index cell, and size corresponding to the number of cells in the neighbourhood. Graph edges depict the number of cells shared between neighbourhoods.

Building cell phylogenies (MEDALT)

We used the python package MEDALT²⁵ to build the phylogenetic tree of cells. This was performed separately for different samples. MEDALT takes the absolute copy numbers as input and calculates distances between cells using a distance metric called the minimum event distance (MED). The basic assumption of MED is that gain or loss of a copy of adjacent genes can happen together. After getting a pairwise distance between all cells (which can be treated as a densely connected distance graph), a rooted minimum spanning tree was built on the distance graph, resulting in the final phylogenetic-like tree of cells. We used Cytoscape v3.9.1 (https://cytoscape.org/) to visualize the graph.

Trajectory analysis

The R packages phateR v1.0.7⁶⁵ and slingshot v1.8.0⁶⁶ were used to build trajectories. First, we used phateR to generate two-dimensional (2D) embeddings (PHATE space), with the 50 harmony embeddings as input. Then we built a spanning tree of different cell types in the 2D PHATE space using slingshot. For differential analysis, we first clustered cells into five clusters along the trajectory, and then the cluster-specific maker genes were found using the FinderAllMarkers() function in Seurat. To test the robustness of the results, we also used the diffusion map implemented in the R package destiny v3.4.0 to generate the 2D embeddings of the 50 harmony embeddings, and used PAGA in scanpy to summarize the k-nearest neighbours (knn). All of these methods and others resulted in trajectories with similar topologies between cell states.

Pathway enrichment analysis

We used the R package enrichR (v3.0)⁶⁷ for the enrichment analysis. For the microglia cell types in Fig. 6a, we first find differentially expressed genes for each of the clusters using Seurat::FindAllMarkers(), with a cutoff of P_adj < 0.05 and fold change > 2. For each of the clusters, we fed the differential gene list to enrichR and compared it against the variety of gene-set libraries implemented in enrichR. Representative enriched terms from the gene-set libraries HDSigDB_Mouse_2021, Mouse_Gene_Atlas and MSigDB_Hallmark_2020 were selected to show in the main figure plot, and the complete results can be found in Supplementary Table 8.

Cell–cell communication

We used CellChat v1.1.3³⁸ for the analysis of cell–cell communication, which was performed separately for each of the four stages of tumorigenesis. We followed the official workflow with default parameters unless otherwise indicated. First, we loaded the normalized counts into CellChat, followed by the preprocessing steps identifyOverExpressedGenes() and identifyOverExpressedInteractions(). We smoothed gene expression by applying a diffusion process on the protein–protein interaction network implemented in projectData() function. We then ran the computeCommunProb() function for communication analysis with the parameter population.size = FALSE to eliminate possible bias due to cell population size. This resulted in a network of communication strength between all cell states for each of the ligand–receptor pairs that passed the filtering steps. We used the aggregation functions computeCommunProbPathway() and aggregateNet() to determine the communication strength between cell states at pathway and global levels, respectively. For each of the pathways (data slot netP), we evaluated the role of different cell states as senders or receivers on the basis of the out-degree or in-degree of the communication network, implemented in the netAnalysis_signallingRole() function.

Statistical modelling

We conducted a comprehensive evaluation of how well the cell phylogenies can be reconstructed by the pipeline we adopted using simulation data (Supplementary Figs. 7 and 8). The evaluation pipeline consisted of three parts. The first part consisted of simulating the cell phylogeny and the corresponding CNA profiles. An important consideration is that we are able to see only the living cells on the phylogeny. We truncated the phylogeny to keep only living cells and their CNA profile for downstream analysis. The second part consisted of reconstructing the cell phylogeny based on the CNA profiles of living cells. Besides MEDALT, we also test hierarchical clustering with three different distance measures as the baselines. The third part consisted of evaluating how well the cell phylogeny has been reconstructed from two aspects: the clone identification and the directionality of cell evolution. Our simulation results showed that all four methods can recover the clones reasonably well on the basis of CNA profile, but MEDALT showed better performance in terms of reflecting the directionality of cell evolution.

Simulating cell phylogenies

We used two models to simulate the cell phylogenies, the classic birth–death process and the birth–death process with immigration. The second model has been proved to better explain cell growth dynamics in the tumour context when there is a proliferative hierarchy involving a slow-cycling stem cell-like population²³.

For the classic birth–death process, the expected total number of cells N(t) can be expressed as:

in which N₀ is the initial number of cells, and b and d are the birth and death rate, respectively. The interval of time Δt between two adjacent events (the length of the branch in the phylogenetic tree) follows an exponential distribution with mean E(Δt) = 1/(b + d). When branching happens, it can be a birth with a probability of b/(b + d), or a death with a probability d/(b + d). In our simulation, we assumed that the birth and death rates do not change during the evolution.

For the birth–death with immigration model (Supplementary Fig. 8), we followed the minimal model of tumour growth proposed in ref. ²³. This model is based on a two-component hierarchy involving transitions from a slow-cycling stem cell-like compartment (S) to a more rapidly cycling progenitor population (P). The simulation is similar to the birth–death model; the difference is that cell division happens in a different way for S cells and P cells. For the P cells, the branching is the same as the classic birth–death model, with an equal probability of being a death or a birth. For the S cells, with high probability (0.8–0.9 for our simulation), it will divide asymmetrically, giving rise to an S cell and a P cell. However, with a small probability, it can also divide symmetrically to self-renewal. In our simulation, we assumed that the birth and death rates for P cells and the probability of symmetrical division for S cells do not change during the evolution.

After simulating the whole phylogeny of all cells, we truncated the phylogeny using the function ape::drop.tips() to keep only living cells on the leaf node, with necessary common ancestors to form a complete dendrogram. This truncated phylogeny of living cells was used as the ground truth for downstream analysis.

Simulating CNAs

The CNAs in a cell are acquired in two ways: inherited from its parent cell; newly acquired during cell division. We simulated the CNA accumulation on the whole phylogenetic tree we generated. For each of the cells, we assume that the number of newly acquired CNAs follows a Poisson distribution with the parameter λ. This means that with probability e^−λ, a cell will not get any new CNA. In our simulation, λ was set between 0.1 and 0.5. CNAs can affect contiguous sites and regions in a chromosome, meaning that a gain (or loss) of copy for adjacent genes can happen together. We assumed that the starting position of each CNA was uniformly distributed across the genome, and the number of genes that one CNA affected followed another Poisson distribution (with a mean of 100–200 in our simulation). After each simulation, CNAs for living cells were selected for downstream analysis.

Reconstructing the cell phylogeny for simulation data

We used four methods to reconstruct the cell phylogeny—MEDALT and hierarchical clustering with three different distance measures of CNAs (namely, the Euclidean distance, the MED and the shortest path distance (SD) on the tree given by MEDALT). MED was calculated using MEDALT. We loaded the tree given by MEDALT into R as an igraph object, visualized it using GGally::ggnet2() and calculated pairwise SD using igraph::distances(mode = ‘all’).

Evaluating the reconstructed tree

We evaluated the reconstructed cell phylogenies from two perspectives: clone identification and directionality of evolution. We used MEDALT and hierarchical clustering with three different distance measures of CNAs, namely the Euclidean distance, the MED and the SD on the tree given by MEDALT. MED was calculated using MEDALT. We loaded the tree given by MEDALT into R as an igraph object, visualized it using GGally::ggnet2() and calculated pairwise SD using igraph::distances(mode = ‘all’).

Analysis of the scATAC-seq data

We ran Cell Ranger ATAC v1.2.0 to process the raw sequencing data. The cellranger-atac count pipeline was used to align reads and generate single-cell accessibility counts for the cells, with mm10 (refdata-cellranger-atac-mm10-1.2.0) as the reference genome. We adopted the R package ArchR v1.0.1⁶⁸ for the downstream analysis of the scATAC-seq data, following the developers’ default recommendations, unless otherwise indicated. We applied a two-round iterative latent semantic indexing (IterativeLSI) on the 500-base-pair tile matrix, resulting in a 30D embedding space of the genome-wide accessibility profiles. Then we used UMAP to generate a 2D visualization of the IterativeLSI embeddings. Clustering of cells was also conducted in the IterativeLSI space by building a shared nearest-neighbour graph followed by Louvain clustering as implemented in the R package Seurat. Gene activity scores were calculated on the basis of the local accessibility of gene regions, which includes the promoter and gene body. As we have the paired scRNA-seq samples for both scATAC-seq samples, annotation of scATAC-seq data was performed using the label transfer approach implemented in ArchR. Subsequently, open chromatin peaks were called on the basis of pseudo-bulk replicates of different cell types using MACS2, and differentially accessible regions were identified using a Wilcoxon test implemented in ArchR (false discovery rate < 0.1 and fold change > 1.5). Transcription-factor-binding motifs were annotated using the CIS-BP database, identified using the peakAnnoEnrichment() function implemented in ArchR (P_adj < 1 × 10⁻¹⁰) and visualized as a network using Cytoscape. We also conducted differential analysis of gene activity scores using the Wilcoxon test and enrichment analysis of the resulting gene set using enrichR_v3.0 ().

Analysis of the mouse scRNA-seq brain injury dataset

The same computational pipeline and method as for the tumorigenesis atlas were used to remove doublets, and integrate and annotate the ten injury samples (Supplementary Fig. 9a). The following cutoffs were used to remove low-quality cells: percentage of mitochondrial >12%; number of genes detected per cell <800; number of UMIs per cell <500. We isolated the microglia and OL cells in separate analyses and recomputed their top differentially expressed markers separately. The abundance fractions of the different cell types in the relevant category of cell types (for example, PC types) in each time point were calculated as the number of cells for each cell type over the total number of cells from all the cell types belonging to that category in each time point, respectively.

Fisher exact test related to Fig. 5c and Extended Data Fig. 8c,d

We used the R package rstatix (v0.7.0) to compute the one-sided Fisher’s exact tests of the proportions of the cell types in the different time points. For each time point, we calculated the number of cells in each cell category to construct the contingency table. The function fisher_test() was used to perform the statistical tests.

Analysis of the spatial transcriptomics datasets

We analysed spot gene expression in Visium data using Scanpy v1.9.3 and Squidpy v1.3.0. For cell type deconvolution within each spot, we used cell2location v0.1.3, using our single-cell mouse datasets as the reference. For mouse Visium samples containing both normal and tumour regions, an initial reference was constructed using the major cell types denoted in Fig. 1c. The resulting analysis successfully delineated tumour from normal regions, corroborating results from both gene expression and H&E imaging. To attain a more detailed deconvolution of cellular subtypes, two distinct references were created for the tumour and normal regions. For the tumour reference, malignant cells and microenvironment cell types (for example, tumour microglia, tumour OLs, microglia, OLs, MDSCs, T cells, macrophages and endothelial cells) were extracted. Subtypes labels were used to build the tumour reference. For the normal reference, normal brain cell types were utilized to construct this reference.

When estimating the reference cell type signatures, the ‘sample ID’ served as the batch_key to mitigate batch effects. Additionally, the mitochondrial gene ratio and total UMI count were incorporated as covariates to counteract potential biases. After these references were established, we performed the deconvolution separately for tumour and normal regions using their respective references. For human samples, as all tissues were sourced from tumour regions, the tumour reference was used for deconvolution.

Inferring CNV from the Visium data

To investigate the clonal distribution of tumour cells in the space, we applied inferCNV to the Visium data. This analysis was performed under the assumption that the cells in each of the Visium spots belong to the same clone (that is, genetic homogeneity within each spot). Given the small number of cells contained in each spot, we consider this as a reasonable assumption.

We selected the reference data for inferCNV with precision, recognizing its crucial role in ensuring result accuracy. For the mouse Visium, we used the normal regions within each slide as the reference, which we can easily distinguish through H&E imaging and scRNA-seq data. We clustered spots on the basis of their RNA profiles. Among the RNA clusters corresponding to normal regions, we earmarked one as the reference and set aside the remaining normal clusters as negative controls. Our analysis using inferCNV revealed no distinct CNVs in these normal region clusters, whereas tumour regions exhibited pronounced CNVs.

For the human Visium samples, which consist of late-stage tumours with almost no normal regions, we turned to the normal brain Visium data from the SpatialLIBD Project. To mitigate potential reference biases, we randomly selected two slides from the project as references, conducted inferCNV independently for each, and compared the results. Both references yielded congruent outcomes.

Analysis of the human GBM stem cell scRNA-seq dataset

The human GBM stem cell dataset (20 samples) from ref. ⁵³ was downloaded from the Single Cell Data Portal and analysed using the same analysis pipeline and parameters as for our tumorigenesis atlas dataset except with human orthologues. The same quality control cutoffs as for our atlas were used; Seurat and Harmony were used to integrate, correct batch effects, and cluster cells together. The same computational pipelines as for our atlas were used for computing the top differentially expressed genes and inferring CNAs.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.