Animals
All animal procedures were carried out in accordance with the Animal Scientific Procedures Act, 1986 and approved by the UCL Animal Welfare and Ethical Review Body (AWERB) in accordance with local ethical and care guidelines and the International guidelines of the Home Office (UK). Mice used in this study were WT C57BL/6NCrl (Charles River), congenic Sarm1tm1Aidi (strain 018069; Sarm1−/−)38 and Sarm+/+ littermates (both a gift of M. Coleman), CRISPR knockout Sarm1em1.1Tftc and Sarm1 WT39 and B6.Cg-Tg(Thy1-YFP)16Jrs/J (Jax Laboratories, 003709)28. Genetically and aged matched animals were used as controls for Sarm1−/− experiments. NOD.CB17-Prkdcscid/NCrCrl (NSG, Charles River) were used for generation of the PDX models through orthotopic injections of patient-derived GBM cell lines. Tissue from Apptm3.1Tcs (strain 5637817) and rTg4510 (strain 024854) mice (gifts from S. Hong and G. Schiavo) were used as positive controls for assessment of proteinopathies. Mice were group-housed (where possible) in individually ventilated cages and maintained under 12 h–12 h light–dark cycles at 20–24 °C, 40–60% humidity, with water and chow available ad libitum. Mice of both sexes were used and, where appropriate, all animal experiments were blinded.
Derivation and culture of cell lines
Cell lines were derived from the CRUK glioma cellular genetics resource (GCGR) with driver mutations shown in Supplementary Table 1 (G. Morrison et al., manuscript in preparation). GBM2 was derived independently as previously described52. Informed consent was obtained from all of the participants. The study was approved by the National Research Ethics Committee (Wales REC 6; reference 20/WA/0251), and all procedures were conducted in accordance with the ethical standards of the approving committee, the Declaration of Helsinki, the Human Tissue Authority and the General Data Protection Regulation. All patient lines were cultured adherently in serum-free GSC medium (N2 (1/200), B27 (1/100) (Life Technologies), 1 mg ml−1 laminin (Merck, L2020), 10 ng ml−1 EGF (Biotechne, NBP2, 35176), 10 ng ml−1 FGF-2 (Biotechne, NBP2, 35152), 1× MEM NEAA (Thermo Fisher Scientific, 12084947), 0.1 mM 2-mercaptoethanol (Thermo Fisher Scientific, 31350010), 0.012% BSA (Thermo Fisher Scientific, 15260-037), 0.2 g l−1 glucose (Merck, G8769), 1,000 U ml−1 penicillin–streptomycin (Merck, P0781). All of the cell lines were mycoplasma negative.
Generation of somatic and orthotopic PDX models
Somatic tumours were generated as previously reported19,20. In brief, plasmids described in Extended Data Fig. 1a were injected into the right ventricle of isoflurane-immobilized pups at postnatal day 2 using an Eppendorf Femtojet microinjector (Eppendorf, 5247000030) followed by electroporation (5 square pulses, 50 ms per pulse at 100 V, with 850 ms intervals). The EF1a-tdTomato-only plasmid (tdTom) was generated by SnaBI and PmeI digestion of npp plasmid to remove Nf1, Pten and Trp53 guide RNAs before religation. piggyBase (hGFAPMIN-SpCas9-T2A-PBase, 1 mg ml−1) and piggyBac vector U6-Nf1,Pten,Trp53-EF1a-tdTomato (npp, 0.564 mg ml−1) or EF1a-tdTomato (0.423 mg ml−1) were diluted in saline (0.9% NaCl) and mixed at a molar ratio of 1:1. Then, 0.1% fast green (Sigma-Aldrich, F7258) was added to the mix to visualize the injection. All experiments using somatic mouse models were performed on a mix of male and female mice. Whole litters of mice were injected with plasmids, regardless of sex.
To pharmacologically inhibit SARM1 protein selectively in neurons from npp tumour initiation, AAV8-Syn-SARM1-CDN-EGFP (a gift from J. Milbrandt, 0.745 × 1013 viral genomes (vg) per ml)43 was added to the piggyBase/npp piggyBac plasmid mix before intraventricular injection and electroporation, as described above. AAV8-Syn-GFP (Addgene, 50465-AAV8 0.745 × 1013 vg per ml) was used as a control. To investigate SARM1 function after tumour initiation, WT mice bearing tumours were injected with either 2.5 µl AAV8-Syn-SARM1-CDN-EGFP (1 × 1013 vg per ml) or AAV8-Syn-GFP (1 × 1013 vg per ml) at the intermediate or late disease stage. Male and female mice were randomized separately into GFP and SarmDN groups. In brief, mice were anaesthetized and mounted onto a stereotaxic frame. A small craniotomy was performed on the tumour ipsilateral right side of the skull 1.7 mm lateral to bregma and −0.5 mm anterior to bregma. Virus was injected through a 5 μl Hamilton syringe attached to a pump (Pump 11 Elite Nanomite, 70-4507, Harvard Apparatus) at a speed of 0.3 µl min−1 to a depth of 2.4 mm. The virus was injected continuously as the needle was introduced and removed. The wound was sutured and the mice were allowed to recover. Then, 4 weeks after injection, the mice were given an intraperitoneal injection of EdU (5 mg per kg) 2 h before brains were collected following transcardial perfusion with 4% paraformaldehyde (PFA) under terminal anaesthesia.
Orthotopic PDX models were generated as previously described in female NSG mice9. All tumour-bearing mice were monitored daily and euthanized at the required timepoints or at terminal stage, as defined by them reaching tumour-associated humane end points, which correlate with a lethal disease stage, specifically, any one or more of the following signs of general pain or distress (including seizures); greater or equal to 15% weight loss, hunched posture, piloerection, inactivity, ocular/nasal discharge, intermittent abnormal respiratory pattern or loss of body conditioning).
Injury induction in tumour-bearing mice
Brain injury experiments were carried out on tumour-bearing mice at 4.5, 8.5 and 12.5 weeks after electroporation with piggyBase/npp piggyBac plasmids. Male and female mice were randomized separately into sham and injury groups. Mice were anaesthetized and mounted onto a stereotaxic frame. A small craniotomy was performed on the tumour-ipsilateral right side of the skull 1.7 mm lateral to bregma and extending from 0 to 1.0 mm anterior of bregma. A 25 G needle with the bore facing to the right was introduced into the brain through the craniotomy to a depth of 2.5 mm. The needle was moved anterior and posterior three times across a 1.0 mm distance to sever the axons of the corpus callosum. Sham mice underwent the same analgesia and anaesthesia protocol but were not mounted onto the stereotaxic frame. Then, 2 weeks after injury, the mice were given an intraperitoneal injection of EdU (5 mg per kg) 2 h before brains were collected following transcardial perfusion with 4% PFA under terminal anaesthesia. Brains were fixed overnight in 4% PFA at 4 °C, transferred to PBS before vibratome sectioning (50-μm sections) and stored in cryobuffer (ethylene glcol:glycerol:PBS 1:1:2). Survival studies were carried out on tumour-bearing WT mice, which underwent brain injury or sham surgery at 8.5 weeks after electroporation with piggyBase/npp piggyBac plasmid mix as described above. Mice were euthanized when they showed any of the humane endpoints described above.
Behavioural assessment
Behavioural neuroscores were determined as follows. A ledge test was carried out observing mice walking along the edge of a cage and lowering themselves into the cage and scored as follows: 0, confident walk and good landing; 1, trips and wobbles while walking; 2, trips and wobbles, slips from ledge but recovers; 3, unable to walk along ledge. A hindlimb clasping test was carried out and scored as follows: 0, hindlimbs consistently pointing outward away from abdomen; 1, hindlimbs pulled in slightly towards body for more than 50% of the time; 2, hindlimbs pointed downwards towards abdomen for more than 50% of the time; 3, hindlimbs entirely retracted and touching the abdomen for more than 50% of the time. A gait test was carried out and scored as follows: 0, mouse moves normally; 1, slight tremor observed, slightly raised pelvis or slight waddle; 2, severe tremor, raised pelvis or pronounced waddle; 3, movements disjointed, stuttering with raised pelvis and severe waddle. A Kyphosis test was carried out and scored as follows: 0, easily able to straighten its spine as it walks; 1, mild kyphosis (curvature of the spine) but mostly able to straighten itself as it walks; 2, unable to straighten spine completely and maintains mild but persistent kyphosis; 3, maintains pronounced kyphosis as it walks or while it sits. Scores from each test were combined to give an overall neuroscore between 0 and 12. Mice were tested at 8 weeks for early disease stages and at 14 weeks and 16 weeks for advanced disease stages for WT and Sarm1−/− mice, respectively. At the late disease stage, mice were tested three times over three different days and neuroscores were averaged to account for the higher variability at advanced disease stage.
Visium ST data generation
For ST, brains from NSG mice with early or terminal PDX tumours or control NSG mice (15 weeks of age) were dissected, snap-frozen in methylbutane cooled to −20 °C in a bath of dry ice and liquid nitrogen, stored at −80 °C before embedding in OCT and sectioning at 10 μm on the Leica cryostat at −13 °C. From each brain, two 10 μm sections were collected at approximately 200 μm intervals from anterior to posterior in the striatal region onto 10x Visium Spatial Gene Expression slide (1000184, 10x Genomics). Slides were processed according to manufacturer’s instructions (10x Genomics) using a tissue permeabilization time of 36 min. RNA libraries were prepared according to the manufacturer’s instructions (library preparation kit 10x Genomics) and sequenced on the NovaSeq system with paired-end 150 bp reads.
Visium ST data analysis
Read selection and mapping
Reads from PDX experiments were aligned to reference genomes GRCh38-2020-A (human) and mm10-2020-A (mouse) using 10x Genomics Space Ranger v.2.0.1 (Supplementary Table 2). To assign reads to either species confidently, reads were mapped three times: (1) to the human genome; (2) to the mouse genome; (3) to a combination of both genomes. Reads mapping consistently to a single species using this procedure were selected and remapped to the combined genome as in point (3) to form the final dataset (Extended Data Fig. 2b). Reads from normal NSG mouse brains were filtered as above and remapped to the mouse genome for downstream analysis.
Selection of tumour-free spots
To define tumour-free Visium spots, we remapped reads from normal mouse brain to the combined genome and calculated for each spot the ratio between UMI counts assigned confidently to human or mouse. We reasoned that, as only mouse reads were present in the dataset, the distribution of these ratios was representing the background distribution (Extended Data Fig. 2c,d). On the basis of these data, we defined tumour-free spots as having a human to mouse ratio <2−4 (Extended Data Fig. 2c,d).
Data filtering and normalization
Spots with total UMI counts below 256 were discarded. Genes with non-zero counts in at least 1% of spots were retained for further analysis. Across tumour spots in each PDX model, human and mouse UMI counts were normalized separately using the posterior mean derived from the bayNorm package with parameter ‘mean_version=TRUE’53. bayNorm normalized counts were used in Fig. 2e and Extended Data Fig. 2g,h–o. For annotation of anatomical regions (Extended Data Fig. 2p), data were normalized using the sctransform package54.
Quantification of tumour density
Tumour density in individual Visium spots was measured using two different methods (Extended Data Fig. 2e,f and Supplementary Table 3): (1) using Visium spots positions on H&E images (Extended Data Fig. 2e (bottom)). For each spot, the area occupied by nuclei was calculated using Squidpy55. The spot area occupied by nuclei divided by the total spot area was then used as a measure of tumour density. (2) Based on sequencing data. For each spot, the total UMI counts from human transcripts divided by the total UMI counts from both species was used as a measure of tumour density.
Annotation of anatomical regions
To generate a reference dataset, normal mouse brain data were first clustered using the BayesSpace package56. Spatial clusters were then compared with annotation from the Allen Mouse Brain Atlas57 database and manually assigned labels (from specific anatomical regions). This reference dataset was used to annotate PDX data where brain morphology is partially disrupted by tumour cells and difficult to annotate manually. Data from each PDX was integrated with normal mouse brain data using the harmony package58. PDX anatomical regions were then predicted using random forest within the Harmony space59. Specifically, each dataset (PDX and normal mouse brain) was normalized using the sctransform package54. PDX and normal mouse brain data were then merged after selection of common features using the SelectIntegrationFeatures function from Seurat60, using RunHarmony in the Harmony package58. This was done in principal component analysis (PCA) space generated using the RunPCA function from Seurat60. For annotation, a random forest model was trained using annotated anatomic regions from normal mouse brain and its low-dimension space data from Harmony and used to predict anatomic regions in PDX data (Fig. 2e, Extended Data Fig. 2p and Supplementary Table 3).
Annotation of myelin high/low spots
We selected five myelination-related genes as markers of myelinated regions (Mbp, Cnp, Plp1, Mog and Mag)9 (Fig. 2a). Visium spots with mean total normalized counts for the five genes over the 70th percentile of their mean expression across spots in each section were labelled as myelinhigh.
Generation of pseudospots
The number of UMIs of human and mouse genes in a spot changes with tumour cell density (by definition) (Extended Data Fig. 2h–o). To control for spurious gene enrichments resulting from variations in number of a species UMI per spots (sequencing depth), a control dataset of randomized pseudospots was constructed as follows. In each PDX cell line, among spots with human/mouse ratio of total UMI counts between 0.5 to 1.5, the spot with the largest number of genes with non-zero UMI counts was selected for further downsampling and creation of pseudospots. For mouse genes, binomial downsampling53 was used on the spot UMI counts to generate pseudospots with UMI numbers corresponding to tumour densities ranging from 0.05 to 0.95. In total, 500 pseudospots were created in each PDX. As the pseudospots were created from a single spot, gene signatures are expected to show no correlation with tumour density.
Deconvolution of Visium data
To estimate cell type composition in each spot, Visium data were deconvoluted using the cell2location package61 and normal mouse brain scRNA-seq data from the Ximerakis study as reference 51 (Fig. 2d).
Normalization of cell type distributions across density bins
Tumour density was first discretized into 20 bins ranging from 0 to 1 with step size 0.05. Then, let xgcij denote the estimated number of cell types g in the ith spot of cth cell line, which lies in jth density bin, where j ∈ {1, …, 20}. In each spot, the proportion of cell types was calculated as \({p}_{gcij}={x}_{gcij}/{\sum }_{g}{x}_{gcij}\). Values in each bin were then summarized by taking the average across spots: \({\bar{x}}_{gcj}={\sum }_{i}{p}_{gcij}/{n}_{ci}\), where nci stands for the number of spots from the cth cell line in the jth bin. For each cell type in each bin, the average was taken across cell lines, and z-score-normalized across bins (the values are shown in Fig. 2f).
Gene signatures selection and enrichment analysis
The AUCell package was used to calculate gene signatures enrichment (area under the curve (AUC value))62. GO term gene lists were retrieved from the msigdbr database using the R package msigdbr63. Mouse GO terms with at least 50 genes were retained for further analysis (n = 1,980). Gene signatures used in Fig. 2b,c,e,g, are described in Supplementary Tables 4 and 6.
Human WM markers were derived using the ST dataset published previously26. In brief, spots from the cortex of samples 242_C, 248_C, 259_C, 265_C, 313_C and 334_C with log2-transformed total UMI counts between 8 and 14 were selected. Spots were combined and total counts were normalized. WM markers were defined as genes significantly that were upregulated in spots annotated as ‘white matter’ compared with spots annotated as ‘vascular’, ‘hyper cellular’, ‘grey matter’, ‘infiltrative’ and ‘necrotic edge’ (adjusted PWilcoxon < 0.01 and AUC value above 0.99 quantile of fitted Gaussian distribution on the AUC values reported from the wilcoxauc function of the R package presto)64.
These human WM markers and gene lists described in Supplementary Table 7 were used in Fig. 2g and Extended Data Fig. 2q.
Comparison of gene expression in WM and GM
Data from sections 1 and 2 of all ST experiments were divided into three groups: (1) normal healthy brain (NSG); (2) early tumours; and (3) terminal tumours (Fig. 2b). Count data from either WM or GM spots within each group were summed up to create pseudobulk RNA-seq datasets. Differentially expressed genes between WM and GM were then identified using DESeq2 (ref. 65) within each group separately.
Differentially expressed genes were selected using Padj < 0.01 and absolute log2-transformed fold change of >0.5 as cut-offs. Differentially expressed genes from the three groups were pooled and k-means clustering was performed on the log2-transformed fold change values with number of clusters set to 6. We applied the enricher function from the R package clusterProfiler66 on each cluster for enrichment analysis. Six GO terms enriched in cluster 2 were selected and the mean log2-transformed fold change of genes associated with each one of them is shown on Fig. 2b (right) (Supplementary Table 4 (geneID column)).
Comparison of WM spots between groups
The WM pseudobulk data from Fig. 2b for each group were used individually as input for DESeq2 (ref. 65) using the other two groups as a reference (Fig. 2c). The R package fgsea67 was used for GO term enrichment analysis on the log2-transformed fold change values from each group. NESs of selected significant GO terms (Padj < 0.1) are shown on Fig. 2c and Supplementary Table 6.
Normalization of gene signature across tumour density bins
Tumour density was first discretized into 20 bins ranging from 0 to 1 with step size 0.05 (Fig. 2e,g and Extended Data Fig. 2). Then, for each gene signature (g), let xgcij denote the AUC value of that gene signature from the ith spot of the cth cell line in the jth bin, where j ∈ {1, …, 20}. The average AUC value of the spots in the jth bin: \({\bar{x}}_{gcj}={\sum }_{i}{x}_{gcij}/{n}_{cj}\), where ncj stands for the number of spots from the cth cell line in the jth bin. Then the mean across cell lines was calculated and as \({\bar{x}}_{gj}={\sum }_{c}{\bar{x}}_{gcj}/{n}_{c}\), where nc stands for the number of cell lines used. Finally, \({\bar{x}}_{gj}\) was z-score normalized across 20 bins such that \({z}_{gj}=\frac{{\bar{x}}_{gj}-{\mu }_{gj}}{{\sigma }_{gj}}\), where μgj and σgj stand for the mean and s.d. of \({\bar{x}}_{gj}\) across 20 bins respectively.
Comprehensive evaluation of gene signatures expression trends as a function of tumour density
A Mann–Kendall trend test (R function mk.test from the R package trend)68,69, which was originally developed for testing monotonic trend in time-series data, was used to explore expression trends of gene signatures as a function of binned tumour densities. Let \({x}_{gcj}={{\rm{median}}}_{i}({x}_{gcij})\) denotes the median of AUC value of gene signature (g) of cth cell line in the jth bin (xgcij as defined above). mk.test with alternative=two.sided was applied to xgc across 20 bins (bins with missing values due to a limited number of spots were not considered) for each cell line and pseudospots. mk.test reports two statistics, S and pval. Positive/negative S values stand for increasing/decreasing trend of gene signature as a function of binned tumour densities (\(S=\mathop{\sum }\limits_{k=1}^{n-1}\mathop{\sum }\limits_{j=k+1}^{n}{\rm{sgn}}({x}_{gcj}-{x}_{gck})\) where sgn is the sign function and n = 20 is the number of bins), while pval indicates whether that trend is significant or not68,69. GO terms with at least one PDX cell line with pval < 0.1 were kept for k-means clustering on S values (Extended Data Fig. 2g and Supplementary Table 5).
Reanalysis of published spatial datasets from human glioblastoma
For the ref. 26 dataset, data were downloaded from https://datadryad.org/stash/dataset/doi:10.5061/dryad.h70rxwdmj. Tumour densities were determined from H&E images using the image based approach used on Extended Data Fig. 2e,f (see above).
For the ref. 27 dataset, Cosmx data were downloaded from https://data.mendeley.com/datasets/wc8tmdmsxm/3. Following the preprocessing steps reported previously27, cells with fewer than 20 total transcripts, fewer than 20 genes detected or more than 3 negative control probes were removed. Filtered data were log-normalized and scaled using Seurat. Clustering of cells was done using PCA space for identifying tumour cells based on marker genes from refs. 7,27. For each cell (including tumour and non-tumour cells), we calculated the proportion of tumour cells present around it within a 55 µm diameter circular area (corresponding to the spot size on the Visium platform). These tumour densities were then discretized into 20 bins. Bins with upper bounds 0.05, 0.8, 0.85, 0.9, 0.95 and 1 were discarded as the number of cells per bin was low (<300).
Tissue preparation and immunohistochemistry
Animals were perfused (4% PFA in PBS; Merck P6148) under terminal anaesthesia, brains were collected, post-fixed overnight at 4 °C in PFA (4%) before transferring to PBS. Vibratome sections (50 µm) were prepared and stored in cryopreservative (glycerol:ethylene glycol; PBS 1:1:2) before immunohistochemistry. For staining, floating sections were permeabilized overnight (1% Triton X-100, 10% serum in PBS) at 4 °C, incubated in primary antibodies overnight (1% Triton X-100, 10% serum in PBS) at 4 °C and for 3 h in secondary antibody (0.5% Triton X-100, 10% serum in PBS) containing DAPI counterstain (Insight Biotechnology, sc3598). The sections were mounted with antifade mounting solution (Prolong gold antifade mountant, Thermo Fisher Scientific, P36934) before imaging on a 3i confocal spinning disk (3i SlideBook Version 2023). For imaging of axonal damage, brain tissue from Thy1-YFP mice were imaged using the Airyscan function of the LSM 880 confocal microscope (Zeiss Zen Black v.2.1).
The following antibodies were used: rabbit anti-Ki-67 (1:250; Abcam, ab16667), goat anti-GFAP (1:1,000, Abcam, ab53554), rat anti-CD68 (1:500, Abcam, ab53444), rabbit anti-Iba1 (1:1,000, Wako, 019-19741, L0159), mouse anti-neurofilament H (1:1,000, Enzo, ENZ-ABS219-0100), mouse anti-MBP (1:1,000, Covance, SMI-99), mouse anti-SMI32 (1:1,000, Enzo, ENZ-ABS219-0010) chicken anti-GFP (1:1,000, Abcam, ab13970), rabbit anti-RFP (1:1,000, ABIN129578), rabbit anti-pMLC2 (1:100, Cell Signalling, 3671), mouse anti-phospho-Tau S202/T205 (1:500, a gift from G. Schiavo), mouse anti-TDP-43 (1:500, Abcam, ab104223), rabbit anti-TOMM20 (1:1,000, ab186735) and mouse anti-amyloid-β (1:100, Merk, MAB348A4), rabbit anti-laminin (1:500, Sigma-Aldrich, L9393), goat anti-CD31 (1:100, BioTechne, AF3628), rat anti-PDGFRB (1:200, gift from I. Kim), donkey anti-mouse IgG 488 (1:500, Thermo Fisher Scientific, A21202). For detection of EdU, the sections were stained using the Click-it EdU Alexa Fluor 647 Imaging Kit (Invitrogen, C10340) according to the manufacturer’s guidelines.
Hypoxic regions were identified by intraperitoneal injection of pimonizadole (60 mg per kg; Hypoxyprobe Omnit Kit HP3-1000Kit) 90 min before brains were collected after transcardial perfusion with 4% PFA under terminal anaesthesia. Brains were fixed in 4% PFA overnight and sectioned (40 μm) on a vibratome. The sections were permeabilized in 0.3% Triton X-100, 10% donkey serum in PBS overnight before incubation in primary antibody (1:1,000 rabbit anti-pimonidazole adducts) and detection with donkey anti-rabbit Alexa Fluor 647 (1:1,000, Thermo Fisher Scientific, A-31573) and counterstained with DAPI.
Computational image analysis
Analysis of tumour cell localization and proliferation was performed in Imaris 10.1.0 on single z plane images from a 3i spinning-disk microscope. Spot segmentation was first performed on tdTomato/GFP channel, before being filtered for intensity median or centre on DAPI to segment tumour cells. Tumour cells were then classified as EdU/Ki-67+/−. For quantitative assessment within WM and GM (Fig. 1b,c and Extended Data Fig. 1b,g), we analysed the striatum because it is an anatomically well-defined brain region that is infiltrated by early tumour cells and contains both GM and WM in discrete bundles. Surfaces were manually drawn for the SVZ, haemorrhagic/necrotic regions, striatum and injury sites, and tumour cells within SVZ and haemorrhagic/necrotic regions were filtered out. WM bundle surfaces were generated using the machine learning function. The percentage of WM area (Extended Data Fig. 1b) was calculated by dividing the area of WM bundles in tumour infiltrated striatum by the total area of tumour infiltrated striatum. For analysis in Extended Data Fig. 6i,k, tdTomato spots were additionally filtered on a surface generated for the virally targeted area using GFP fluorescence.
Analysis of Thy1-YFP (Fig. 3b and Extended Data Fig. 3a,c) and neurofilament (Extended Data Fig. 3e) mean fluorescence intensity, as well as GFAP+ cell density (Extended Data Fig. 4e) and CD68 integrated density (Extended Data Fig. 4g) was performed in ImageJ on maximum-intensity projection (MIP) images from a 3i spinning-disk confocal microscope using a custom script. Individual bundle ROIs were manually drawn and tdTomato+ and GFAP+ cells manually counted. ROI area and mean fluorescence intensity was measured using the Measure function. Mean fluorescence intensity was normalized to the average of mean fluorescence intensities in contralateral bundles (≥5 bundles per animal). CD68 integrated density was measured by first thresholding CD68 channel with Li autothreshold, and integrated density (IntDen) was quantified using the AnalyzeParticles function. Analysis of distance of axonal varicosities to a tumour cell body or tumour cell process (Fig. 3f) was performed in ImageJ on single z-plane images. Individual varicosities (n = 111) within tumour-involved WM were manually selected, and the distances between varicosities and tumour cells were measured using the Measure function. Varicosities located within a distance of <5 µm from the tumour cell body or cell process were categorized accordingly, while those at a distance of >5 µm classified as ‘No tumour cell’.
Analysis of GFAP area (Fig. 4e and Extended Data Figs. 5e and 6e) and CD68 intensity (Fig. 4f and Extended Data Figs. 5f and 6f) was performed in ImageJ on MIP images from the 3i spinning-disk confocal microscope using a custom script. Triangle threshold was used on the tdTomato image to generate a tdTomato ROI, which was used for ‘sham’. ROIs for the injury site were manually drawn in ImageJ. The injury site ROI was generated from the overlap of tdTomato and injury site ROIs. The ‘injury (excluding injury site)’ ROI was generated from the tdTomato ROI excluding the injury site ROI. The injury site was excluded from the analysis to avoid confounding effects of elevated neuroinflammation in this region after wounding. The GFAP area was calculated by thresholding the GFAP channel using the ImageJ Triangle threshold, and measuring the area covered within each ROI using the ImageJ AnalyzeParticles function. CD68 analysis was calculated by thresholding CD68 channel with Triangle or Li autothreshold, and integrated density (IntDen) was quantified using the AnalyzeParticles function. For both, measurements were normalized to their own sham control at each timepoint. Analysis of injury responses in Fig. 4 and Extended Data Figs. 5 and 6 was carried out across all areas occupied by tdTomato+ tumour cells.
Analysis of GFAP area and CD68 intensity for time course in npp WT mice (Extended Data Fig. 4b,c) and PDX (Extended Data Fig. 4m,n) was performed as above for Sham mice, and measurements normalized to control (non-tumour-bearing brains) or contralateral, respectively. Analysis of vascular phenotypes (Extended Data Fig. 8e–h,j,l) was performed in ImageJ on MIP images from the 3i spinning-disk confocal microscope. Images were converted into RGB images, and the Vessel Analysis plug-in was used to produce thresholded vasculature images and derive the percentage of CD31+ area, vascular length (measured as the vascular length density) and the mean vascular diameter. Furthermore, the Skeletonize3D and AnalyzeSkeleton plugins were used on the thresholded images produced by the complete vessel analysis to derive number of branches. Colocalization of laminin or PDGFRB with CD31 was measured using thresholded images of CD31 and either laminin or PDGFRB, colocalization was determined using the Image Calculator AND function, and the percentage colocalization was calculated using the Measure Area function on CD31 and CD31 and laminin or CD31 and PDGFR, respectively. Analysis of IgG area (Extended Data Fig. 8n) was performed in ImageJ on MIP images from the 3i spinning-disk confocal microscope. The Triangle threshold was used on the tdTomato image to generate a tdTomato tumour ROI. IgG-positive areas were drawn manually and the area covered within each tdTomato ROI was measured using the ImageJ AnalyzeParticles function. To produce the rendered image in Fig. 3g, the z-stack confocal microscopy image was imported, 3D reconstructed and processed in Imaris v.10.1.0. 3D. Surfaces were constructed from the tdTomato/GFP channels using the Surfaces segmentation tool, combining automatic and manual segmentation. This resulted in 3D surfaces representing co-localized tumour cells and axon with axonal varicosities, which were subsequently animated in 3D alongside the three-channel tdTomato/GFP/mitoBFP confocal microscopy images (Supplementary Video 1).
Image quantification of H&E images
For H&E image quantification, the watershed method was applied to grey scale smoothed images using the Python package squidpy55 for segmentation of nuclei. The function skimage.measure.regionprops_table from the Python package scikit-image70 was used to count the number of cells.
Atomic-force microscopy
Thy1-YFP brains bearing intermediate npp tumours were snap-frozen in liquid nitrogen before sectioning at 10 μm on the Leica cryostat. Atomic-force microscopy measurements were performed using an MFP-3D BIO Inverted optical atomic-force microscopy (Asylum Research) mounted on a Nikon TE2000-U inverted fluorescence microscope and placed onto a vibration-isolation table (Herzan TS-150). Silicon nitride cantilevers with a nominal spring constant of 0.06 N m−1 and a borosilicate glass spherical tip with 5 μm diameter (Novascan Tech) were used. Cantilevers were calibrated using the thermal fluctuation method. Frozen sections were equilibrated to room temperature by immersion in PBS for 5 min before mounting. TdTomato and Thy1-YFP fluorescence were identified in the same section, the cantilever was placed in the corresponding regions and the specimens were indented at a 2 μm s−1 loading rate. The Young’s moduli of the samples were determined by fitting force curves with the Hertz model using a Poisson ratio of 0.5.
Single-cell RNA preparation
Mouse brains were collected into ice-cold HBSS medium and dissected into 1 mm coronal sections using a brain matrix (WPI, RBMS200C). Tumour regions were dissected out and mechanically dissociated into small pieces. Cells were isolated by papain dissociation (as above) and RNA libraries prepared using Chromium Next GEM Chip G Single Cell Kit (10x genomics; 1000127) and sequenced on Nova Seq X Plus PE 150.
scRNA-seq data analysis
Read selection and mapping
Reads were preprocessed and mapped to the mm10-2020-A mouse genome using 10x Genomics Cell Ranger v.7.0.1 (Supplementary Table 8 and 9). The tdTomato sequence, expressed by transformed cells, was added to the reference genome.
Cells and genes filtering
Cells with zero UMI counts for the 4 red blood cell markers Hbb-bs, Hba-a1, Hba-a2 and Hbb-bt and with either tdTomato expression ≤ 2 (microenvironment cells) or tdTomato ≥ 5 (tumour cells) were retained for further analysis. For all analysis, cells from WT mice were downsampled so that each genotype had ~20,000 cells. Cells with a proportion of mitochondrial genes of below 0.25 and log2-transformed total counts of between 9 and 16 were retained for further analysis. Genes with non-zero UMI counts in at least 0.5% of cells were retained for further analysis. As a result, the dataset presented in this study consists of 19,939 and 21,206 cells for the WT and Sarm1−/− samples, respectively, and of a total of 14,842 genes.
Identification of high-confidence tumour cells
Two rounds of data filtering were used to identify high-confidence tumour cells. First, the Harmony package58 was used to integrate the WT and Sarm1−/− datasets from this study with scRNA-seq datasets from normal mouse brain51,71, TAMs72 and from a mouse GBM model, which contains annotated tumour cells23. Specifically, data from each study were normalized using the sctransform package54 and merged after selection of common features using the SelectIntegrationFeatures function from Seurat60. Then, the batch correction function RunHarmony from the Harmony package58 was applied to the data in PCA space (generated with the RunPCA function in Seurat)60.
High-confidence tumour cells were defined as either (1) expressing at least 5 tdTomato UMI count; (2) predicted to be aneuploid using the copyKat package73; (3) predicted to be tumour cells using the integrated dataset annotation from ref. 58 and a random-forest approach in Harmony space. Specifically, labels from refs. 23,51,71,72 were used for training a random-forest model59, which was then applied to predict cell labels in the scRNA-seq data from this study. Finally, tumour cells with UMI counts for the Ptprc (Cd45) and Cd68 genes of >0 were discarded as these are considered to be immune-cell-specific markers.
In a second round of filtering, tumour cells were again integrated and clustered using the Harmony package58 and the Louvain approach74 both in Harmony space (same procedure as above, but this time the integration was done using scRNA-seq from this study only). Cells identified to be TAMs or endothelial cells based on markers from refs. 51,71 were removed from the high-confidence list.
Identification of high-confidence non-tumour cells
Cells with tdTomato UMI counts ≤ 2 and predicted to be diploid using the copyKat package73 were called high confidence.
Cell type annotation of high-confidence tumour and non-tumour cells
High-confidence tumour and non-tumour cells were integrated (between genotypes) and clustered separately using Harmony and Seurat58,60. Clustering was performed using the Louvain approach in Harmony space with resolution = 0.2 (ref. 74). Clusters were annotated using lineage markers and the gene enrichment analysis package fgsea67. Cluster annotation was finally checked manually for accuracy.
Cell type annotation of unassigned cells
Unassigned cells are cells that are not part of the two high-confidence lists. First, all the cells from this study were integrated using Harmony as above58. Second, a random-forest model was trained using high-confidence tumour and non-tumour cells (training dataset) and used to identify and annotate tumour cells in the list of unassigned cells. Third, a new round of clustering was applied on tumour or non-tumour cells separately. Cell type labels were then assigned using random forest and the cell type annotation from high-confidence tumour or non-tumour cells (Supplementary Table 9). TAMs were reclustered separately and macrophage/microglial markers from two studies44,45 were used for annotation (Extended Data Fig. 9a).
Proportion test
To perform proportion tests on equal numbers of cells in both genotypes tumour and non-tumour cells, were downsampled to 12,054 tumour cells and 4,000 cells respectively. The prop.test function from R was used to test the significance of difference in proportion of cell types between two the genotypes. P < 0.05 and absolute proportion difference above 0.1 was considered to be significant (Fig. 5g,h and Extended Data Fig. 9d).
Differential gene expression and gene enrichment analysis
Differentially expressed genes between Sarm1−/− and WT cells in each cell type (PWilcoxon < 0.01 and AUC value above the 0.99 quantile of the fitted Gaussian distribution on the AUC values reported by the wilcoxauc function of the R package presto)64 were analysed for GO enrichment using the enricher function of the R package clusterProfiler66,75 (Supplementary Table 10).
Ligand–receptor analysis
The cellphoneDB method76 from the LIANA package77 was used to identify significant ligand–receptor pairs in each genotype (Extended Data Fig. 9e–h).
Flow cytometry analysis
Brains were collected into ice-cold HBSS medium and dissected into 1 mm coronal sections using a brain matrix (World Precision Instruments, RBMS200C). Tumour regions were dissected out and mechanically dissociated into small pieces, followed by enzymatic dissociation using Liberase TL (Roche, 05401119001) supplemented with DNase I (Merck, 11284932001) for 30 min at 37 °C. After addition of EDTA to stop the enzymatic reaction, cells were washed with PBS and filtered through a 70 μm cell strainer (Falcon, 352350) to remove large debris. The samples were blocked on ice for 20 min (BioXCell blocking buffer, BE0307) before incubation in antibodies and fixable viability dye eFluor780 (eBioscience, 65-0865-18, 1:1,000) at 4 °C for 20 min. To detect immune cells within the tumour population, the following antibodies were used rat anti-LY6G-BUV563 (1:100, IA8, BD, 612921), rat anti-CD11b-BUV661 (1:400, M1/70, BD, 612977), rat anti-MHC Class II-BB700 (1:800, M5/114.15.2, BD, 746197), mouse anti-CD45-BUV805 (1:400, 30-F11, BD, 748370), mouse anti-CD64-BV421 (1:100, X54-5/7.1, BioLegend, 139309), mouse anti-CX3CR1-BV510 (1:400, SA011f11, BioLegend, 139309), rat anti-LY6C-BV605 (1:200, AL-21, BD, 563011), rat anti-CD19-BV650 (1:50, ID3, BD, 563235), hamster anti-CD11C-BV785 (1:100, N418, BioLegend 117336), rat anti-CD49d-APC (1:200, R1-2, BioLegend, 103622), rat anti-F4/80-AF700 (1:100, BM8, BioLegend, 123130), mouse anti-Ki67-BUV395 (1:100, B56, BD, 564071), rat anti-CD3-BUV737 (1:300, 17A2, BD564380), rat anti-CD206-AF488 (1:100, C068C2, BioLegend, 141710). Flow data were acquired using BD FACSymphony, FACS DIVA version 9.1. Data were analysed using BD FlowJo Software (v.10.8.1). Data were compensated (using ArC reactive and negative beads (Invitrogen, A10346 A and B) for viability dye, and UltraComp eBeads Compensation Beads (Invitrogen, 01-2222-42) for all other fluorophores), fluorescence minus one controls were generated, and only viable singlets were used for downstream analysis.
Targeted EM
Brains were perfused with electron microscopy (EM)-grade 4% formaldehyde immersion fixed overnight, embedded in 4% agarose and sectioned on a vibrating microtome (100 µm).
Sections were stained with DAPI and imaged using confocal microscopy (×20 objective) to map the tdTomato+ tumour cells and identify regions of interest. These regions were prepared for electron microscopy by processing, ultrathin sectioning and imaging on a scanning electron microscope (SEM)78,79. All EM analysis was conducted on ≥50 axons per bundle in ≥3 bundles (n = 4 mice). Degenerating axons were identified as those exhibiting any of the following features of axonal pathology: condensed/dark axoplasm, organelle accumulation, axonal swelling, vacuolization (Fig. 3c,d). For quantitative analysis of demyelination, inner diameter, outer diameter, myelin thickness and corresponding g-ratios of myelinated axons were semiautomatically calculated using the software program MyelTracer (v.1.3.1)80. Feret diameters were used to account for the imperfect circularity of axons9.
Statistical analysis and data visualization
Statistical analysis was performed in Prism 10 or R (v.4.3.2). Significance was calculated as indicated in the figure legends. All t-tests were two-tailed. All data are expressed as mean ± s.d. unless otherwise stated. Exact P values are provided in the source data and on the figures or in the legends. No statistical method was used to predetermine sample size. Sample size was determined based on existing literature and our previous experience. Data visualization was done using the ggplot2 package in R81. Heat maps was generated using the ComplexHeatmap package82.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
