Spatial and temporal heterogeneity in the lineage progression of fine oligodendrocyte subtypes

Background Oligodendrocytes are glial cells that support and insulate axons in the central nervous system through the production of myelin. Oligodendrocytes arise throughout embryonic and early postnatal development from oligodendrocyte precursor cells (OPCs), and recent work demonstrated that they are a transcriptional heterogeneous cell population, but the regional and functional implications of this heterogeneity are less clear. Here, we apply in situ sequencing (ISS) to simultaneously probe the expression of 124 marker genes of distinct oligodendrocyte populations, providing comprehensive maps of the corpus callosum, cingulate, motor, and somatosensory cortex in the brain, as well as gray matter (GM) and white matter (WM) regions in the spinal cord, at postnatal (P10), juvenile (P20), and young adult (P60) stages. We systematically compare the abundances of these populations and investigate the neighboring preference of distinct oligodendrocyte populations. Results We observed that oligodendrocyte lineage progression is more advanced in the juvenile spinal cord compared to the brain, corroborating with previous studies. We found myelination still ongoing in the adult corpus callosum while it was more advanced in the cortex. Interestingly, we also observed a lateral-to-medial gradient of oligodendrocyte lineage progression in the juvenile cortex, which could be linked to arealization, as well as a deep-to-superficial gradient with mature oligodendrocytes preferentially accumulating in the deeper layers of the cortex. The ISS experiments also exposed differences in abundances and population dynamics over time between GM and WM regions in the brain and spinal cord, indicating regional differences within GM and WM, and we found that neighboring preferences of some oligodendroglia populations are altered from the juvenile to the adult CNS. Conclusions Overall, our ISS experiments reveal spatial heterogeneity of oligodendrocyte lineage progression in the brain and spinal cord and uncover differences in the timing of oligodendrocyte differentiation and myelination, which could be relevant to further investigate functional heterogeneity of oligodendroglia, especially in the context of injury or disease. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-022-01325-z.


Background
Classically, oligodendrocytes (OLs) have been considered a homogeneous population in the central nervous system (CNS). Recent advances in single-cell RNA sequencing (scRNA-seq) have allowed systematic studies of the transcriptome in individual cells, which has changed the paradigm of cellular classification [1]. At least twelve distinct OL populations can be distinguished in the CNS [2]. If these transcriptionally distinct populations are indeed different cell types or rather states remains a recurring question in the field [3], and their spatial distribution could give further insights into the heterogeneity of OLs.
The OL lineage presents itself as a continuum of oligodendrocyte precursor cells (OPCs), committed OPCs (COPs), newly formed OLs (NFOLs), myelin-forming OLs (MFOLs), and mature OLs (MOLs) [2]. In our previous study, focusing on a limited number of markers by applying RNASCOPE in situ hybridization, we have shown that OPC, MOL2, and MOL5/6 take spatial preference in the CNS, and the abundance changes with age [4]. More precisely in the neocortex and corpus callosum, the OPC population decreases with age while the MOL5/6 population increases in both regions [4]. Also, in the spinal cord, the MOL5/6 population increases with age, and it shows a spatial preference for gray matter, while the MOL2 population shows a spatial preference for white matter [4].
Most recently, novel spatial transcriptomics methods with high multiplexy have been developed that can in principle map the location of each gene and each cell [5,6]. Multiple ex situ and in situ spatial transcriptomics technologies are now available. In situ sequencing (ISS) is a targeted, image-based approach that allows the precise mapping of genes in their natural context and offers cellular resolution [7]. Here, we applied probabilistic cell typing by in situ sequencing (pciSeq) [8] to target OLs for the first time, leveraging scRNA-seq data to assign detected transcripts to segmented cells and, subsequently, cells to cell types. We used the full advantages of ISS by probing 124 marker genes that characterize all twelve previously described populations of the OL lineage and the vascular and leptomeningeal cells (VLMCs) [2]. We investigated multiple brain and spinal cord tissue sections and found region-and age-related differences in the oligodendrocyte lineage progression, uncovering the differences in the timing of oligodendrocyte differentiation and myelination in these areas.

In situ sequencing allows the detection of major OL populations
We performed in situ sequencing (ISS) in postnatal (P10), juvenile (P20), and young adult (P60) mouse brain sections as well as thoracic spinal cord sections, and we quantified the OL lineage populations in their white matter (WM) vs. gray matter (GM), respectively. In the first batch, we mapped 6 P20-and 6 P60-hemisected coronal brains including the neocortex (CTX) and corpus callosum (CC) from 8 mice in total (4 biological replicates each). In the second batch, we mapped 4 P10-hemisected coronal sections from an additional 4 mice. We also mapped 4 P20 and 4 P60 spinal cord sections. The gene maps of CC (WM) in a P60 brain section are shown in Fig. 1A, with the colored symbols representing unique transcripts. The 124 genes were selected from single-cell RNA sequencing data from [2] based on the median expression values to target major markers for the oligodendrocyte populations (Fig. 1C). The color of the assigned gene symbols represents the cell population with the highest expression of the respective gene; however, we used a probabilistic gene to cell type assignment strategy. Applying pciSeq [8], we assigned detected transcripts to segmented cells and, subsequently, based on their genetic signature cells to cell populations. Thereby, each cell is represented as a pie chart, with the angle of each slice being proportional to the cell type probability. The highest probability in the pie chart defines the final cell population. The mean probability, i.e., the mean certainty, for the OL assignments was 77%, with similar probabilities for the different cell populations (Additional file 1: Fig. S1A). Only cells that express Plp1 as a PAN marker for OLs were considered for downstream analysis. Example cells that are cells with the median count of genes per population are shown in Fig. 1B. On average, 13.4 ± 4.1 transcripts and 9.4 ± 0.4 distinctive genes per cell were measured in the OL populations in P10 brains, 14.6 ± 2.1 transcripts and 10.0 ± 0.1 distinctive genes per cell were measured in P20 brains, and 17.4 ± 2.7 transcripts and 9.8 ± 0.1 distinctive genes per cell were measured in P60 brains (Additional file 1: Fig.  S1B, C). Heatmaps of the log2 expression per gene and the most differentially expressed genes per cell population are shown in Additional file 2: Fig. S2, emphasizing that many OL-associated genes are shared between OL populations but that the composition of the genes per OL population is unique, allowing us to map these populations with pciSeq.
Keywords: Oligodendrocytes, Lineage progression, Spatial transcriptomics, In situ sequencing, Cortex, Corpus callosum, Spinal cord Fig. 1 Oligodendrocyte types in the neocortex (CTX), grey matter (GM), corpus callosum (CC), and white matter (WM). a Maps of 124 genes in coronal mouse brain sections (P60) targeted by in situ sequencing, here highlighting the CC. Three different zoom levels are shown with white dashed lines denoting the outlines for CTX and CC. The red boxes mark the positions of the zoom-ins. The scale bars are 500 μm (left), 100 μm (middle), and 15 μm (right). b Individual example cells of P60 CC and their gene assignments to cells are shown. The cells are classified based on single-cell RNA sequencing (scRNA-seq) data and represented as pie charts where the angle of each slice is proportional to the likelihood of the cell being of a certain type and the size of the pie chart is indicative of the number of transcripts. Examples shown are median cells for each cell population, i.e., cells with the median count of the number of transcripts being assigned within each type. The most differentially expressed genes of the example cells are shown in the insets, and the probability for each putative cell type from pciSeq (probabilistic cell typing by in situ sequencing) is listed below the example cells. The scale bar is 5 μm. c The 124 genes are grouped based on the expression values in the scRNA-seq data. The color of the assigned gene symbols represents the cell population with the highest expression of the respective gene. The individual genes in the heatmaps are ranked by their expression values Representative cell maps of OL populations are shown in Fig. 2A (P60) and Additional file 3: Fig. S3. The size of each pie chart in Fig. 2A is indicative of the number of assigned transcripts. We also generated the mean pie charts counting how often a certain OL population was called with other OL populations (Fig. 2B). The mean pie charts per OL population over the brain tissues (n = 6 tissue sections) reflect the progression of OL lineage very well, with OPCs often being shared with COPs, NFOLs with MFOLs, and MFOLs with MOLs. The relatedness of the six MOL populations can also be seen in the mean confusion matrix of the cell populations during the probabilistic assignment (Additional file 4: Fig. S4).

Myelination advances earlier in GM (CTX) than in WM (CC) in the brain
In total, we captured 2693 OLs in P10 CTX, 2127 OLs in P10 CC, 4209 OLs in P20 CTX, 3879 OLs in P20 CC, 9129 OLs in P60 CTX, and 7139 OLs in P60 CC (n = 4 P10 tissue sections, n = 6 P20 tissue sections, n = 6 P60 tissue sections, Fig. 2C). For P10, we captured a tissue area of 3.2 ± 0.3 mm 2 of CTX and a tissue area of 0.8 ± 0.1 mm 2 of CC, yielding an OL density of 207.2 ± 16.3 cells/mm 2 in CTX and an increased OL density of 648.4 ± 134.6 cells/mm 2 in CC (p = 0.002, n = 4 tissue sections). For P20, we captured a tissue area of 3.6 ± 0.3 mm 2 of CTX and a tissue area of 0.9 ± 0.1 mm 2 of CC, yielding an OL density of 193.5 ± 31.8 cells/mm 2 in CTX and an increased OL density of 750.4 ± 82.8 cells/mm 2 in CC (p = 0.0003, n = 6 tissue sections). For P60, we captured a tissue area of 8.9 ± 0.6 mm 2 of CTX and a tissue area of 2.0 ± 0.1 mm 2 of CC, with an OL density of 170.1 ± 48.2 cells/mm 2 in CTX and an increased OL density of 608.6 ± 190.4 cells/mm 2 in CC (p = 0.05, n = 6 tissue sections), as expected for gray and white matter regions.
In order to compare the abundances of the OL populations between GM (CTX) and WM (CC) at P10, P20, and P60, we calculated both the cell occurrences per Plp1+ cells (%, Fig. 2D-F) and the cell occurrences per tissue area (cells/mm 2 , Additional file 5: Fig. S5A-C). The bar plots of the OL abundances show that at P10 and P20 stages, when myelination is actively ongoing, OPC, COP, NFOL1, and MFOL2 are most abundant in CTX and CC ( Fig. 2D, E), while at P60, MOL5 are most abundant in CTX and CC and MOL1 and MOL4 populations are a minority (Fig. 2F). Calculating the rate of change of cells/ mm 2 for each cell population further highlights that the spatial density in P10 CTX is in general higher than in P20 CTX, except for MOL5 and VLMC and some populations that remain stagnant (NFOL2, MOL2, MOL3) (Fig. 3A, left). In P10 CC, the cell densities of NFOL2 and MFOL1 are decreased compared to P20 CC as well as the cell densities of MOL5 and VLMC (Fig. 3A, right). Interestingly, comparing P20 CTX vs. P60 CTX shows that the cell densities of MOL5 and MOL6 strongly increase, while the other densities largely remain stagnant (Fig. 3B, left). Similarly, comparing P20 CC vs. P60 CC shows that the cell densities of MOL5 and MOL6 strongly increase, together with the cell density of MOL2, while the other densities largely remain stagnant (Fig. 3B, right). This could suggest that at the early stages from P10 to P20 (especially in CTX), the decrease in cell density of OPC, COP, NFOL1, and MFOL2 could reflect their differentiation into MOLs, whereas at the later stages, the cell densities for early-stage OLs remain mostly unchanged, which to some extent could reflect their ability to proliferate and to replace those cells that have differentiated.
To further understand the relative changes of OL populations in respect to the overall number of OLs and to compare the abundances across ages and regions, we normalized the number of cells for each OL population by the number of Plp1+ cells. We then calculated how much the OL population in the more prevalent region (P10, P20, P60 CTX and CC) is elevated compared to the non-prevalent region. In P10, we found that MFOL1 (59% ± 7% increased in GM versus WM) and MFOL2 (36% ± 12%) are more abundant in GM (CTX) compared to WM (CC), while COP (39% ± 7%) and NFOL1 (22% ± 4%) are increased in WM compared to GM (p = 0.009, p = 0.01, p = 0.03, p = 0.04, respectively; n = 4 tissue sections) (Fig. 3C, left). In P20, we found that NFOL2 (25 ± 8%), MFOL1 (25 ± 7%), and MFOL2 (22 ± 8%) are more abundant in GM (CTX) compared to WM (CC), while there is an enrichment of COPs (31 ± 8%) and VLMC (56 ± 4%) in CC compared to CTX (p = 0.03, p = 0.004, p Oligodendrocyte map and abundances in P10, P20, and P60 CTX (GM) and CC (WM). a Spatial map of Plp1+ cell types in a P60 coronal mouse brain section including neocortex (CTX) and corpus callosum (CC). The scale bar is 500 μm. b The mean pie charts for P60 (CTX and CC joined) generated by counting how often a certain OL population was called with other OL populations. Note: Considering the top pies within the pie charts (caption below), the charts highlight the oligodendrocyte lineage from OPC to MOL6. c The number of captured oligodendrocytes as well as the relative distribution of the oligodendrocyte types per age condition and brain region (n = 4 P10 tissue sections, n = 6 P20 tissue sections, n = 6 P60 tissue sections). = 0.02, p = 0.05, p = 0.02, respectively; n = 6 tissue sections) (Fig. 3C, middle). Interestingly, OPCs were not differentially distributed between these two regions at P10 and P20. At P60, comparing the relative occurrences of OLs over all six tissue sections, we found MFOL1 to be more abundant in WM (CC) compared to GM (CTX) (45 ± 7%), whereas mature OL populations such as MOL1 (45 ± 4%), MOL3 (9 ± 3%), and MOL6 (52 ± 4%) were more enriched in CTX (p = 0.004, p = 0.0004, p = 0.02, p = 0.0002, respectively; n = 6 tissue sections) (Fig. 3C, right). These data suggest that myelination might be still ongoing in the CC (WM) at P60, while being more advanced in the CTX (GM), since it is more populated with mature OLs. Ongoing oligodendrogenesis and myelination can be observed until late adulthood [9,10]. , and P60 (right) (n = 4 P10 tissue sections, n = 6 P20 tissue sections, n = 6 P60 tissue sections). d Change of relative abundances of oligodendrocytes between P10 and P20 as well as P20 and P60 for CTX. e Same as (d) for CC. Statistical comparisons were determined using ANOVA and post hoc test with Tukey correction (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

Enrichment of specific OL populations at distinct areas and layers in the CTX
Specific populations of neurons and even fine transcriptomic cell types are differentially distributed in six distinct layers and areas in the CTX [12,13]. We thus investigated whether the same applies to OLs. First, we plotted the kernel density estimates for each of the twelve OL populations on top of the cell maps, highlighting the spatial preference of the OL populations ( Fig. 4A-C). The heatmaps are normalized to the highest cell density per section and reflect the relative abundances; cell maps of raw data are plotted in Additional file 3: Fig. S3. We further segmented the sections according to the Allen Brain Atlas into the primary somatosensory cortex (S1), primary motor cortex (M1), secondary motor cortex (M2), cingulate cortex area1 (Cg1), and cingulate cortex area2 (Cg2) [14].
At P20, we observed a clear transition of early myelinating populations from lateral towards medial regions of the CTX. We found NFOL1, MFOL1, and MFOL2 being more abundant in M2, Cg1, and Cg2. OPC and MOL1 presented a similar distribution in P20 as P10 being decreased in Cg2 (OPC: p = 0.004; NFOL1: p = 0.003, p = 0.002; MFOL1: p = 0.005, p = 0.005; MFOL2: p = 0.005, p = 0.008; MOL1: p = 0.002; MOL5: p = 0.002, p = 0.003, respectively; n = 6 tissue sections per condition, Fig. 4E, top). Mature OL populations were slightly more abundant, and equally distributed in the different areas, with the exception of MOL5, which was more prevalent in S1 than in other regions and less prevalent in M2 than in other regions (p = 0.002, p = 0.003, respectively). Thus, our ISS analysis indicates that the different areas in the CTX are in different stages of OL lineage progression in the postnatal and juvenile brain.
The layer analysis in P20 also showed that there are fewer OLs in layer I compared to deeper layers (I vs. II/ III: p = 0.001; I vs. IV: p = 0.0008; I vs. V: p = 0.0006; I vs. VI: p = 0.0003, n = 6 tissue sections per condition), in accordance with the literature [15,16]. In addition, there were less OLs in layer II/III and layer IV compared to layer V and layer VI, respectively (II/III vs. V: p = 0.001; II/III vs. VI: p = 0.004; IV vs. V: p = 0.03; IV vs. VI: p = 0.05, n = 6 tissue sections per condition). OLs associated with early stages of differentiation and myelination, such as NFOL1, NFOL2, MFOL1, and MFOL2 in P20 brains exposed an increased occurrence in the infragranular layers, peaking at layer V and VI (p = 0.01, p = 0.01, p = 0.02, p = 0.03, p = 0.03, respectively; n = 6 tissue sections per condition, Fig. 4E, bottom). These data could suggest that the cellular environment at layers V and VI might actively be involved with the induction of OL differentiation and myelination. Interestingly, MOL1 was also enriched in this layer at P20 (p = 0.03; n = 6 tissue sections, Fig. 4E, bottom). OL differentiation and OL stress/apoptosis are closely connected [17], and our data suggests that the local environment at the cortical layer V at P20 might be amenable for both these processes.
There is a substantial increase in the numbers of MOLs from P20 to P60. However, in contrast to P10 and P20, we did not observe any enrichment in the different cortical areas at P60 (Fig. 4F, top), suggesting that the process Fig. 4 Density of oligodendrocyte types within P10, P20, and P60 CTX (GM). a Heatmaps of the density of cells for each oligodendrocyte type in P10 CTX (GM), normalized to the highest cell density per section. The scale bar is 500 μm. b Same as (a) for P20 CTX. c Same as (a) for P60 CTX. d Top: the mean normalized occurrences of P10 oligodendrocyte populations in primary somatosensory cortex (S1), primary motor cortex (M1), secondary motor cortex (M2), cingulate cortex area1 (Cg1), and cingulate cortex area2 (Cg2) (n = 4 P10 tissue sections). Bottom: the mean cortical depth profiles with the transparent shades representing the standard error of the mean. Dashed lines denote layers I, II/III, IV, V, and VI. The X-axis (top) and Y-axis (bottom) respectively represent the occurrences after smoothening non-binned data with a kernel statistical comparisons were determined using ANOVA and post hoc test with Tukey correction (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001). e Same as (d) for P20 CTX (n = 6 P20 tissue sections). f Same as (d) for P60 CTX (n = 6 P60 tissue sections) of differentiation, myelination, and subtype specification is very advanced at this stage (although it is possible that specified MOL can still form new myelin sheaths later in late adulthood). Given that OPC and COP are being more abundant in the medial regions M2 and Cg1 at P10, while NFOL1 (M1 and S1), MFOL1 (M1), and MOL5 (S1) are more abundant in the lateral regions (Fig. 4A), a burst of differentiation and myelination might occur prior to P20 in the lateral regions. These early differences in oligodendrocyte lineage progression between the lateral and medial regions might be compensated later on within the medial regions after P20, to achieve a homogenous distribution of MOLs at P60. Similar to P20, there were less OLs in layer I compared to deeper layers (I vs. II/III: p = 0.02; I vs. IV: p = 0.02; I vs. V: p = 0.001; I vs. VI: p = 0.001, n = 6 tissue sections per condition) as well as less OLs in layer II/III and layer IV compared to layer V and layer VI, respectively (II/III vs. V: p = 0.03; II/III vs. VI: p = 0.02; IV vs. V: p = 0.05; IV vs. VI: p = 0.05, n = 6 tissue sections per condition). In P60 cortical brain sections, the early differentiating/myelinating OLs and MOL1 presented a more uniform distribution when compared to P20. Moreover, there was an increased abundance in deeper regions, particularly again at layer V/VI, of mature OL populations such as MOL3, MOL5, and MOL6 populations (p = 0.05, p = 0.008, p = 0.04, respectively; n = 6 tissue sections, Fig. 4E, bottom). Thus, our data suggest that layers V and VI are particularly enriched for myelinating OL at P20, especially in the medial regions, while at P60, these layers are enriched in specific mature OL populations.

Differential dynamics and abundances in brain and spinal cord WM and GM
Next, we investigated the occurrence of the OL populations in the P20 and P60 spinal cord. In total, we captured a tissue area of 1.1 ± 0.2 mm 2 of P20 GM and a tissue area of 1.0 ± 0.1 mm 2 of P20 WM. For P60, we captured a tissue area of 1.4 ± 0.2 mm 2 of GM and a tissue area of 1.2 ± 0.2 mm 2 of WM. The gene map of a thoracic P60 spinal cord section is shown in Fig. 5A, left, with colored symbols representing the unique transcripts, and the segmented gray matter (GM) and white matter (WM) cell maps are shown in Fig. 5A, right. Similar to the brain, we also plot kernel density estimates for each of the twelve OL populations on top of the cell maps, highlighting the spatial preference of the OL populations (Fig. 5B). Note that we split the heatmaps for GM and WM to allow direct comparisons and highlight the cellular densities in the dorsal funiculus region of WM. The cell maps of raw data are shown in Additional file 6: Fig. S6.

Oligodendroglia populations' neighboring preferences are altered from the juvenile to the adult CNS
As OL populations show different spatial preferences in the brain and spinal cord, the question arises if certain OL populations are organized in distinct neighborhoods. Therefore, we performed a nearest neighbor analysis, where each cell is represented as a node in a graph. For each node, there is a corresponding region consisting of all points of the plane closer to that node than to any other. This mathematical approach partitions the plane into regions and is called Voronoi tessellation (see shades in Fig. 7A, D). Each node is then connected to all its neighboring nodes (Delaunay triangulation), and Fig. 6 Change of relative oligodendrocyte abundances (cells/Plp1+ cells) in P20 and P60 GM and WM for brain vs. spinal cord. a Change of relative abundances of oligodendrocytes between the brain and spinal cord for P20 GM (left) and P20 WM (right) (n = 6 P20 tissue sections, n = 6 P60 tissue sections). Statistical comparisons were determined using ANOVA and post hoc test with Tukey correction (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001). b Change of relative abundances of oligodendrocytes between the brain and the spinal cord for P60 GM (left) and P20 WM (right). the distance between the points can be measured using the Euclidean distance. The neighbor with the closest Euclidean distance is the nearest neighbor. We calculated the nearest neighbor distances for each OL population for P20 GM and WM as well as P60 GM and WM and visualized the distribution of relative abundances of the nearest neighbors in heatmaps (Fig. 7B-F). Therefore, the absolute abundances of the nearest neighbors were normalized for each cell type, i.e., the relative abundances represent the percentages of instances in which a certain cell type is the nearest neighbor.
To assess the arrangement of the OLs, we applied neighborhood enrichment analysis similar to [18]. This analysis utilizes a permutation-based test and compares the cell type labels with a random configuration of labels by maintaining the positional information. In the P20 CTX, we found that NFOL2 had the highest neighboring preference with NFOL2, MFOL2 with MOL4, MOL5 with MOL5, and MOL6 with MOL3 (p = 0.05, p = 0.04, p = 0.01, p = 0.01, respectively; n = 6) (Fig. 7B, top, red boxes). These neighboring preferences were statistically different from the randomized positioning in the simulated data. In the P20 CC, we found the neighboring of MFOL1 with NFOL2, MOL2 and MOL4 with MOL2, and MOL6 with MOL4 significantly different (p = 0.04, p = 0.05, p = 0.006, p = 0.03, respectively; n = 6) (Fig. 7B, bottom, red boxes). The average Euclidean distances between OL populations, however, were not significantly different in the P20 CTX, being on average between 40 and 50 μm. In the P20 CC, the distances were even more variable, ranging from 15 μm for MOL2-MOL2 to nearly 30 μm for MFOL2-MOL3.
Interestingly, VLMCs mainly show preference to other VLMCs or to OPCs (P20 CTX VLMC-VLMC: p = 0.04, P20 CC VLMC-OPC: p = 0.01, P60 CTX VLMC-VLMC: p = 0.02, P20 GM VLMC-OPC: p = 0.05; n = 6) (Fig. 7B-F). VLMCs are associated with the vasculature and can thereby be a proxy for vasculature proximity at the CNS parenchyma. OPCs have been reported to be closely associated with the vasculature during embryonic development and during their migration [19]. While OPC migration in the adult CNS is limited, it is reactivated upon injury, alongside OPC association with the vasculature [20]. Thus, our nearest neighbor analysis supports the association of OPCs with vasculature.

Discussion
In this study, we analyzed simultaneously the expression of 124 genes expressed in different OL lineage populations in three regions of the juvenile and adult CNS. In situ sequencing (ISS) and pciSeq analysis allowed us to determine the spatial distribution of OLs with an unprecedented resolution. The pciSeq cell calling algorithm yields probabilistic readouts based on prior scRNA-seq classification. The algorithm performs well on different scRNA-seq datasets, e.g., in our previous study [21], we have successfully applied pciSeq on human cortical tissue sections, leveraging two independent scRNA-seq datasets. However, since oligodendrocytes are more in a developmental continuum compared to neurons and to minimize the risk of potential false classification, we only considered cells expressing Plp1 for our analysis, as Plp1 expression is particularly high in OL lineage cells [2] (see Fig. 1C).
In the present study, the distribution of OL lineage cells in the forebrain WM and GM at P10, P20, and P60 suggests that differentiation and myelination advance earlier in CTX (GM) than in CC (WM) in the brain. We observed an enrichment of COPs in WM at P10 as well as at P20, while more intermediate stages such as NFOL2 and MFOL1/2 were enriched in the GM (Fig. 3A). This data suggest that differentiation is actively ongoing at this stage in the CC, while in CTX, OLs are rather in a myelination phase. Indeed, using lineage tracing with the NG2creERTMBAC:ZEG mice, Zhu and colleagues showed that OL differentiation plateaus between P10 and P60 in the CTX, while it is still occurring until P60 in the CC [22]. Interestingly, these intermediate myelin-forming states are enriched in the CC WM at P60, while the CTX GM has already transitioned to fully mature OLs (Fig. 3A). This is consistent with previous lineage tracing findings with various transgenic mouse lines, indicating that OPCs in the cortical GM have a longer cell cycle time than in WM CC at P20 and P60 [23,24] and that differentiation and myelination of the CC WM can extend until at least 240 days in the mouse, but to a lesser extent in the CTX GM [22,[25][26][27]. An interesting aspect here would be to include more cell cycle/proliferation markers into future probe panel designs to further elucidate the spatial differences between OPCs and MOLs in GM and WM in terms of cell cycle time/proliferation and differentiation rates as shown by [23].
In the P60 WM and GM in the brain, the composition of mature OL populations is quite similar, with a marked increase in MOL6 and MOL5, which become the prevalent populations in these regions (Figs. 2 and  3). In the spinal cord, the spatial preferences of oligodendroglial populations are even more pronounced, which are maintained from P20 to P60 (enrichment of OPC/ COP, MOL1, and 5/6 in GM and MOL2/3 in WM). Here, it would be interesting to further investigate whether such clear spatial preference is also established later on in adulthood in the brain. Alternatively, it is possible that other regions in the brain (for instance, in ventral regions) exhibit more pronounced spatial preferences of certain OL populations.
In this study, we observed a decrease in MOL1 consistent with our previous RNASCOPE dataset, where there was a tendency for a decrease in the MOL1 population with age [4], which becomes clear with ISS. MOL1 is a population we previously observed to be dramatically increased at the injury site upon spinal cord injury [4]. The observed developmental decrease in MOL1 occurs as populations such as NFOLs, and MFOLs are also reduced, not only as a result of the dynamics of differentiation between the two time points, but also due to programmed cell death [17,[28][29][30]. The presence of MOL1 in niches during development and pathology that are associated with cell death could suggest that this mature OL state might be associated with stress and even apoptosis. Indeed, MOL1 is characterized by the expression of genes such as Egr1, Egr2, Fos, and Arc, among others [2]. Many of these genes can be upregulated upon neuronal activation but also upon stress and apoptosis [31].
However, stress and apoptosis can be induced by the processing of tissue in single-cell transcriptomic approaches [32,33]. Thus, it is challenging to determine whether stress/activated related signatures are truly biological or artifactual. Nevertheless, the detection of MOL1 in situ in the CNS by technologies such as RNASCOPE [4] and ISS (this study) strongly suggests that this population represents a physiological state of mature oligodendrocytes, although it remains unclear if it is a stress-related transient state, which might become stabilized in certain developmental or pathological niches.
One of the major advantages of technologies such as ISS is the spatial resolution, and the possibility to examine population distributions in anatomical areas. We, therefore, superimposed our data with previously defined cortical areas. Interestingly, at P10, we observed newly formed (NFOL1) and myelin-forming (MFOL1) oligodendrocytes present mainly in the deep layers of the somatosensory cortex and motor cortex 1, while at P20, we observed a higher abundance in the cingulate and motor cortex 2 areas. At P60, we observed a rather uniform distribution of the OL lineage populations. Our data, thus, indicates a lateral-to-medial gradient of OL lineage progression in the CTX and possibly a deep-tosuperficial layer gradient. While most of the mature OL lineages have been previously shown to accumulate in the deeper layers of the cortex [9,15,16], the lateral-tomedial gradient was unexpected and suggests that OL differentiation and myelination occur at different time points in different areas in the cortex. Such arealization of myelination has been shown in the mouse corpus callosum between P3 and P21, with sensory pathways being myelinated first, followed by motor pathways and association pathways [16], and also in humans [34,35]. Interestingly, specific mature OL populations have been shown to myelinate inhibitory neurons, excitatory neurons, or both [36], although it is unclear if these populations have any correspondence to the transcriptional states as MOL1-6. It is possible that neurons and other cells such as microglia present in layers V-VII might have a particular role in the differentiation, myelination, and final subtype specification of MOLs. Neuropilin-1 has been recently shown to regulate OPC proliferation and to be expressed only in activated microglia in the developing WM (CC), but not in GM (CTX) or adult CC [37]. It is also possible that distal neuronal connections modulate the timing of OL lineage progression in the cortex, as many of these cortical areas project to distinct and remote areas of the CNS.
We also observed that the maturation into MOLs is more advanced at P20 in both WM and GM of the spinal cord, when compared to the brain (Fig. 6). This is consistent with previous reports of an posterior-to-anterior gradient of myelin glycoprotein (Mog) expression in the post-natal rat CNS [38], autopsy studies in humans [35], longer cell cycle time of OPCs in the P21 and P60 spinal cord WM compared to the corpus callosum [23], and reduced capacity of differentiation of OPCs at the P30 in the spinal cord when compared to the brain [22,[25][26][27].
In our previous study, based on RNASCOPE, MOL2 cells could barely be detected in CTX and CC and were prominent in the spinal cord [4]. Using several markers for MOL2, we could confirm the prevalence of this population in the spinal cord WM versus brain WM (Fig. 6B), but could now detect a higher number of this population in the brain and an increase in these brain regions from P20 to P60 (Fig. 3B). This suggests that MOL2 might also have relevant functions not only in the spinal cord, but also in the brain, which could be particularly relevant given that we observed that MOL2 is particularly affected in the context of injury [4].
Another important aspect to consider is that our neighborhood analysis is centered on the nucleus of each cell. OPCs and OLs exhibit extensive processes, which sense the environment and neighboring cells. Thus, while the cell body of certain cells might be distant, it is possible that their processes are in close contact. Interestingly, our nearest neighbor analysis indicates that the closest neighbors of OPCs are VLMCs or OPCs themselves (Fig. 7). However, OPCs have been shown to self-avoid in the mouse adult (P60-P90) somatosensory CTX, thereby maintaining distinct territories [39]. This apparent paradox most likely reflects the high level of area coverage of OPCs in the CNS, as observed in Figs. 4 and 5. Cell death and differentiation of OPCs are compensated by proliferation and migration of neighboring OPCs, maintaining thus a constant and high density [39,40] and remodeling [41], and response to inflammatory insults in the context of injury or disease [42].

Conclusions
Our ISS dataset provides the first detailed analysis of oligodendrocyte lineage spatial organization in the mouse juvenile and adult WM and GM. While confirming previous literature findings regarding the caudal-rostral timing of myelination and association of OPCs with the vasculature, our data extends these findings by giving a detailed view of the spatial organization of distinct OL lineage populations. Furthermore, it uncovers the timing of OL differentiation and myelination in specific areas of the cortical GM and OL lineage neighborhoods. Future integration of our study with other emerging spatial transcriptomics investigations, such as targeting in detail subtypes of neurons, astrocytes, and microglia, among others, is bound to lead to a better understanding of the cellular architecture of the CNS and the principles by which it is established.

Animals
All animals (P10: Sox10-Cre/RCE, P20 and P60: Pdgfrα::CreERTM-loxP-RCE (Z/EG) [27] with mixed C57BL/6NJ and CD1 background) were free from the most common mouse viral pathogens, ectoparasites, endoparasites, and mouse bacterial pathogens harbored in research animals. The battery of screened infective agents met the standard health profile established in Karolinska Institutet animal housing facilities. Mice were kept with the following light/dark cycle: dawn 6:00-7:00, daylight 7:00-18:00, dusk 18:00-19:00, and night 19:00-6:00, and housed to a maximum number of five per cage in individually ventilated cages (IVC Sealsafe GM500, Tecniplast). Cages contained hardwood bedding (TAPVEI), nesting material, shredded paper, gnawing sticks, and a cardboard box shelter (Scanbur). The mice received a regular chew diet (either R70 diet or R34, Lantmännen Lantbruk, or CRM-P, 801722, Special Diet Services). 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. Briefly, consistent relative air humidity of 50%, 22 °C, and the air quality is controlled with the use of stand-alone air handling units supplemented with high-efficiency particulate air-filtered air. Husbandry parameters were monitored using ScanClime (Scanbur) units. Water was provided in a water bottle, which was changed weekly. Cages were changed every other week. All cage changes were done in a laminar airflow cabinet. Facility personnel wore dedicated scrubs, socks, and shoes. Respiratory masks were used when working outside of the laminar airflow cabinet.
In situ sequencing assay P10, P20, and P60 mouse brain tissues and spinal cord tissues were snap-frozen in OCT mounting medium, coronally sectioned in 10-μm-thick cryosections, and stored at − 80 °C until fixation. Sixteen coronal brain sections (bregma 1.10 mm) and 8 thoracic spinal cord sections were studied. The brain sections were hemisected. In situ sequencing (ISS), a targeted multiplexed mRNA detection assay employing padlock probes, rolling circle amplification (RCA), and barcode sequencing, was applied. Padlock probes targeting a set of 124 genes and the CartaNA NeuroKit (1010-01, CartaNA AB, now part of 10x Genomics) were used (similar to [43]). For a step-by-step protocol of the in situ sequencing chemistry, see [44]. In short, 10-μm-thick cryosections were fixed by 3% (w/v) paraformaldehyde (PFA) in DEPCtreated PBS (DEPCPBS) at room temperature (RT) for 5 min, and they were washed in DEPC-PBS followed by 0.1N HCl treatment at RT for 5 min. The sections were washed in DEPC-PBS again, and they were subjected to reverse transcription, probe ligation, and rolling circle amplification reactions. The RCA products were detected by hybridization of AF750-labeled detection oligo and decoded through four cycles of barcode sequencing, each involving a base-specific incorporation of fluorescence dyes (A = Cy5, C = Texas Red, G = Cy3, T = AF488), imaging, and removing of the incorporated dyes.

Imaging setup and image analysis
Images were acquired using a Zeiss Axio Imager Z2 epifluorescence microscope (Zeiss Oberkochen, Germany), equipped with a × 40 objective. A series of images (10% overlap between two neighboring images) at different focal depths was obtained, and the stacks of images were merged into a single image thereafter using the maximum-intensity projection (MIP) in the Zeiss ZEN software. After exporting images in .tif format, images were aligned between cycles and then stitched together using the MIST algorithm. The stitched images were then re-tiled into multiple smaller images, henceforth referred to as tiles. Tiled DAPI images were segmented with standard watershed segmentation. The tiled images with RCA products were top-hat filtered, followed by initial nonlinear image registration, spot detection, fine image registration, crosstalk compensation, gene calling, and cell calling using the pciSeq pipeline (https:// github. com/ acycl iq/ full_ coron al_ secti on [8]).

Genes and cell typing
The analysis pipeline further assigns genes to cells and then cells to cell types. It leverages the area, the global x and y coordinates, and the unique cellular identifier of the segmented cells as well as the global x and y coordinates of the genes. The assignment is done using a probabilistic framework based on single-cell RNA sequencing data (probabilistic cell typing by in situ sequencing (pciSeq)). Here, our 124 targeted genes for the ISS experiments were selected from the single-cell RNA sequencing data from [2] that was also used for cell typing. The gene panel selection algorithm and the cell typing algorithm have been described in [8] and showed that for hippocampal interneurons, around 70 genes were needed for fine classification. For the oligodendrocytes, we performed a similar analysis; however, more genes were needed for a fine classification in pciSeq to reach the same levels of accuracy. The ranking of the oligodendrocyte genes was based on the median gene expression values, i.e., the median expression values per gene per cell type from the single-cell RNA sequencing experiments from [2] were calculated. We also applied differential gene expression analysis to the single-cell RNA sequencing data using Scanpy [45]. Some genes were manually curated, analogous to [8], i.e., genes that have been used in classical immunohistochemistry experiments have been included. For example, Plp1 was added as a PAN marker gene, and we only considered cells expressing Plp1 for analysis (to minimize the risk of potential false classification). We opted to include 8-10 genes per cell type; however, since many OL-associated genes are shared between OL populations, we did not always find distinctive genes.

Data analysis
The highest probability in the pie charts was used to define the final oligodendrocyte type (winner-takes-all principle), and only cells containing Plp1 were kept for downstream analysis. Out of the 4820 OLs in P10 tissue, only 58 had a probability smaller than 50%, and only 5 OLs had a probability smaller than 40%. Out of the 8088 OLs in P20 tissue, only 117 had a probability smaller than 50%, and only 13 OLs had a probability smaller than 40%. Out of the 16,268 OLs in P60 tissue, only 109 had a probability smaller than 50%, and only 5 OLs had a probability smaller than 40%. The downstream analysis was carried out with custom-made scripts in MAT-LAB and Python (https:// github. com/ Moldia).
Comparisons of cell counts between conditions (CTX and CC or P20 and P60) were calculated by normalizing the cell counts for each sample first. We calculated both, the cell occurrences per Plp1+ cells (%) and the cell occurrences per tissue area (cells/mm 2 ). For comparisons of cell densities, we calculated the rate of change of the oligodendrocyte densities as ratio = 100/(condition X for cell_type A/condition Y for cell_type A). For comparisons of relative cell occurrences, the cell counts for each oligodendrocyte population were divided by the total number of Plp1+ cells per sample in order to account for the differences in sample size and cell densities (values between 0 and 100). The ratios of the relative cell counts were calculated as ratio = 100/(condition X for cell_type A/condition Y for cell_type A) with condition X denoting the larger cell count and condition Y denoting the smaller cell count (X, Y = {CTX, CC, P20, P60}). The ratio was subtracted from 100 to express increases, i.e., if condition X is 30 and condition Y is 10, the heatmap shows a value of 66.66%, meaning that the increase from condition Y to condition X is 20 (or 66.66% of condition X). Bar plots show the mean ± SEM (standard error of mean). Density plots were generated using the MATLAB function ksdensity and registered to Allen Mouse Brain Atlas using a custom-made script similar to [46].
Nearest neighbor analysis was performed by first partitioning the cell maps into regions using the MATLAB