Animals
This study followed some applicable aspects of the PREPARE89 planning guidelines checklist, such as the formulation of the in vivo study, dialogue between scientists and the animal facility, and quality control of the in vivo components in the study. All animals were born, bred and housed at the Karolinska Institutet, Comparative Medicine Biomedicum animal facility (KMB). Mouse brain tissues (postnatal days P0, P2, P5, P7, P10 and P21) were obtained from a mouse line generated by crossing Sox10:cre animals (The Jackson Laboratory, 025807) on the C57BL/6j genetic background with RCE:loxP (eGFP) animals (The Jackson Laboratory, 32037-JAX) on a C57BL/6xCD1 mixed genetic background. Females with a hemizygous cre allele were mated with males lacking the cre allele, while the reporter allele was kept in hemizygosity in both females and males. In the resulting Sox10:cre-RCE:loxP (eGFP) progeny, the entire OL lineage was labelled with eGFP.
None of the experimental animals in this study were subjected to previous procedures before enrolment in the study. All of the animals were free from mouse viral pathogens, ectoparasites, endoparasites and mouse bacterial pathogens. Mice were kept under the following light–dark cycle: dawn, 6:00–7:00; daylight, 7:00–18:00; dusk, 18:00–19:00; night, 19:00–6:00; a maximum number of 5 mice were housed per cage in individually ventilated cages (IVC Sealsafe plus GM500, Tecniplast). General housing parameters such as relative humidity, temperature and ventilation follow the European Convention for the Protection of Vertebrate Animals used for experimental and other scientific purposes treaty ETS 123, Strasbourg 18.03.1996/01.01.1991. In brief, consistent relative air humidity of 55 ± 10% was maintained at 22 °C and the air quality was controlled with the use of stand-alone air handling units supplemented with HEPA filtrated air. Monitoring of husbandry parameters was done using ScanClime (Scanbur) units. Cages contained hardwood bedding (TAPVEI), nesting material, shredded paper, gnawing sticks and card box shelter (Scanbur). The mice received a regular chow diet (CRM(P) SDS and CRM(P), SAFE). Water was provided by using a water bottle, which was changed weekly. The cages were changed every other week. Cage changes were done in a laminar air-flow cabinet (NinoSafe MCCU mobile cage changing unit) with a HEPA H14 EN1822 filter (0.3 μm particle size). Facility personnel wore dedicated scrubs, socks and shoes. Respiratory masks were used when working outside of the laminar air-flow cabinet. Both sexes were included in the study. Randomization was performed to assign samples to time-point groups.
All experimental procedures on animals were performed following the European directive 2010/63/EU, local Swedish directive L150/SJVFS/2019:9, Saknr L150, Karolinska Institutet complementary guidelines for procurement and use of laboratory animals, Dnr. 1937/03-640 and Karolinska Institutet Comparative Medicine veterinary guidelines and plans (version 2020/12/18). The procedures described were approved by the regional committee for ethical experiments on laboratory animals in Sweden (Stockholm Norra Djurförsöksetiska Nämnd, 1995/2019 and 7029/2020).
Tissue and slide preparation
For all postnatal collection points except for P21, we performed a decapitation and extracted the brain, immediately placing the tissue in Tissue Tek O.C.T. compound over a bath of dry ice and 70% ethanol. For all mice older than P21 we anaesthetized the mice and then performed a transcardiac perfusion with gassed, ice-cold artificial cerebral spinal fluid. Next, we decapitated the mice and removed the brain and followed the same embedding procedure as for the neonates.
All tissues were stored at −80 °C until further usage. Fresh 50 × 75 mm (Ted Pella, NC1811932) 0.01% poly-L-lysine-coated (Sigma-Aldrich, P1524-100) slides were prepared and stored at 4 °C and used within 1 week for mounting tissue. For brain tissue sectioning, the tissue was incubated for 30 min in the cryostat chamber at −20 °C. The slides were simultaneously precooled (poly-L-lysine coated for DBiT and charged Superfrost Plus (Eprendia, 4951PLUS-001) for CODEX) and kept in the cryostat for the entire procedure. The sections were cut using an antiroll plate, to a thickness of 10 µm. Tissue collection was done by placing the cold slide in the correct position on top of the tissue then placing over a gloved backhand until the tissue was fully bound. The slide was immediately placed back in the cryostat chamber until sectioning was complete, at which point all slides were stored at −80 °C.
LPC mouse model (1% LPC injection)
Sox10:cre-RCE:loxP (eGFP) mice of mixed sex, at around P85, were used in the LPC study. In brief, under an aseptic technique, mice received an intraperitoneal injection of a non-steroidal anti-inflammatory analgesic (Carprofen; 5 mg kg−1), and were then deeply anaesthetized with isoflurane; an ophthalmic ointment was then applied in both eyes to prevent corneal desiccation, and the mouse was placed over a warm pad for the entire procedure to prevent hypothermia. The site of the incision was shaved and sanitized with 2% chlorhexidine before cutting a 0.5 cm incision with a scalpel to expose the cranium. Using a stereotaxic frame, the rough coordinates of −0.8 mm posterior and 0.8 mm lateral were measured for a skull drilling perforation. After completing, the 1% LPC solution was loaded into a pulled-glass micropipette and the coordinates were recalculated. The micropipette was then slowly inserted 1.3 mm deep (CC below cingulum bundle). The 2 µl injection was performed using an injection speed of 5 nl s−1. After the injection was complete the micropipette was kept in place for an additional 5 min to prevent efflux and then slowly pulled out of the brain. Monitoring of reflexes, respiratory rate and breathing pattern was performed during the entire surgical procedure. The skin was sutured with a non-absorbable suture (Vicryl, Ethicon), and the mouse was observed until regaining full consciousness and mobility. Post-operative care was administered daily for the next 72 h. Mice were euthanized at 5, 10 and 21 days after surgery using the sampling procedure and sample collection method as for the P21 mice (above).
Retro-AAV-eGFP injection
P1 mice were injected with 0.5 μl of retro-AAV-eGFP virus (ssAAV-retro/2-CAG-eGFP-WPRE-SV40p(A); V24-retro; produced at the Universität Zürich Viral Vector Facility) into either layer II/III of the CTX (1.5 mm anterior of lambda, 0.7 mm lateral of midline and 0.5 mm deep from skin) or the thalamus (1.5 mm anterior of lambda, 1.6 mm lateral of midline and 2.8 mm deep). In brief, neonates received a subcutaneous injection buprenorphine (0.5 mg kg−1) and, after 30 min, the neonates were put under general anaesthesia using isoflurane. The bregma coordinates were found, and a pulled-glass capillary preloaded with virus and mineral oil was inserted to the desired depth before injection of virus into the parenchyma at a rate of 100 nl min−1. After the injection, the glass capillary was kept in place for an additional 5 min before it was slowly withdrawn. The neonate’s foot was tattooed, and the neonates were then acclimatized back to the mother by rubbing bedding on the pup to prevent maternal rejection. The neonates were monitored to verify that they successfully came out of anaesthesia. Daily monitoring and pain assessments were performed until euthanization at P10. Mice were perfused with cold PBS followed by 4% PFA (3–5 ml each) and tissue was incubated overnight in 4% PFA and then transferred to 30% sucrose for 2–3 days. Finally, the brain tissue was cryosectioned at a thickness of 20 µm, stained with anti-rat MBP (EMD Millipore, MAB 386) antibodies and both eGFP and Alexa Fluor 555 were imaged using a Zeiss LSM900-Airy using z stacking. Blinded imaging analysis was performed by first creating a maximum-projection image for each channel; the MBP area was drawn for each image then the percentage area of MBP, retro-AAV-eGFP and the dual overlap percentage area was measured. To compare the extent of myelination for each labelled tract, we calculated the area under the curve across the entire CC (percentage area of the CC that was eGFP and MBP positive) then performed a two-tailed t-test using the total mean area of MBP and standard error values that were generated through the area under the curve calculation. Experimental blinding was implemented during retro-AAV-eGFP image analysis.
Tris-mediated retention of in situ hybridization signal during clearing DISCO
Three Sox10:cre-RCE:loxP (eGFP) mice (aged 8–10 weeks) were stereotactically injected with 1% LPC and euthanized 5 days after the procedure, as described above except that a 4% PFA perfusion was included directly after the artificial cerebrospinal fluid. Each brain was carefully extracted and then stored overnight in 4% PFA at 4 °C before processing through the full TRIC-DISCO protocol73. In brief, a methanol gradient (60%, 80%, 100%, 100%) was done and then stored overnight at −20 °C. Successive steps to delipidate (100% dichloromethane), wash and bleach (100% methanol + 5% hydrogen peroxide) were done with overnight incubations at 4 °C. For whole-brain in situ hybridization chain reaction, the samples were incubated with 50% formamide wash buffer until brain sinking, followed by incubation for 3 days at 37 °C in 50% hybridization buffer with Csf1r DNA probe (1.6 pmol of probe in 400 µl). The samples were then washed three times in 5× saline-sodium citrate (5× SSC) buffer with Tween-20 (5× SSCT), incubated in amplification buffer overnight (room temperature), and then washed three times in 5× SSCT buffer. Tris-mediated signal retention involved three wash steps (for 1 h at room temperature) in 500 mM of Tris-HCl, sample dehydration in 100% methanol (1 h at room temperature) then dehydration and delipidation in 66% DCM and 33% methanol (3 h) and finally two washes with 100% DCM. To promote tissue transparency an overnight dibenzl ether (DBE) incubation was performed with one DBE solution change. Whole-brain light-sheet microscopy using LaVision Ultramicroscope II (MiltenyBiotec) provided images that were used in the study.
Human samples and slide preparation
De-identified human tissue samples from the second trimester were collected at Zuckerberg San Francisco General Hospital (ZSFGH), with the acquisition approved by the UCSF Human Gamete, Embryo and Stem Cell Research Committee (10-05113). All of the procedures adhered to protocol guidelines, and informed consent was obtained before sample collection and use for this study. The de-identified third-trimester and early postnatal tissues without known neurological disorders were obtained from the University of Maryland Brain and Tissue Bank through NIH NeuroBioBank under the bank’s Institutional Review Boards oversight. For all samples, coronal cryosections of 10 μm were prepared on poly-L-lysine-coated glass slides (63478-AS, Electron Microscopy Sciences) or 50 × 75 mm poly-L-lysine-coated glass slides, and stored at −80 °C for subsequent use.
Multiplex immunofluorescence imaging (CODEX)
Immunofluorescence imaging was conducted using the Akoya Biosciences PhenoCycler system, which is integrated with a Keyence BZ-X800 epifluorescence microscope. The protocol for antibody–barcode conjugation, staining and imaging fresh-frozen tissue, as described in the manufacturer’s user manual, was strictly adhered to. Antibodies and detailed information are listed in Supplementary Tables 6 and 13. As CODEX images vary in magnification, the group’s post-acquisition exposure settings are uniform for each magnification within a figure. All post-acquisition processing was performed using QuPath90. Batch-to-batch (S1 and S2) variation correction for the exposure settings for each antibody channel was adjusted by performing mean intensity measurements over one hemisphere using ImageJ91 until the mean intensities aligned between the exposure settings.
DNA barcodes sequences, DNA oligos and other key reagents
DNA oligos used for transposome assembly, PCR and library construction are shown in Supplementary Table 1. All of the DNA barcode are provided in Supplementary Tables 2 and 3, and all other chemicals and reagents were listed in Supplementary Table 6. The hydrogel-based DBiT devices preloaded with 220 barcodes were purchased from AtlasXomics (AXO-0457).
Spatial ARP-seq
The fresh-frozen tissue slide was thawed for 10 min before fixing with 0.2% formaldehyde with 0.05 U μl−1 RNase inhibitor (RI) and quenching with 1.25 M glycine for 5 min. After three washes with 0.5× DPBS-RI, the tissue section was permeabilized with lysis buffer for 15 min, followed by another two washes with 0.5× DPBS-RI and a final wash with deionized H2O. The ADT cocktails (diluted 20 times from the original stock) from BioLegend and self-conjugated antibodies (including intracellular and surface proteins for mouse brain) were added and incubated for 3 h at 4 °C. The mouse ADT protein panel encompasses 119 antibodies targeting cell surface antigens, 9 isotype control antibodies specific to immune cells and, optionally, 8 self-conjugated antibodies for intracellular markers representative of canonical mouse brain cell types (Supplementary Table 4). Similarly, the human ADT panel includes 154 antibodies against unique cell surface antigens, covering key lineage antigens, complemented by nine isotype control antibodies for immune cell profiling (Supplementary Table 5). The ADT cocktail was removed by washing three times with 0.5× DPBS-RI. The tissue was then washed for 5 min with a lysis buffer containing 0.001% Digitonin, 10 mM Tris-HCl (pH 7.4), 3 mM MgCl2, 0.01% NP-40, 1% BSA, 10 mM NaCl, 0.05 U μl−1 RNase inhibitor and 0.01% Tween-20. Subsequently, it was washed twice for 5 min each with a wash buffer composed of 0.1% Tween-20, 1% BSA, 10 mM Tris-HCl (pH 7.4), 3 mM MgCl2 and 10 mM NaCl. A transposition mix was prepared consisting of 5 μl homemade transposome, 50 μl 2× Tagmentation buffer, 33 μl 1× DPBS, 1 μl 1% Digitonin, 0.05 U μl−1 RNase inhibitor, 10 μl nuclease-free H2O and 1 μl 10% Tween-20. This mixture was then added to the sample and incubated at 37 °C for 30 min. Next, 40 mM EDTA containing 0.05 U μl−1 RNase inhibitor was added and incubated for 5 min at room temperature to halt the transposition process. Subsequently, the tissue section was washed three times with 0.5× PBS-RI for 5 min each. The RT mixture consisting of 3.1 μl 10 mM dNTPs each, 4.5 μl RNase-free water, 0.4 μl RNase inhibitor, 12.5 μl 5× RT buffer, 6.2 μl Maxima H Minus reverse transcriptase, 0.8 μl SuperaseIn RNase inhibitor, 10 μl RT primer and 25 μl 0.5× PBS-RI was added to the tissue. The tissue was then incubated for 30 min at room temperature and subsequently at 42 °C for 90 min in a humidified chamber. After the reverse transcription reaction, the tissue was washed for 5 min with 1× NEBuffer 3.1 containing 1% RNase inhibitor.
For the first in situ ligation of barcodes (barcodes A), the PDMS chip was positioned over the region of interest (ROI) on the tissue slide. To ensure proper alignment, a bright-field image was captured using a ×10 objective on the Thermo Fisher Scientific EVOS FL Auto microscope (AMF7000) with EVOS FL Auto 2 Software Revision v.2.0.2094.0. The PDMS device and tissue slide were then securely fastened using a homemade clamp. Initially, barcode A was combined with ligation linker 1: 10 μl of 100 μM ligation linker, 10 μl of 100 μM each barcode A and 20 μl of 2× annealing buffer (2 mM EDTA, 100 mM NaCl, 20 mM Tris, pH 7.5–8.0) were thoroughly mixed. For each channel, a ligation master mixture was prepared, containing 2 μl of ligation mix (11 μl T4 DNA ligase, 72.4 μl RNase-free water, 27 μl T4 DNA ligase buffer, 5.4 μl 5% Triton X-100), 2 μl of 1× NEBuffer 3.1 and 1 μl of each annealed DNA barcode A (A1–100, 25 μM). The ligation master mixture was loaded into 100 channels of the device using a vacuum. The device was then incubated at 37 °C for 30 min in a humidified chamber. After the incubation, the PDMS chip and clamp were removed after a 5-min wash with 1× NEBuffer 3.1. The slide was subsequently washed with water and then air-dried.
For the second in situ ligation of barcodes (barcodes B), the second PDMS chip was placed over the slide. A bright-field image was captured using a 10× objective for alignment purposes. The PDMS device and tissue slide were securely clamped together using a clamp. The annealing of barcodes B (B1–100, 25 μM) and the preparation of the ligation master mix were performed identically to as described for barcodes A. The device was then incubated at 37 °C for 30 min in a humidified chamber. After incubation, the PDMS chip and clamp were removed after a 5-min wash with 1× DPBS containing SUPERase In RNase inhibitor. The slide was then rinsed with water and air-dried. A bright-field image was taken afterwards to assist with further alignment.
The ligation of barcodes A and B for the 220 × 220 barcode chips was performed according to the manufacturer’s protocol (AXO-0457).
After the barcoding process, the ROI on the tissue was treated with a reverse-crosslinking mixture containing 1 mM EDTA, 0.4 mg ml−1 proteinase K, 50 mM Tris-HCl (pH 8.0), 1% SDS and 200 mM NaCl. This mixture was incubated at 58 °C for 2 h in a humidified chamber. The lysate was then transferred to a 1.5 ml tube and incubated at 65 °C overnight.
To separate DNA and cDNA, the lysate was purified using the Zymo DNA Clean & Concentrator-5 kit and then eluted into 100 μl of RNase-free water. The Dynabeads MyOne Streptavidin C1 beads (40 μl) were washed three times using 1× B&W buffer containing 0.05% Tween-20. The beads were then resuspended in 100 μl of 2× B&W buffer with 2.5 μl of SUPERase In RNase inhibitor. This bead suspension was then mixed with the lysate and allowed to bind at room temperature for 1 h with agitation. Finally, a magnet was used to separate the beads from the supernatant in the lysate.
The supernatant was processed for ATAC–seq library construction and purified using the Zymo DNA Clean & Concentrator-5 kit, then eluted into 20 μl of RNase-free water. For the PCR, 2.5 μl of 25 μM indexed i7 primer, 25 μl of 2× NEBNext Master Mix and 2.5 μl of 25 μM P5 PCR primer were added and thoroughly mixed. The initial PCR program was set as follows: 72 °C for 5 min; 98 °C for 30 s; followed by 5 cycles of 98 °C for 10 s, 63 °C for 30 s and 72 °C for 1 min. To determine the need for additional cycles, 5 μl of the pre-amplified mixture was combined with a quantitative PCR (qPCR) solution consisting of 0.5 μl of 25 μM new P5 PCR primer, 5 μl 2× NEBNext Master Mix, 3.76 μl nuclease-free water, 0.24 μl 25× SYBR Green and 0.5 μl of 25 μM indexed i7 primer. The qPCR reaction was then performed using the following program: 98 °C for 30 s; followed by 20 cycles of 98 °C for 10 s, 63 °C for 30 s and 72 °C for 1 min. On the basis of the qPCR results, the remaining 45 μl of pre-amplified DNA was further amplified for additional cycles until it reached one-third of the saturation signal. The final PCR product was then purified using 45 μl of 1× Ampure XP beads and eluted in 20 μl of nuclease-free water.
The beads were used to construct cDNA libraries. Initially, they were washed twice with 400 μl of 1× B&W buffer containing 0.05% Tween-20 and then once with 10 mM Tris (pH 8.0) with 0.1% Tween-20. The Streptavidin beads, with cDNA molecules bound, were resuspended in a TSO solution composed of 44 μl of 5× Maxima RT buffer, 22 μl of 10 mM dNTPs, 44 μl of 20% Ficoll PM-400 solution, 5.5 μl of RNase inhibitor, 5.5 μl of 100 μM template switch primer, 88 μl of RNase-free water, and 11 μl of Maxima H Minus reverse transcriptase. The beads were incubated at room temperature for 30 min and then at 42 °C for 90 min with gentle shaking. Next, the beads were washed once with 10 mM Tris containing 0.1% Tween-20 and once with water, then resuspended in a PCR solution consisting of 91.9 μl of RNase-free water, 8.8 μl of 10 μM each of primers 1 and 2, 0.5 μl of 1 μM primer 2-citeseq and 110 μl of 2× Kapa HiFi HotStart Master Mix. PCR thermocycling was conducted with an initial denaturation at 95 °C for 3 min, followed by 5 cycles of 98 °C for 20 s, 65 °C for 45 s and 72 °C for 3 min. After the initial five cycles, the Dynabeads MyOne Streptavidin C1 beads were removed from the PCR solution and 25× SYBR Green was added at a 1× concentration. The samples were then subjected to further qPCR with the programme set at 95 °C for 3 min, cycled at 98 °C for 20 s, 65 °C for 20 s, and 72 °C for 3 min for 15 cycles, with a final extension at 72 °C for 5 min. The reaction was terminated once the qPCR signal began to plateau.
To differentiate and purify cDNAs derived from mRNA and ADT, we used 0.6× SPRI beads according to the standard protocol. We mixed 120 µl of SPRI beads with 200 µl of the PCR product solution and allowed the mixture to incubate for 5 min. The supernatant, which contained the ADT-derived cDNAs, was then transferred to a 1.5 ml Eppendorf tube. The beads remaining in the original container were washed with 85% ethanol for 30 s and eluted in RNase-free water over 5 min. The mRNA-derived cDNAs were subsequently quantified using Qubit and BioAnalyzer. For further purification of the supernatant, we added 1.4× SPRI beads and incubated for 10 min. These beads were cleaned once with 80% ethanol and resuspended in 50 µl of water. An additional purification step was conducted by adding 100 µl of 2× SPRI beads and incubating for another 10 min. After two washes with 80% ethanol, the ADT-derived cDNAs were finally eluted with 50 µl of RNase-free water.
The sequencing libraries for the two types of cDNA products were constructed separately. For the mRNA library, we used the Nextera XT Library Prep Kit. We began by diluting 1 ng of purified cDNA in RNase-free water to a final volume of 5 μl. We then added 10 μl of Tagment DNA buffer and 5 μl of Amplicon Tagment Mix, and incubated the mixture at 55 °C for 5 min. Next, 5 μl of NT buffer was added and the mixture was incubated at room temperature for another 5 min. The PCR master solution, consisting of 15 μl PCR master mix, 1 μl of 10 μM P5 primer, 1 μl of 10 μM indexed P7 primer and 8 μl RNase-free water, was then added. The PCR was then conducted with the following programme: an initial denaturation at 95 °C for 30 s; followed by 12 cycles of 95 °C for 10 s, 55 °C for 30 s, 72 °C for 30 s and a final extension at 72 °C for 5 min. The PCR product was purified using 0.7× Ampure XP beads to obtain the final library.
For the ADT cDNAs, the library was constructed using PCR. In a PCR tube, 45 µl of the ADT cDNA solution was mixed with 10 µM 2.5 µl of P5 index, 2.5 µl of a 10 µM customized i7 index and 50 µl of 2× KAPA HiFi PCR Master Mix. The PCR conditions were set as follows: an initial denaturation at 95 °C for 3 min, followed by six cycles of denaturation at 95 °C for 20 s; annealing at 60 °C for 30 s; extension at 72 °C for 20 s; and a final extension at 72 °C for 5 min. The PCR product was then purified using 1.6× SPRI beads.
The size distribution and concentration of the library were analysed using the Agilent Bioanalyzer High Sensitivity Chip before sequencing. Next-generation sequencing (NGS) was then performed using the Illumina NovaSeq 6000 sequencer in paired-end 150 bp mode.
Spatial CTRP-seq
Tissue slide preparation and antibody incubation were performed according to the same protocol as used in the spatial ARP-seq method. After these initial steps, the tissue was washed twice with wash buffer (0.5 mM spermidine, 150 mM NaCl, one tablet of protease inhibitor cocktail, 20 mM HEPES, pH 7.5) and was subsequently incubated in NP40-Digitonin wash buffer (0.01% NP40, 0.01% Digitonin, 0.5 mM spermidine, 150 mM NaCl, one tablet of protease inhibitor cocktail, 20 mM HEPES pH 7.5) for 5 min, followed by three washes with 0.5× DPBS-RI. The section was then washed with wash buffer (0.5 mM spermidine, 150 mM NaCl, one tablet of protease inhibitor cocktail, 20 mM HEPES pH 7.5). The primary antibody was diluted 1:50 in antibody buffer (2 mM EDTA, 0.001% BSA, in NP40-Digitonin wash buffer) and applied to the tissue, which was then incubated overnight at 4 °C. The secondary antibody (guinea pig anti-rabbit IgG), diluted 1:50 in NP40-Digitonin wash buffer, was then added and incubated for 30 min at room temperature. The tissue was washed with wash buffer for 5 min. Next, a 1:100 dilution of the pA-Tn5 adapter complex in 300-wash buffer (20 mM HEPES pH 7.5, 300 mM NaCl, one tablet of protease inhibitor cocktail, 0.5 mM spermidine) was added and incubated at room temperature for 1 h, followed by a 5 min wash with 300-wash buffer. Tagmentation buffer (10 mM MgCl2 in 300-wash buffer) was then added and incubated at 37 °C for 1 h. Finally, to stop the tagmentation, 40 mM EDTA with 0.05 U μl−1 RNase inhibitor was added and incubated at room temperature for 5 min. The tissue was washed three times with 0.5× DPBS-RI for 5 min for further processing.
For the reverse transcription, two ligations and beads separation, the procedures were performed as described for those established in the spatial ARP-seq protocol. For constructing the CUT&Tag library, the supernatant was purified using Zymo DNA Clean & Concentrator-5 kit and eluted into 20 µl of RNase-free water. The PCR mix, consisting of 2 µl of 10 µM each of P5 PCR primer and indexed i7 primer along with 25 µl of NEBNext master mix, was added and thoroughly mixed. The PCR programme included an initial incubation at 58 °C for 5 min, 72 °C for 5 min and 98 °C for 30 s, followed by 12 cycles of 98 °C for 10 s and 60 °C for 10 s, with a final extension at 72 °C for 1 min. The PCR product was then purified using 1.3× Ampure XP beads according to the standard protocol and eluted in 20 µl of nuclease-free water. The construction of the cDNA libraries for mRNA and ADT was carried out according to the spatial ARP-seq protocol.
Before sequencing, the size distribution and concentration of the library were assessed using the Agilent Bioanalyzer High Sensitivity Chip. Subsequently, NGS was performed using the Illumina NovaSeq 6000 sequencer in paired-end 150 bp mode.
Quality control metrics for spatial ARP-seq and spatial CTRP-seq
For spatial ARP-seq, across all samples from mouse development, the LPC model and human brain V1 development, we observed a median of 11,635 unique fragments per pixel for ATAC modality. Among these, 19% were enriched at transcription start sites (TSS) and 18% were located in peaks (Supplementary Figs. 15a and 16a). For the RNA portion, a total of 23,824 genes were detected with an average of 1,230 genes and 2,392 unique molecular identifiers (UMIs) per pixel (Supplementary Fig. 16a). For the proteins, the average protein count per pixel was 59, and the protein UMI account per pixel was 394.
Spatial CTRP-seq for H3K27me3 was performed on 5 and 21 d.p.l. mouse brains. We obtained a median of 9,102 unique fragments per pixel, of which 10% of fragments overlapped with TSS regions, and 12% were located in peaks (Supplementary Figs. 15a and 16a). For the RNA data, a total of 23,397 genes was detected with an average of 1,318 genes per pixel and 2,392 UMIs per pixel (Supplementary Fig. 16a). For the proteins, the average protein count per pixel was 88, and the protein UMI counts per pixel was 793.
The insert size distributions of spatial ARP-seq and CTRP-seq (H3K27me3) fragments were consistent with the captured nucleosomal fragments in all tissues (Supplementary Fig. 15b). The correlation analysis between replicates showed high reproducibility (r = 0.99 for ATAC, r = 0.98 for RNA and r = 0.99 for protein in spatial ARP-seq; r = 0.98 for CUT&Tag (H3K27me3), r = 0.95 for RNA and r = 0.99 for protein in spatial CTRP-seq, where r is the Pearson correlation coefficient; Supplementary Fig. 16b,c).
Data preprocessing
For ATAC and RNA, as well as CUT&Tag and RNA data, linker 1 and linker 2 were used to filter read 2, and the sequences were converted to Cell Ranger ARC format (v.2.0.2, 10x Genomics). The genome sequences were in the newly formed read 1, barcodes A and barcodes B were included in newly formed read 2. The human reference (GRCh38) or mouse reference (GRCm38) genome was used to align the fastq files.
For protein data from ADTs, the raw FASTQ file of read 2 was reformatted in the same manner as for RNA. Using the default configurations of CITE-seq-Count v.1.4.2, we quantified the ADT UMI counts associated with each antibody at various spatial locations.
Data clustering and visualization
Initially, we determined the locations of pixels on the tissue sections from brightfield images using MATLAB 2020b. Subsequently, we generated spatial data files by using codes available at GitHub (https://github.com/di-0579/spatial_tri-omics).
Seurat (v.4.1)92 was loaded in R (v.4.1) to construct Seurat objects for each sample using RNA matrices. Attribute information was then added to each object, and the objects were merged using the merge function. Pixel normalization, logarithmic transformation of counts and gene scaling were performed, followed by PCA, retaining the top 50 principal components. To reduce sample heterogeneity, the RunHarmony function was used for integration, and RunUMAP and FindClusters were applied to generate RNA clustering information. Finally, metadata were extracted to enable the spatial visualization of RNAcluster for each sample.
We used Library Signac (v.1.8)93 in R (v.4.1) to read the ATAC data matrices. To integrate multiple datasets, functions from the GenomicRanges package were used to establish a common peak set. Fragment objects were then created for each sample, peaks were quantified across datasets and these quantified matrices were used to create a Seurat object for each dataset, with the Fragment object stored within each respective assay. The datasets were merged using the merge function. To minimize sample heterogeneity, we used FindIntegrationAnchors to identify integration anchors, using LSI embeddings for integration. Subsequently, we performed ATAC clustering and extracted metadata for spatial visualization of ATACcluster.
We used the SpatialGlue94 package to integrate and cluster ATAC and RNA data, obtaining spatial domain information. The ArchR95 package was used to load both ATAC and RNA data, incorporating the spatial domain information from SpatialGlue. This enabled us to identify differential RNA expression and the GAS across various regions. We used the addImputeWeights() function in ArchR to compute smoothing weights based on low-dimensional embeddings derived from the GeneExpressionMatrix within the ArchR project, which was selected for better spatial domain resolution. We then extracted RNA expression matrices, ATAC GAS expression matrices and CUT&Tag CSS expression matrices from ArchR objects for spatial gene visualization.
Spatiotemporal data analysis
To better identify the spatiotemporal patterns of gene regulation using our developed spatial multi-omics technology, we implemented a computational framework that consists of the following three main steps:
-
(1)
Spatiotemporal regression model fitting
For each gene in RNA with expression in more than 2% of pixels across timepoints, we fit a negative binomial generalized additive model (NB-GAM) to account for the effects of library size, time and spatial location on gene expression levels and smooth terms were used to represent the effects of time and spatial location on gene expression, as well as their interaction96. Specifically, for the RNA counts Ygc for a given gene g and pixel c with gene-specific means μgc and dispersion parameters ϕg, we have
$${Y}_{{gc}} \sim {NB}({\mu }_{{gc}},{\phi }_{g}),$$
$$\log ({\mu }_{{gc}})={\beta }_{g0}+{\beta }_{g1}{N}_{c}+{f}_{g1}(tc)+{f}_{g2}(sc)+{f}_{g3}(tc,sc),$$
where ~ indicates ‘distributed as’, N indicates the sequencing depth; tc indicates the timepoint information (P0, 2, 5, 7, 10, 21) of pixel c; sc indicates the spatial location of pixel c. In the analysis of the CC, the spatial location was defined by manually dividing the CC region into ten bins after manually delineating the region. In the analysis of cortical layers, the spatial location was defined using the SpatialGlue results. For a given g, fg1() is a smooth function capturing temporal variations, fg2() is a smooth function capturing spatial variations and fg3() is a tensor product smooth function capturing the interaction between time and spatial location. We used the cubic regression spline as the marginal basis of the smooth function. The model is fitted using gam(), with ti() as the smooth function and family as nb() from R package mgcv.
For ATAC data, we first derived the gene score matrix by aggregating the peak accessibility for peaks falling within a window ±50 kb of the TSS around each gene and expressed in more than 1% pixels across five timepoints, using the function getDORCScores() from R package FigR97. We smoothed the sparse gene score matrix per pixel per gene using its four nearest neighbours based on the first 50 principal components from RNA data, using the function smoothScoresNN(). We then fitted the GAM on the log-scaled smoothed gene score matrix with the Gaussian family.
This model enables us to identify genes with significant associations with either time, spatial location or their interaction. We determined the significance using Wald tests for the model coefficients and used the Benjamini–Hochberg procedure to adjust the P values. Genes with an adjusted P value of less than 0.01 in any of the spatial, time or time and spatial interaction terms were selected for further analysis.
-
(2)
Gene spatiotemporal clustering
We performed a two-step procedure to group patterns of genes that show significant changes across time and space in either RNA or ATAC. First, we derived pixel-level estimates μgc from our statistical model and aggregated them by their respective timepoints and spatial location regions, calculating RNA and ATAC separately. These aggregated profiles were then scaled. Next, we concatenated the RNA and ATAC aggregated profiles and then performed hierarchical clustering using one minus Pearson correlation coefficient as the distance metric and Ward’s minimum variance as a linkage method98, with the number of clusters set as 20 in the analysis of the CC and set as 30 in the analysis of cortical layers. This procedure aims to jointly capture both ATAC and RNA patterns. The number of clusters was intentionally set higher to overestimate the number of clusters, enabling us to capture as many joint patterns of RNA and ATAC as possible.
As we observed that some patterns in RNA and ATAC are similar, we subsequently combined the clusters on the basis of their similarity in RNA and ATAC profiles, respectively. First, we aggregated the RNA and ATAC profiles by their gene clustering labels from the previous step. We then performed hierarchical clustering with a smaller number of clusters: 6 for RNA and 4 for ATAC in the CC analysis (integrated S1 and S2 data), and 15 for RNA and 10 for ATAC for S2 data, and 13 for RNA and 7 for ATAC for S1 data in the cortical layers analysis. The number of clusters was determined by assessing whether the main patterns belonged to separate clusters and by evaluating the interpretability of the results. For the previous overclustered results in the previous step, we assigned labels to represent their RNA and ATAC patterns. If clusters shared the same RNA and ATAC pattern labels, we combined these clusters. Each cluster is then labelled based on its RNA and ATAC pattern, for example, R1-A1 for a cluster where the RNA pattern belongs to R1 and the ATAC pattern belongs to A1. This step summarizes the RNA and ATAC patterns and improves the interpretability of each cluster. By characterizing genes with similar RNA and ATAC patterns separately, we can also jointly examine clusters with different RNA and ATAC patterns across space and time. For example, clusters labelled R1-A1 and R1-A3 share the same RNA pattern but have different ATAC patterns, while clusters labelled R1-A1 and R2-A1 share the same ATAC pattern but have different RNA patterns.
-
(3)
GO enrichment analysis of gene spatiotemporal clustering
We performed GO enrichment analysis to further characterize the identified clusters, focusing on the orthogonal ontology of biological processes. The analysis was performed using the function enrichGO() from the R package clusterProfiler99 with minGSSize set to 20 and maxGSSize set to 200.
CODEX data clustering
Whole-cell segmentation was performed with Cellpose100 using the Cytoplasm model cyto. After obtaining the cell-level protein intensity matrix, we used Seurat v.4.1 to perform the subsequent analysis. For each sample, the data are first normalized using arcsinh transformation of the intensity matrix that scaled with a cofactor of 150 (ref. 101) and then scaled using ScaleData() and RunPCA() with all informative features. We removed the sample effect using the FindIntegrationAnchors() function, with reciprocal PCA (rpca) set as the dimensional reduction method to identify anchors and the first 20 principal components used as the dimensions. After integration, we performed clustering using the FindClusters() function on the neighbour graphs built based on the first 20 PCs, with the resolution set to 0.4. Each cluster was then manually annotated on the basis of their highly expressed markers.
Cell type and niche deconvolution
We performed cell2location (v.0.1.3)102 to deconvolute the cell types of our spatial transcriptomics data using public references in Python v.3.9. We combined all samples in mouse development or in mouse LPC model together. We set the model ‘cell2location.models. Cell2location()’ with the expected average cell abundance as 5 and trained with full data with maximum epochs set as 30,000. The estimated cell abundance is based on the posterior mean using export_posterior().
To map cell types and spatial niches from MERFISH reference data17 to our human visual cortex spatial transcriptomics data, we implemented Seurat’s (v.4.3) label-transfer pipeline. For each developmental stage, we first identified transfer anchors using canonical correlation analysis on the top 30 principal components. We then transferred both cell type and spatial niche annotations from reference to query data using the TransferData() function. To ensure annotation reliability, pixels with maximum prediction scores below the 25th percentile threshold were excluded and isolated or sparsely distributed pixels were eliminated.
Selection of lesion-like compartments
To characterize microglial responses to demyelination in both primary lesion and distal regions, we implemented a stepwise strategy to identify spatially resolved lesion-like compartments. We first used cell2location-predicted cell type abundances, combined with spatial expression patterns of canonical microglial markers Cx3cr1, P2ry12 and Tmem119 across RNA and ATAC modalities, to define high-confidence microglia pixels. To remove random pixels not located in primary and distal regions, we applied a k-nearest neighbour (k-NN) filter, retaining core pixels of which the average distance to their five nearest neighbours was within the 85–90th percentile. Core regions were then expanded by including surrounding pixels within a 10–15-neighbour k-NN window (around 80–100 µm) to capture broader lesion-like compartment.
Multimodal clustering and differential analysis of lesion-like compartments
To further dissect the heterogeneity of microglia within primary and distal lesion-like compartments, we performed joint RNA and ATAC analysis using ArchR95 for all the LPC samples. Dimensionality reduction was performed using iterative latent semantic indexing (LSI) separately on ATAC and RNA matrices. Harmony batch correction was applied to each modality using the sample identity as the grouping variable. We then integrated both modalities using the addCombinedDims() function to obtain a joint reduced embedding (LSI_Combined). UMAP visualization and clustering were performed on this combined space using a resolution of 1.0. The resulting multimodal clusters were visualized using cluster-specific UMAP plots and assigned labels LC1–LC5 for downstream interpretation. To identify marker genes for each subset (LC1–LC5), we applied getMarkerFeatures() to both the gene activity (GeneScoreMatrix) and RNA expression (GeneExpressionMatrix) using the Wilcoxon test. GO enrichment analysis was then performed on the cluster-specific marker genes to uncover functional distinctions across clusters. To compare protein-level differences across subsets, we applied the ArchR derived LC1–LC5 cluster labels to the ADT data. Expression levels of selected markers were visualized using the DotPlot() function from Seurat, based on normalized ADT assay values.
Gene module scoring of lesion-like compartments
To assess microglia in primary lesion and distal compartments, we constructed gene modules based on two previously published studies75,76 (Supplementary Table 12). We computed module scores independently for RNA and ATAC data using the addModuleScore() function in ArchR95, applied to the GeneExpressionMatrix and GeneScoreMatrix, respectively. To enhance spatial coherence of gene module scores, we applied nearest-neighbour smoothing using the addImputeWeights() function, with weights computed from LSI embeddings derived from the corresponding GeneExpressionMatrix and GeneScoreMatrix. We then extracted RNA and ATAC module score matrices from ArchR objects for spatial visualization.
Subclustering of microglia and identification of differentially expressed genes
We selected microglial-specific pixels based on cell2location predictions, retaining only those in which microglia showed the highest abundance among all cell types. Using these microglial-specific pixels, we performed multimodal subclustering by applying low-dimensional embedding (LSI and Harmony) and clustering on the combined RNA–ATAC data. This approach identified three distinct subclusters (MC1, MC2 and MC3). Differential analysis between subclusters was performed for both RNA expression (GeneExpressionMatrix) and chromatin accessibility (GeneScoreMatrix) using the getMarkerFeatures() function in ArchR95, with thresholds of FDR ≤ 0.1 and log2[FC] ≥ 0.1. Volcano plots were generated to visualize the log2 fold changes and statistical significance (−log10[FDR]) for each modality. Genes exhibiting significant differences in both RNA and ATAC were highlighted to support joint interpretation of transcriptional and epigenetic states across the microglia subtypes.
Cell–cell communication analysis using CellChat
To explore intercellular communication among microglial subclusters (MC1–MC3) and other surrounding cell types, we applied CellChat103 to spatial RNA data from 5 d.p.l. and 10 d.p.l. samples. The analysis focused on SCT-normalized RNA expression values from lesion-like compartments. Each pixel was assigned a specific cell type identity based on cell2location predictions, with microglia pixels labelled as MC1–MC3 according to their respective subclusters. We constructed a spatial CellChat object (datatype = “spatial”) with spatial filtering enabled (distance.use = TRUE, interaction.range = 250, scale.distance = 0.5, contact.dependent = TRUE). Communication probabilities were computed using the truncated mean method (10% trimming), and significant interactions were filtered using min.cells = 2 and adjusted P < 0.05. All significant interactions (P < 0.05) across cell types were visualized using netVisual_bubble(), highlighting signalling associated with MC1–MC3 microglia.
Ligand–receptor analysis using NICHES
We applied the NICHES104 framework to lesion-like compartments at 10 d.p.l. using the ALRA-imputed expression matrix (alra_SCT) as input to denoise sparse spatial expression. NICHES was run with the CellToCellSpatial, CellToNeighborhood and NeighborhoodToCell modes enabled, using the FANTOM5 ligand–receptor database, cell annotation and a neighbourhood size of k = 9.
To analyse cell–cell signalling patterns based on the NeighborhoodToCell_SCT_9_alra assay (produced by NICHES), we performed unsupervised clustering including variable feature selection, scaling, PCA (100 components), neighbourhood graph construction and Louvain clustering (resolution = 0.5). A UMAP embedding was generated using the top 30 principal components. Cluster-specific ligand–receptor signatures were identified using FindAllMarkers with a minimum expression threshold (min.pct = 0.5) and log-transformed fold change threshold (logfc.threshold = 0.25). The top 10 and top 30 marker genes per cluster were ranked by average log2-transformed fold change and exported as CSV files. SpatialFeaturePlots were generated for all top marker genes across spatial coordinates to visualize the spatial distribution of key ligand–receptor signals.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
