Combined single-cell and spatial transcriptomics reveals the molecular, cellular and spatial bone marrow niche organization

Identification and characterization of BM resident cell types by scRNAseq

Frequencies of BM cell types differ by several orders of magnitude¹¹, imposing a challenge to scRNAseq approaches. To capture both highly abundant as well as extremely rare BM resident cells, we performed droplet-based scRNAseq¹² of cells from total mouse BM, followed by progressive depletion of highly abundant cell types or enrichment of rare populations from undigested BM or enzymatically digested bones (Figure 1a). In total, our dataset comprises 7497 cells with a median detection of 1999 genes per single cell, which formed 32 clusters corresponding to distinct cell types or stages of differentiation (Figure 1b, S1). Importantly, this map is not quantitative with regard to the relative size of the different cell populations, since dissociation rates largely differ between cell types¹¹. As detailed below, cell type annotation was performed based on marker gene expression (Table S1, Figure S2,3), gene ontology analyses (Table S1), and by quantifying the enrichment of cluster gene signatures in transcriptomic data of previously described bulk populations using the CIBERSORT algorithm^13,14 (Figure S4, Supplementary Note). We used SOUP¹⁵ to confirm that the mesenchymal populations described in the following are primarily distinct cell types, while transition states between clusters exist (Figure S5).

As expected, scRNAseq of total bone marrow identified the major haematopoietic and immune cell types, including dendritic cells, neutrophils, monocytes, T cells, and distinct developmental B cell stages (Figure 1b, S2, S4a). Upon depletion of these major immune populations, scRNAseq primarily yielded erythroid progenitors that exhibited low expression of the pan-haematopoietic marker CD45 and displayed erythroid markers such as CD71¹⁶ (Figure S2). An additional 2% of these cells were non-haematopoietic (Figure 1a, S1). To efficiently capture non-haematopoietic cells in depth, we subsequently depleted cells expressing lineage markers, CD45 or CD71 from non-digested bone marrow and enzymatically digested bone chips (Figure 1a). This allowed us to interrogate rare cell populations, including neural Schwann cells (Mog, Mag), smooth muscle cells (Tagln, Acta2), putative myofibroblasts, nine different Pdgfra-positive mesenchymal cell populations and two endothelial cell clusters (Cdh5, Pecam) (Figure S3, Table S1 for a list of signature genes). The endothelial populations comprised Sca-1 (Ly6a)-expressing arterial and Emcn-expressing sinusoidal endothelial cells, respectively (Figure 2a, b)¹⁷. Among the mesenchymal cell populations, we identified chondrocytes (Acan, Sox9), osteoblasts (Osteocalcin/Bglap, Col1a1), as well as several less well described cell types. First, we observe three distinct fibroblast-like populations, which we further annotate based on their differential localization below. Second, we observe two additional populations that showed high transcriptomic similarity to previously described SCF-GFP^{+ 5} and Cxcl12-GFP^{+ 9} cells, also called Cxcl12-abundant reticular (CAR) cells (Figure 2c)^4,18. Remarkably, these populations differentially expressed adipocyte and osteo-lineage genes (Figure 2d), while exhibiting similar overall transcriptomic profiles. As a consequence, we termed these previously undescribed subpopulations Adipo- and Osteo- CAR cells, respectively. Adipo-CAR cells expressed high levels of leptin receptor (Lepr), and showed high transcriptomic similarity to previously described LepR-Cre cells^5,8,19 (Figure 2d). In contrast, Osteo- CAR cells expressed higher osterix (Sp7, Figure 2d, S3n) and lower Lepr levels. Third, we identified a cluster of Ng2 and Nestin-expressing mesenchymal cells that show transcriptomic similarity to previously described Ng2⁺Nestin⁺ mesenchymal stem cells (MSCs), which we therefore termed Ng2⁺ MSCs (Figure 2a,c). Ng2⁺Nes⁺ MSCs are clearly distinct from Nes⁺ endothelia²⁰ or smooth muscle cells²¹. Pseudotime analysis using RNA-Velocity²² placed Ng2⁺ MSCs at the apex of a differentiation hierarchy with osteoblasts, CAR cells, chondrocytes, and fibroblasts being downstream (Figure 2e).

To comprehensively cover haematopoietic stem and progenitor cell (HSPC) populations, we additionally performed scRNAseq of Lin^-cKit⁺ cells. This revealed a differentiation continuum^23–26 spanning megakaryocyte-erythrocyte and lympho-myeloid branches, as well as a separate cluster of eosinophil/basophils progenitors (Figure 1b, S2).

In order to investigate the impact of the isolation strategies of BM resident cells on the cell type recovery in scRNAseq experiments, we compared flushing of BM versus crushing of whole bones to release the marrow, and evaluated the impact of enzymatic digestion. This demonstrated that the majority of cell populations are found both in flushed and crushed BM, but many populations are only released efficiently upon intense physical treatment or enzymatical digestion (Figure S6, Figure 2f). In particular, the fibroblast populations were found more abundantly in crushed if compared to flushed bones. While we have confirmed the presence of these subpopulations in the diaphyseal BM using imaging and spatial transcriptomics (see Figure 3, 4c-f), similar cell types deriving from the cortical bone, epiphysis or periosteum might also be present in scRNAseq datasets from whole bones. Only the myofibroblast and Schwann cell populations were exclusively present in the crushed bone samples, suggesting that they are either firmly attached to the bones or may derive from distal regions, such as the periosteum or the epiphysis.

In summary, our dataset constitutes the most comprehensive scRNAseq study of the homeostatic BM to date^21,27 (Figure S6e), and spans almost all BM resident cell types described previously as well as several novel cell types. Osteoclasts, neurons and mature megakaryocytes are not covered by our dataset, likely due to cell size limitation of the scRNAseq approach. The full dataset can be interactively browsed at http://nicheview.shiny.embl.de

Spatial allocation of BM resident cell types by combined single-cell and spatial transcriptomics

While single-cell transcriptional profiling provides a powerful tool to characterize the identity and molecular makeup of BM resident cell types, information about their spatial distribution is lost in such experiments. To gain spatial information on the cell types identified above, we considered integrating our scRNAseq dataset with spatially resolved transcriptomics data of the BM. Recently described spatial transcriptomic approaches require high image and RNA quality^28–30 or rely on unfixed tissue material^31,32, and have therefore not been successfully adapted to the adult BM. Here we developed a robust and significantly improved version of laser-capture microdissection-sequencing (LCM-seq)³³, based on random priming that copes with fixed material and low input RNA quality, and is therefore compatible with fixed bone marrow sections (Figure S7a, see methods). Using this approach, we generated full-length, high quality transcriptomic data from LCM-dissected areas of fixed BM sections containing 200-300 cells in a single cell layer. We applied LCM-seq to 78 microdissected regions collected from the central diaphyseal bone marrow, based on presence or absence of sinusoidal and arteriolar blood vessels, or based on the distance from the endosteum, in order to characterize the endosteal, sub-endosteal, arteriolar, sinusoidal, and non-vascular niche composition (Figure 3a, S7b). Spatially resolved regions from the endosteum and sub-endosteum were collected solely based on the distance to the bone lining, and independent on the presence of blood vessels. Due to the high vascularization of bone marrow, ‘non-vascular’ niches can be expected to be in close proximity but not directly adjacent to sinusoids¹¹.

To evaluate our approach, we selected marker genes of cell types known to be specifically present in the respective niches, and compared their expression in scRNAseq and corresponding spatial transcriptomics data (Figure 3b, Table S2). As expected, osteoblast genes were selectively enriched at the endosteum, genes specific for sinusoidal endothelial cells were enriched at regions with high abundance of sinusoids, and arterial endothelial genes were enriched on arterioles. Marker gene sets for haematopoietic populations, Schwann cells and myofibroblasts were not significantly associated with any of the defined niches, suggesting that these cell types are either relatively evenly distributed across niches (haematopoietic populations), or insufficiently covered in the LCM-seq data (Schwann cells and myofibroblasts) (Figure 3c, Figure S7d,e). In contrast, marker gene expression of the remaining 12 BM resident populations differed significantly across analysed niches.

To systematically assess the preferential localization of these BM resident cell types to candidate niches, we computationally estimated the frequencies of cell populations defined by scRNAseq in the spatially resolved transcriptomics data using the CIBERSORT algorithm ^13,14 (Figure 3d). We extensively validated the ability of CIBERSORT to decompose bulk transcriptomes using a single-cell reference, and evaluated its performance on assembled pools of 100 cells with known composition that were processed using the LCM-seq protocol (see Supplementary Note and Figure S8). As expected, osteoblasts and chondrocytes were found exclusively at endosteal niches (Figure 3e,f). In addition, a specific fibroblast population localized preferentially to the endosteum and was therefore tentatively termed endosteal fibroblasts. In contrast, arterial endothelial cells, smooth muscle cells, and a distinct fibroblast population localized specifically to arteriolar niches (Figure 3e,f). According to its localization we tentatively termed this fibroblast population arteriolar fibroblasts. Sinusoidal endothelial cells were found in sinusoidal niches, but were also present in regions not selected based on presence of vascular subtypes (i.e. (sub)-endosteal niche), in accordance with a widespread web of sinusoids spanning the entire bone marrow (Figure 2b, S9a). Adipo-CAR cells were also found predominantly in areas with high sinusoidal occurrence, in line with their similarity to LepR-Cre cells and the reported peri-sinusoidal localization of these cell types ^4,5. In contrast, previously not described Osteo-CAR cells preferentially localized to arteriolar or non-vascular niches, suggesting that the two newly described CAR cell populations occupy distinct vascular niches. Ng2⁺ MSCs could not be unanimously placed to a niche, potentially due to additional heterogeneity of this population (Supplementary Note, Figure S8i). In summary, systematically integrating spatial transcriptomics with single-cell transcriptomics data allowed us to localize the majority of known and newly defined BM resident populations to distinct endosteal, sinusoidal, arteriolar and non-vascular niches.

Validation of cell type localization

To confirm the spatial relationships of BM cell types identified by LCM-Seq, we determined marker combinations specific to the individual populations and performed immunofluorescence imaging of bone sections.

Gene expression analyses suggested that differential expression of alkaline phosphatase (Alpl) and Cxcl12 permits the discrimination of Osteo-CAR cells (Cxcl12⁺Alpl⁺) and Adipo-CAR cells (Cxcl12⁺Alpl^neg) (Figure 4a). In contrast, osteoblasts, MSCs and arterial ECs expressed Alpl, but only low levels of Cxcl12 (Cxcl12^negAlpl⁺). To confirm the in situ localization of Adipo- and Osteo-CAR cells, we performed whole-mount immunofluorescence imaging of long bones from Cxcl12-GFP mice⁴. Co-staining with the sinusoidal maker Emcn revealed that Cxcl12⁺Alpl^neg Adipo-CAR cells in the central BM predominantly ensheathed sinuosoids, in line with results from LCM-seq (Figure 4a,b, S9a). In contrast, Cxcl12⁺Alpl⁺ Osteo-CAR cells in central BM typically showed a non-sinusoidal localization and a highly reticular morphology (Figure 4a,b S9a). By co-staining with the arterial maker Sca1, we also found numerous instances of Cxcl12⁺Alpl⁺ Osteo-CAR cells in immediate vicinity of Sca1⁺ arterioles (Figure 4c, S9b). In some instances, these formed highly reticulate structures (Figure S9b). We further noticed that GFP⁺ protrusions from Osteo-CAR cells covered Sca1⁺ arterioles, whereas endothelial cells contributed only faint levels of Cxcl12-GFP (Figure 4c, S9b). Together, these observations confirm the predominantly arteriolar and non-vascular localization of Osteo-CAR cells, as predicted from LCM-seq. Importantly, ubiquitous Alpl staining at the endosteum and sub-endosteum (possibly derived from osteoblasts) prevented analysis of CAR cell populations in these regions.

Next, we used CD31, SM22, Pdpn, Pdgfr, Col1a1 and Sca1 as markers to specifically identify smooth muscle cells (SM22⁺Pdpn^neg), fibroblast populations (Pdpn⁺Pdfgr⁺), osteoblasts (Pdpn^negCol1a1⁺), and arterial endothelial cells (Pdpn^negCD31⁺Sca1⁺) to localize them in situ (Figure 4d-f). As suggested by LCM-seq, CD31/Sca1-expressing arterioles were enveloped by SM22⁺Pdpn^neg smooth muscle cells and Pdpn⁺ arteriolar fibroblasts (Figure 4d). Arteriolar fibroblasts appear to be the cellular source for the collagen layer of the tunica externa surrounding arterioles (Figure S7c) and are likely overlapping with previously described peri-arteriolar Pdpn-expressing stromal cells³⁴. Moreover, immunofluorescence confirmed the existence of Pdpn⁺Col1a^low fibroblasts localizing to the bone-facing side of the endosteal lining made up of Pdpn^negCol1a^high osteoblasts (Figure 4e). Notably, Pdpn⁺ cells were also found at the cortical bone and periosteal regions, making it possible that fibroblast-like cells similar to those in the central marrow also derive from distal regions of the BM (i.e. cortical bone, epiphysis or periosteum) (Figure 4f).

Together, these data validate the ability of our approach to identify novel cell types, such as Adipo- and Osteo-CAR cells, and spatially allocate in the BM. Moreover, our data demonstrates that Cxcl12 is mainly synthesized at sinusoidal surfaces by Adipo-CAR cells, but also by Osteo-CAR cells at arteriolar surfaces, and in some instances in non-vascular regions (see also Figure 6g).

Spatial relationships of BM resident cell types can be accurately predicted based on scRNAseq data and cell adhesion molecule expression

How spatial relationships of cell types are established and maintained within complex organs, such as the bone marrow, remains poorly understood. It has been suggested that the expression of cell adhesion molecules represents an important mechanism that translates basic genetic information into complex three-dimensional patterns of cells within tissues³⁵, but this hypothesis has never been investigated systematically in complex systems, due to a lack of spatial and molecular tissue maps. To investigate whether cell type-specific localization of BM populations can be predicted by the differential expression of cell adhesion molecules, we compiled a comprehensive list of well-annotated cellular adhesion receptors and their cognate plasma membrane or ECM-bound ligands and developed the RNA-Magnet algorithm (Figure S10a, Table S3, Methods). RNA-Magnet predicts potential physical interactions between single cells and selected ‘attractor’ populations, based on expression patterns of cell surface receptors and their cognate surface expressed binding partners. For each single cell, RNA-Magnet provides scores for the strength of attraction (RNA-Magnet adhesiveness) and a direction, indicating the attractor population the cell is most attracted to (RNA-Magnet location). To investigate whether RNA-Magnet is capable of recapitulating the spatial relationships of BM cell types, we introduced four anchor populations representing the following niches: osteoblasts for the endosteal niche, sinusoidal endothelial cells for the sinusoidal niche, as well as arterial endothelial and smooth muscle cells both representing arteriolar niches (Figure 5a). Predicted adhesiveness of BM populations to distinct niches strongly correlated with their degree of differential localization as measured by spatial transcriptomics (Figure 5b). Strikingly, cell type-specific localization was also recapitulated with high accuracy for almost all populations (Figures 5c, S10b-d), including the differential localization of Adipo-CAR and Osteo-CAR cells to sinusoidal and arteriolar endothelia, respectively. Interestingly, smooth muscle cells were most attracted to arterial endothelial cells, whereas arteriolar fibroblasts adhered to smooth muscles, recapitulating the consecutive layering observed in blood vessels with the tunica intima (endothelial cells) surrounded by the tunica media (smooth muscle cells), and the tunica externa (ECM produced by arteriolar fibroblasts) ³⁶. Together, these observations demonstrate the ability of RNA-Magnet to predict spatial localization from single-cell gene expression data and highlight the importance of cell adhesion proteins for tissue organization in the BM.

Cellular and spatial sources of cytokines and growth factors in the BM

Key biological processes occurring in the bone marrow are thought to be mediated by the coordinated action of a diverse set of cytokines and growth factors. However, the identity of cytokine-producing cells and their organization into spatial and functional BM niches remain poorly understood. In particular, the cellular and spatial sources of haematopoietic stem cell (HSC) maintenance factors, such as Cxcl12 and Scf (Kitl), remain controversial ^6–9 with Lepr-Cre, NG2-Cre, CAR cells, Nestin-dim, osteoblasts and endothelial cells each being separately considered as potential main sources of these factors. To identify cells serving as the source of HSC maintenance factors, we initially examined our scRNAseq dataset. Our data show that Cxcl12 and Scf are indeed expressed by arterial endothelial ³⁷ and some mesenchymal cell types ⁹. However, their expression was several orders of magnitude higher in the Adipo- and Osteo-CAR populations (Figure 6a, see also Figure 4). In order to confirm Cxcl12 expression from Adipo- and Osteo-CAR cells at protein level, we developed FACS marker strategies to discriminate these cells types, and confirmed them by FACS-based index single-cell RNAsequencing (Figure 6b). Comparative analyses demonstrated that Adipo-CAR cells are CD45^negCD71^negTer119^negCD41^negCD51⁺VCAM1⁺CD200^midCD61^low, whereas Osteo-CARs and NG2⁺ MSCs expressed CD200 and CD61 at high levels (Figure 6c, S11). Intracellular flow cytometric analyses confirmed that both CAR populations are the main producers of Cxcl12 protein, whereas endothelial cells produced detectable but significantly lower levels (Figure 6d, see also Figure 4a-c). CAR cell populations were also among the main producer of key cytokines required for B cell and myeloid lineage commitment, such as Il7 and m-CSF (Csf1) (Figure 6a). Intriguingly, among all BM cell types, CAR cell populations produced the highest numbers of distinct cytokines and growth factors, and attributed the highest proportion of transcriptional activity to cytokine production, suggesting that they act as ‘professional cytokine-producing cells’ (Figure 6e,f). Together, these observations suggest a model in which differential localization of professional cytokine producing cells to cellular scaffolds results in the establishment of specific micro-niches. In line with this, spatial transcriptomics revealed that the five niches investigated showed unique production patterns of cytokines and growth factors (Figure 6g,h). Importantly, CAR cell-derived cytokines were predominantly produced in arteriolar and sinusoidal niches, and net production of growth factors and cytokines was significantly higher in vascular if compared to non-vascular niches, in line with the preferential localization of Osteo- and Adipo-CAR cells to their respective endothelial scaffolds (Figure 6g). Together, these observations suggest that specific localization of professional cytokine producing BM resident cells (such as Adipo- and Osteo-CAR cells) results in the establishment of unique niches with both arteriolar and sinusoidal niches being key sites for the production HSC maintenance and differentiation factors.

Systems-level analysis of intercellular signaling interactions of BM resident cell types

To gain a systems-level overview of potential intercellular signaling interactions, we applied RNA-Magnet to soluble signaling mediators (e.g. cytokines, growth factors etc.) and their receptors (Figures 7a, S10e). Unlike previous approaches for the reconstruction of signaling networks from single cell data^38–40, RNA-Magnet incorporates information on surface receptors with low mRNA expression; is based on a highly curated list of ligand-receptor pairs (Table S3); and specifically identifies the enrichment of signaling interactions between pairs of cell types, with an improved runtime compared to currently available packages³⁸ (see also Methods).

The network obtained from RNA-Magnet analyses formed two disconnected signaling clusters consisting of either mature immune or non-haematopoietic cells (i.e. endothelial, mesenchymal and neuronal cells), suggesting that immune and non-immune cells preferentially communicate within their respective groups. For example, many signals potentially sensed by osteoblasts⁴¹ (e.g. Bmp-, PTHrP-, FGF- signalling, see Figure S10f) are released by the newly identified endosteal fibroblasts described above. In contrast to mature immune cells, HSPC populations frequently received signals from non-haematopoietic cells, suggesting that HSPCs gradually switch from a mesenchymal to an immune signaling niche upon lineage commitment. Importantly, both CAR cell populations were the most important source of signals sensed by all myeloid and lymphoid progenitors (Figure 7a). In accordance with the specific localization of Osteo- and Adipo-CAR populations to arteriolar or sinusoidal scaffolds, an analysis of the net signaling output of distinct local niches implied that lymphoid and myeloid progenitors receive strong input from cytokines produced in vascular and especially sinusoidal niches (Figure 7b, c). Together, these analyses support a concept where distinct biological BM processes are mediated by specific combinatorial signaling input from different local niches^42,43 and provide a systems-level view of signaling interactions in the BM.

Mouse experiments

Mice were purchased from the distributors Janvier and Envigo, and housed under specific pathogen-free conditions at the central animal facility of the German Cancer Research Center. All animals used were 8-12 weeks old C56Bl/6J females. All animal experiments were performed according to protocols approved by the German authorities (Regierungspräsidium Karlsruhe).

Tissue harvesting and processing

Femurs, tibiae, hips and spines were dissected and cleaned from surrounding tissue. For all cell sorting and flow cytometry analyses, bones were crushed in cell suspension medium (RPMI 1640 (Sigma) containing 2% fetal bovine serum) using a mortar and a pestle. Dissociated cells were filtered through a 40 μm filter, spun down at 1500 rpm for 5 min, then incubated with 5 ml ACK lysis buffer (Thermo Fisher Scientific) for 5 min at room temperature for red blood cell lysis. Neutralization was achieved with 20 ml cell suspension medium. Cells were then lineage-depleted using Dynabeads Untouched Mouse CD4 Cells Kit (Thermo Fisher Scientific) according to manufacturer’s recommendations, using a home-made lineage cocktail.

For cell extraction from bones, crushed bone chips were washed four times each with 10 ml of cell suspension medium, then incubated with 10 ml of digestion medium (1 mg/ml of each Collagenase II and Dispase in HBSS, all from Gibco) for 30 min at 37 °C in a water bath. Cell suspension was then removed and filtered through a 40 μm filter and the digestion reaction was stopped by adding 40 ml of cell suspension medium. From this point on, cells were treated exactly the same as bone marrow cells above.

Flow Cytometry

Lineage-depleted bone and bone marrow cells obtained following crushing and digestion were stained with FACS antibodies (Table S4) for 30 min on ice, then washed with cell suspension medium. For intracellular CXCL12 staining, cells were fixed/permeabilized using BD Cytofix/Cytoperm (BD) for 30 min on ice, washed with BD Perm/Wash (BD) and incubated overnight with anti-CXCL12 antibody (1:10), then washed with BD Perm/Wash. All flow cytometric analyses were performed using BD Fortessa flow cytometers. Cell sorting was done using BD Aria I, Aria II and Aria Fusion sorters.

scRNA sequencing – 10x Genomics

For single-cell RNA sequencing using 10x genomics, bone and bone marrow cells were processed as described above. In addition to FACS markers, cells were stained with a DNA dye (Vybrant™ DyeCycle™ Violet, Thermo Fisher Scientific) to exclude debris and ensure that only cells are sorted for droplet-based scRNA-seq. For this purpose, cells were incubated with 2.5 μl/ml Vybrant dye in cell suspension medium for each 3x10⁶ cells at 37 °C for 30 min in a water bath. Following the incubation, cells were placed on ice and for each experiment, a total of 1.0-1.5x10⁴ events were sorted immediately into 15 μl PBS containing 2% FBS. Cell numbers were confirmed using LUNA™ Automated Cell Counter (Logos Biosystems). 33.8 μl of cell suspension were used as input without further dilution or processing, with final concentrations around 100-200 cells/μl. Reverse transcription and library construction were carried out following the Chromium Single Cell 3’ Reagent v2 protocol (10x genomics, Pleasanton, CA) according to manufacturer’s recommendations. Total cDNA synthesis was performed using 14 amplification cycles, with final cDNA yields ranging from approximately 2 ng/μl to 10 ng/μl. 10x genomics sequencing libraries were constructed as described and sequenced on an Illumina Next-Seq500, with read length 26+58 or 26+98.

FACS-indexed scRNA sequencing

Lineage-depleted bone marrow cells were obtained by crushing and stained with the following antibodies on ice for 30 min: CD41, CD45, CD51, CD61, CD71, CD200, Ter119 and VCAM1 (Table S4). Indexed single-cells were sorted into 5μl of Smart-Seq2 lysis buffer (2μM Oligo-dT30VN primer, 2mM dNTP mix (10mM each, NEB), 1:50 RNAse inhibitor (promega) and 1:125 Triton X-100 10% (Sigma-Aldrich)) and immediately snap frozen in an ethanol and dry ice bath. Plates were kept at -80 °C until processing. cDNA amplification was performed using a modified Smart-Seq2 protocol by adding, after 3 min at 72 °C, 5μl of RT mix containing 1× SMART First Strand Buffer (Clontech), 2 mM dithiothreitol (Clontech), 2 μM template switching oligo (Exiqon), 10 U μl-1 SMARTScribe (Clontech) and 10 U μl-1 RNASin plus (Promega). Transcriptome amplification was performed using 1X KAPA HiFi HS MM and 0.1μM ISPCR primer, with 21 PCR enrichment cycles. Libraries were constructed using in house produced Tn5⁵⁴ at 1:100 dilution and sequenced on an Illumina Next-Seq 500 sequencer, with 75 cycles single end sequencing.

Bone preparation and sectioning for immunofluorescent staining

Femurs were dissected and cleaned from muscle, and placed immediately in 4% PFA at 4°C for 30 min. Subsequently, femurs were washed three times in 1X PBS before incubating with 15% sucrose for 2 h, followed by 30% sucrose for another 2 hours. All incubations were performed at 4 °C. Femurs were placed in OCT and frozen at -80 °C until sectioning.

For sectioning, a cryotome (ThermoFisher) was used to generate 12 μm sections at -20 to -22 °C, which were then transferred to slides using the CryoJane Tape-Transfer System (Leica Biosystems). Slides were post-fixed with 4% PFA for 1 min to enhance the adherence of sections to the slide, then washed for 2 min by dipping in PBS. Sections were stained with antibodies in blocking buffer (PBS containing 10% goat serum and 0.2% Triton X100, see table S4 for antibodies used) at 4°C overnight, then washed by dipping into PBS for 1 min. Secondary antibody staining was done in blocking buffer at room temperature for 2 h. Sections were imaged using an LSM710 microscope (Zeiss) equipped with a polychromatic META detector (Lasers: 458, 488, 514, 561, 594, and 633nm), and an Olympus FV3000 Confocal laser scanning microscope with 4 GaAsP spectral detectors, FRAP and FRET (Lasers: 405, 445, 488, 514, 561, 594 and 640 nm). Imaris software (v8.41) was used for data analysis and representation.

Whole mount imaging of immunostained mouse femurs

The employed whole mount staining protocol has been previously described in detail¹¹. Briefly, femurs were isolated, the surrounding connective and muscle tissue carefully removed and fixed in 2% paraformaldehyde in PBS (6h, 4°C). After dehydration in 30% sucrose in PBS (72h, 4°C) the femurs were embedded in OCT medium, snap-frozen and bi-sectioned using a cryotome. The remaining thick bone marrow slice was incubated in blocking solution overnight (0.2% Triton X-100, 1 % bovine serum albumin, 10% donkey serum, in PBS) at 4°C and afterwards stained with primary and secondary antibodies (see table S4) diluted in blocking solution for 2-3 days each, including three one hour washing steps with PBS in between. The samples were immersed in RapiClear 1.52 for 6 hours to increase optical transparency. Images were acquired on an SP8 Leica confocal microscope system. Imaris software (v8.41) was used for data analysis and representation.

Bone sectioning for LCM-seq

Femurs for LCM-seq were carefully processed following guidelines of good RNA work practice. Femurs were harvested and cleaned as quickly as possible, then placed immediately in ice-cold 4% PFA, fixed for 30 min, and dehydrated in 15% and 30% sucrose solutions, prepared using RNase-free sucrose powder (Acros Organics), for 2 h each. All incubation steps were carried out on ice. Femurs were then flash-frozen in a 2-Methylbutane and dry ice bath and stored overnight at -80°C. Bone sectioning for RNA retrieval was performed in a sterilized cryotome, blades were wiped with RNaseZAP (Sigma) and slides were stored on dry ice until staining (on the same day). For staining, we adapted a shortened immunostaining protocol⁵⁵ starting with thawing the slide quickly at room temperature, then incubating with the primary antibody (1:20 to 1:40, depending on the antibody) on an aluminum rack on ice for 10 min, washing in ice-cold PBS for 30 s, then adding the secondary antibody for 5 min on ice. Higher antibody concentrations may be needed for antibodies of lower quality. Final wash and dehydration were done as follows: dipping for 30 s in each ice-cold PBS, RNase-free H₂O, 70%, 95% and 100% ethanol, in this order.

LCM-seq

Bone sections were processed using the Zeiss PALM MicroBeam Axio Observer Z1 (Zeiss), with a monochromatic Axiocam 506 mono camera, shooting laser at 355 nm and filter sets FS18 C adv (Dapi), FS44 C adv (FITC) and FS45 C adv (mCherry). Image acquisition and sample isolation were carried out using a LD Plan-NEOFLUAR 20X/0.4 objective with the adjustment ring set to 1 (Zeiss). Cutting energy was set to 45 and focus to 67, while laser pulse catapult (LPC) energy was used with delta of 12. For cutting and shooting we used the “CloseCut + AutoLPC” option. No more than 4 sections were processed and scanned in parallel, keeping collection time after staining under 30 min. A 3-channel colour image was acquired for each sample, as well as before/after LPC images and metadata including day of collection, slide number and distance of the area of interest from the bone lining. Areas of 14500 μm² on average, corresponding to around 200-300 cells, were isolated from different bone marrow districts and collected separately in 200 μl AdhesiveCap opaque Eppendorf tubes. After LPC, each collection lid was covered with 15μl of a 1:16 dilution of proteinase K in PKD buffer (Qiagen), incubated at room temperature for 5 min and subsequently snap frozen in dry ice and stored at -80°C overnight. For reverse crosslinking, samples were thawed for 5 min at room temperature, spun for 30 s to collect all the liquid at the bottom and incubated for 1h at 56°C on a PCR block⁵⁶. After incubation, samples were resuspended in 100 μl TRI Reagent (Sigma-Aldrich) under a laminar flow hood and stored in 8-PCR-strips at -80°C until RNA extraction. RNA extraction and library construction were performed in batches of 16. In brief, after thawing and spinning, each sample was transferred to a 1.5 ml Eppendorf tube and 20 μl chloroform were added. Phase separation was achieved at room temperature after vigorous shaking and spinning at 12500 rpm for 5 min, 40 μl of aqueous phase were collected from each sample and added to 75.5 μl of isopropanol and glycoblue (Invitrogen) diluted 1:150 as co-precipitant. RNA was then dehydrated at -80°C for at least 24-36 h and precipitated by centrifugation at 4°C and maximum speed. Supernatant was removed, the pellet washed once with EtOH 70%, air dried and resuspended in 8 μl nuclease-free water. Library preparation followed the SMARTER® stranded Total RNA-Seq Kit v2 – Pico Input Mammalian (Takara Bio, Japan). Due to the degraded nature of input material we omitted the fragmentation step, and 16 enrichment cycles were used in the final RNA-seq amplification. Libraries were then eluted in a minimal volume of 12 μl and sequenced on Illumina Next-Seq500, with 75bp single end reads.

Bioinformatic data analysis – Single cell data

Raw sequencing data were processed using the CellRanger pipeline (10x genomics, Pleasanton, CA), or kallisto⁵⁷ (for indexed scRNAseq). Count tables were loaded into R and further processed using the Seurat R package⁵⁸. We removed all cells with less than 500 distinct genes observed, or cells with more than 5% of UMIs stemming from mitochondrial genes. PCA was then performed on significantly variable genes, and the first 16 PCs were selected as input for clustering and t-SNE, based on manual inspection of a PC variance plot (“PC elbow plot”). Clustering was performed using the default method from the Seurat package, with resolution parameter set to 5. While lower resolution parameters caused biologically distinct groups with a low number of cells to be merged into single clusters (e.g. sinusoids and arterioles were merged into a single cluster), this relatively large parameter resulted in groups with a high number of cells to be split into an undesirable number of subgroups. We therefore computed the mean scaled gene expression values for each cluster and performed hierarchical clustering of the means using a correlation distance. Clusters with correlations of greater 0.8 were then merged together to result in the final clusters displayed in figure 1b. Marker genes for each population were identified using the FindMarkersAll function and ROC-based test statistics. For mapping of cells to a reference, the Seurat label transfer routines⁵³ were used (Figure S6), except in the context of single-cell index-sorting (Figure 6c, S11), where the cell number was insufficient for this method and scmap⁴⁹ was used instead. All mapping results were confirmed by the analysis of marker gene expression (Figure S11b and data not shown). For figure 2e, RNA velocity was run with a neighbourhood size of 50 and a linear velocity scale. The result was insensitive to the choice of parameters, except that the relative arrow lengths varied (data not shown).

Bioinformatic data analysis – LCM-seq and bulk RNA-seq data

Reads were aligned to version 38.73 of the mouse genome using STAR⁵⁹, and reads falling on exons of genes were counted using htseq-count⁶⁰. Samples containing less than one million reads on feature, as well as two outlier samples identified by PCA, were removed from the data set. Differential expression between niches was determined using the limma/voom workflow⁵⁰ while accounting for batch (i.e. slide number). Cell type proportions in the different samples were then estimated using a custom R implementation of the CIBERSORT algorithm¹³, run on raw count data. CIBERSORT requires the specification of a population-specific gene expression ‘signature matrix’; for that purpose, average gene expression profiles were computed for each cell type and 1571 population-specific genes (Figure S8a, genes were defined by specificity to a given population of 0.8 or greater, as quantified from areas under the ROC curve). To simplify analyses, we merged the highly similar HSPC subtypes into one population for CIBERSORT. Beyond the selection of genes, the CIBERSORT algorithm does not have any free parameters. Further considerations underlying the analyses using CIBERSORT, as well as a discussion on the impact of selecting different sets of genes, are detailed in the Supplementary Note.

For reanalysis of published RNA-seq data (Figure 2c, S4), count tables were created from raw sequencing data available in GEO (GSE89811 and GSE48764), or count tables were downloaded from GEO (GSE109125). For published microarray data, raw expression matrices were downloaded from GEO (GSE33158, GSE43613, GSE57729). It has previously been shown⁶¹ that CIBERSORT can be reasonably applied to decompose microarray data using a RNA-seq reference.

Data visualization

All plots were generated using the ggplot2 (v. 3.1.0) and pheatmap (v. 1.0.10) packages in R 3.4.1. Boxplots are defined as follows: The middle line corresponds to the median; lower and upper hinges correspond to first and third quartiles. The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data beyond the end of the whiskers are called “outlying” points and are plotted individually⁶².

RNA-Magnet

We developed RNA-Magnet with the goal of predicting potential physical and signalling interactions between single cells and cell populations, based on expression patterns of cell surface receptors and their cognate binding partners. Potential physical interactions were scored based on receptors that bind to surface molecules expressed on a second cell (e.g. Selectin P ligand-Selectin P, or homophilic interactions of cadherins), or based on receptors binding to structural extracellular matrix components (e.g. Integrin α1β1-Collagen). Signaling interactions were scored based on receptors binding to secreted ligands (e.g. CXCL12-CXCR4).