# 22 10 05 Combined single-cell and spatial transcriptomics reveal the molecular, cellular and spatial bone marrow niche organization

## 1. Introduction
Bone marrow (BM) niches are specialized microenvironments where distinct mesenchymal cells, the vasculature and differentiated haematopoietic cells interact to regulate the maintenance and differentiation of haematopoietic and mesenchymal stem cells (HSCs and MSCs).

- Classically, BM niches have been studied using genetic approaches that involve labelling cell types and ablating candidate niche factors based on the expression of individual reporter genes.
- However, single markers are used to define cell types in such approaches, probably resulting in the **labelling of heterogeneous populations.**
- These limitations have resulted in a controversial **debate** about the importance of **distinct cell types** and factors, and their **localization** to sinusoidal, arteriolar or endosteal niches

**Main approach**: We then used spatially resolved transcriptomics in combination with computational tools to allocate cell types to different BM niches, determine the molecular mediators of inter-cellular interactions, and identify the cellular and spatial sources of niche factors.

## 2. Results

### 2.1 Identification and characterization of BM-resident cell types by scRNAseq.
To capture both highly abundant and extremely rare BM-resident cells, we performed droplet-based scRNAseq of cells from total mouse BM, followed by progressive depletion of the abundant cell types or enrichment of rare populations from undigested BM or enzymatically digested bones

***Figure 1. a**, Overview of the FACS strategy (left). Five consecutive scRnAseq runs were performed in total. T-distributed stochastic neighbour embedding (t-SnE) of all cells, colour-coded according to the respective experiment on the left (right). **b**, 7,497 cells, which formed 32 clusters corresponding to distinct cell types or stages of differentiation. Cell-type annotation was performed based on the expression of marker genes. t-SnE projection of all cells with the clusters colour-coded.*

https://nicheview.shiny.embl.de/

- Major **haematopoietic cell types**, including dendritic cells, neutrophils, monocytes, T cells and distinct developmental stages of B cells. Erythroid progenitors that exhibited low expression of the panhaematopoietic marker CD45 and displayed erythroid markers such as CD71.
- To capture the **non-haematopoietic cells** in depth, we subsequently depleted cells expressing the lineage markers, CD45 or CD71 from non-digested BM and enzymatically digested bone chips. ->

**This is how they come up with the endothelials and mesenchymals in Figure 2.**

- **Endothelial populations** com-prised Sca1 (Ly6a)-expressing arterial ECs and Emcn-expressing sinusoidal ECs.
- **Mesenchymal lineage:** chondrocytes (Acan, Sox9), osteoblasts (Osteocalcin/Bglap, Col1a1) and several less well-described cell types:
- Three distinct fibroblast-like populations (annotated on localization).
- two additional populations that showed a high transcriptomic similarity to previously described SCF–GFP.
- Cxcl12–GFP+, Cxcl12-abundant reticular (CAR). These two populations expressed adipocyte and osteo-lineage genes differentially (Fig. 2d), despite similar overall transcriptomic profiles. We therefore termed these subpopulations ‘**Adipo-CAR**’ and ‘**Osteo-CAR**’ cells, respectively.
- Ng2- and Nestin+

***Figure 2.** **a**, Gene expression levels of Sca1 (Ly6a), endomucin (Emcn), nestin (Nes) and ng2 (Cspg4), with the relevant populations colour-coded. MAGIC44 was used for the imputation of drop-out values. **b**, Deep imaging of a BM section immunostained with antibodies against Sca1 and Emcn. Scale bar, 20 μm. **c**, Enrichment of the gene expression signatures of Adipo-CAR cells (left), Osteo-CAR cells (middle) and ng2+ MSCs (right) in previously published transcriptomes of relevant genetically labelled populations. Individual data points represent n= 3 replicates; the bars denote the mean. **d**, Boxplots of the scaled expression levels of selected genes in single cells from the Adipo-CAR (n= 528), Osteo-CAR (n= 57) and osteoblast (n= 53) populations. The mean expression of all genes annotated with the gene ontology (GO) term ‘ossification’ was computed for each cell (right). ,The endothelial populations com-prised Sca1 (Ly6a)-expressing arterial ECs and Emcn-expressing sinusoidal ECs. **e**, Projection of all mesenchymal cell types using PHATE45 with the time derivatives of the gene expression state, as determined by RnA velocity20, highlighted as arrows; n= 1,784 cells **f**, Comparison of the cell-type frequencies obtained through distinct cell isolation methods used for scRnAseq*

> **RNA velocity:** separating the different stages in mRNA transcription and maturation we can get a sense of whether the gene is being upregulated or downregulated in the next 2-3h, which corresponds to 1st derivative of gene expression. For each gene we have a direction of whether it's oing up or down. If we do it for every single gene in the cell, we estimate the "direction" of the cell. This is RNA velocity.
> 

- The **fibroblast** populations, myofibroblasts and Schwann cells were found to be more **abundant in crushed bones** when compared with flushed bones.

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

**single-cell transcriptional profiling** provides a powerful tool to characterize the **identity and molecular make-up** of BM cell types, information regarding their **spatial distribution is lost**.

- -> laser-capture **microdissection coupled with sequencing (LCM-seq):** BM sections containing 200–300 cells in a single-cell layer.
- Marker genes of cell types known to be specific to the respective niches and **compared their expression in scRNAseq and corresponding spatial transcriptomics data.**
- marker gene expression of the **remaining 12 populations differed significantly across niches**.
- -> To assess the localization of these BM cell types to candidate niches, we estimated the frequencies of the cell populations defined by scRNAseq in the spatially resolved transcriptomics data using the CIBERSORT algorithm.

***Figure 3.a,** Schematic of the experimental design. Bone sections (12 µm) were stained for arterioles (Sca1) and sinusoids (Emcn). Areas of approximately 14,500 µm2 surrounding arteries (dark yellow box), sinusoids (yellow box) and the endosteum (blue box), as well as areas with no vessels (grey box) and sub-endosteal areas (dotted blue box) were collected using laser-capture dissection and subjected to RNA-seq. A confocal image is shown for illustrative purposes. Scale bar, 100 µm.**b,** Expression of osteoblast-, sinusoid- and arteriole-specific genes in the scRNAseq data (10x) and spatial transcriptomics from different niches (LCM, LCM-seq data); n, number of cells or niches as indicated on the graph. **c,** Enrichment of population marker genes (Supplementary Table 1) among genes with differential expression between niches (Supplementary Table 2). The P values were derived using a re-sampling-based test; n = 25 cell types are shown; all HSPC subtypes were summarized as a single population. **d,** Schematic outline of the computational data-analysis strategy used. **e,** Estimated abundance of different cell types in microscopically distinct niches. The means of the CIBERSORT estimates are shown; the error bars indicate the standard error of the mean. Arteriolar niches, n = 28; endosteal niches, n = 12; sinusoidal niches, n = 14; non-vascular niches, n = 11; and sub-endosteal niches, n = 11. **f,** Heatmap summarizing the data from e.*

Using LCM-Seq they can relate the previously described populations to known niches.
- osteoblasts and chondrocytes were found exclusively at the end-osteal niches
- a specific fibroblast popula-tion localized preferentially to the endosteum -> endosteal fibroblasts
- Arterial ECs, smooth muscle cells and a distinct fibroblast population localized specifically to the arteriolar niches
- Sinusoidal ECs were found in the sinusoidal niches but were also present in the (sub)-endosteal niches, in accordance with a widespread web of sinusoids spanning the entire BM.
- Adipo-CAR cells were 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.
- Previously undescribed Osteo-CAR cells preferentially local-ized to arteriolar or non-vascular niches, suggesting that the two CAR cell populations occupy distinct niches.
- Ng2+ MSCs could not be unanimously assigned to a niche, potentially due to addi-tional heterogeneity within this population.

### 2.3 Validation of cell-type localization
**marker combinations specific to the individual populations** and performed immunofluorescence imaging of bone sections.

- alka-line phosphatase (Alpl) and Cxcl12 allows discrimination between Osteo-CAR (Cxcl12+Alpl+), Adipo-CAR cells (Cxcl12+Alpl−), and Cxcl12−Alpl+ cell types, such as osteoblasts, Ng2+ MSCs and arterial ECs.
- confirm the in situ localization of the Adipo- and Osteo-CAR cells, we performed whole-mount immunofluorescence imaging of the long bones of Cxcl12–GFP mice.
- Cxcl12+Alpl− Adipo-CAR cells in the central BM predominantly ensheathed sinusoids, similar to LCM-SEq
- Cxcl12+Alpl+ Osteo-CAR cells in the central BM typi-cally showed a non-sinusoidal localization and a highly reticular morphology.
- Cxcl12+Alpl+Sca1low Osteo-CAR cells were also frequently observed in the immediate vicinity of Sca1+GFPdim arterioles, with GFPhigh protrusions of Osteo-CAR intensively covering the arteriolar blood vessels
- **non-sinusoidal localization of Osteo-CAR cells** and qualitatively confirm that these cells localize to **both arterioles and non-vascular regions**.

***Figure 5 a,** Single-cell gene expression levels of SM22 (Tagln) and Pdpn with the relevant cluster identity colour-coded (left). Immunofluorescence staining of SM22, Pdpn and CD31 in a BM arteriole (centre). Immunofluorescent staining of Pdpn, Sca1 and Pdgfra in a BM arteriole (right). Scale bars, 20 μm. **b**, Single-cell gene expression levels of Col1a1 and Pdpnwith the relevant cluster identity colour-coded (left). Immunofluorescence staining of Col1a1, Pdpn and DAPI at the endosteal surface (right). Scale bar, 20 μm. **c**, Immunofluorescence staining of Pdpn at the endosteum, cortical bone and periosteum. Scale bar, 50 μm.*

Then they use markers for different known populations of cells:
- smooth muscle cells (SM22+Pdpn−)
- fibroblast populations (Pdpn+Pdfgr+)
- osteo-blasts (Pdpn−Col1a1+) and
- arterial ECs (Pdpn−CD31+Sca1+)

***Figure 5. a,** Single-cell gene expression levels of SM22 (Tagln) and Pdpn with the relevant cluster identity colour-coded (left). Immunofluorescence staining of SM22, Pdpn and CD31 in a BM arteriole (centre). Immunofluorescent staining of Pdpn, Sca1 and Pdgfra in a BM arteriole (right). Scale bars, 20 µm. **b,** Single-cell gene expression levels of Col1a1 and Pdpn with the relevant cluster identity colour-coded (left). Immunofluorescence staining of Col1a1, Pdpn and DAPI at the endosteal surface (right). Scale bar, 20 µm. **c,** Immunofluorescence staining of Pdpn at the endosteum, cortical bone and periosteum. Scale bar, 50 µm. All imaging experiments were repeated independently two times with similar results.*

This way they further characterize the unknown subgroups:
- As with LCM-seq data, CD31/Sca1-expressing arterioles were enveloped by SM22+Pdpn− smooth muscle cells and Pdpn+fibroblasts -> **Arteriolar fibroblasts** seem to be the **cellular source for the collagen layer** of the tunica externa **surrounding arterioles**.
- **Pdpn+Col1a-low fibroblasts localizing to the bone-facing side of the endosteal lining** made up of Pdpn−Col1ahigh osteoblasts.
- **Pdpn+cells** were also found at the **cortical** bone and periosteal regions, which raises the possibility that fibroblast-like cells similar to those in the central marrow also derive from distal regions of the BM.

This validates the ability of the approach to identify previously unknown cell types, such as **Adipo- and Osteo-CAR cells**, and spatially allocate them in the BM.

### 2.4 The spatial relationships of the BM-resident cell types can be accurately predicted based on scRNAseq data.

**spatial relationships** of cell types are established and maintained in complex organs such as the BM??

- the expression of cell adhesion molecules rep-resents an important mechanism that translates basic genetic infor-mation into complex three-dimensional patterns of cells in tissues.
- List of well-annotated cellular adhesion receptors and their cognate plasma membrane or extracellular matrix bound ligands -> RNA-Magnet algorithm to investigate whether the cell-type-specific localizations of BM populations can be predicted by the differential expression of cell adhesion molecules.
- RNA-Magnet **predicts the poten-tial physical interactions** between single cells and selected attrac-tor populations (anchors) based on the **expression patterns of the cell-surface receptors and their cognate surface-expressed binding partners**.
- RNA-Magnet adhesiveness:scores for the strength of attraction
- RNa-Magnet location: a direction indicating the attractor population the cell is most attracted to

- **four anchor populations ** representing the following niches: osteoblasts for the endosteal niche, sinusoidal ECs for the sinusoidal niche, as well as arterial ECs and smooth muscle cells to represent the arterio-lar niches.

***Figure 6. a,** t-SNE highlighting the cell type each cell is most likely to physically interact with (RNA-Magnet location, indicated by colour) and the estimated strength of adhesion (RNA-Magnet adhesiveness, indicated by opacity). n = 7,497 cells. **b,** Comparison of the estimated strength of adhesion (RNA-Magnet score) with the degree to which each cell type is differentially localized between niches (spatial transcriptomics; see also Fig. 3). n = 25 cell types are shown, all HSPC subtypes were summarized as a single population. The Pearson’s correlation (R) and its associated P value are shown. **c,** Heatmap depicting a summary of the inferred localization based on RNA-Magnet (left). The fraction of cells assigned to a certain niche is colour-coded according to the legend. Pearson’s correlation between the RNA-Magnet estimate of localization and the LCM-seq estimate of localization across n = 3 niches (right.

- **predicted adhesiveness** of BM populations to distinct niches c**orrelated strongly with their degree of differen-tial localization**, as measured by spatial transcriptomics.
- Localization was also recapitulated for all populations, including the differential localization of Adipo- and Osteo-CAR cells to the sinu-soidal and arteriolar endothelia, respectively.
- smooth muscle cells were most attracted to arterial ECs, whereas arteriolar fibroblasts adhered to smooth muscles, **recapitulating the consecu-tive layering observed in blood vessels** with the tunica intima (endo-thelial cells) surrounded by the tunica media (smooth muscle cells) and the tunica externa (extracellular matrix produced by arteriolar fibroblasts).

### 2.3 Cellular and spatial sources of cytokines and growth factors in the BM
biological processes occurring in the BM are thought to be mediated by the **coordinated action** of a diverse set of cyto-kines and growth factors. -> dentity of the cytokine-producing cells and their organization into spatial and functional BM niches remain poorly understood.

- Cxcl12 and Kitl (Scf) are indeed expressed by arterial ECs36 and some mesenchymal cell types9. However, their expression was several orders of magnitude higher in the Adipo- and Osteo-CAR populations

To confirm Cxcl12 expression from Adipo- and Osteo-CAR cells at the protein level, we developed fluorescence-activated-cell-sorting (FACS) marker strategies to discriminate these cells types and confirmed them using FACS-based index-scRNAseq

Comparative analyses demonstrated that:
- Adipo-CAR cells are CD45−CD71−Ter119−CD41−CD51+VCAM1+CD200midCD61low, whereas
- Osteo-CARs and NG2+ MSCs expressed CD200 and CD61 at high levels

**Figure 7. a**, Contribution of different cell types to distinct cytokine pools. The mean gene expression across all cells constituting each cell type is compared. **b,** Single-cell gene expression levels of Cd200 and Itgb3 (CD61) from scRNAseq data in CD45−Lin−CD71−Vcam1+ cells. MAGIC44 was used to impute the drop-out values. **c,** Surface marker levels of CD200 and CD61 from index-scRNAseq data (n = 91) in CD45−Lin−CD71−Vcam1+ cells (see Methods). The colour indicates the most similar cell type from the main dataset, as identified using scmap46. **d,** Intracellular FACS analyses of Cxcl12 expression in total BM, lineage-negative BM, Lin−CD31+ ECs and CD61/CD200 subpopulations of CD45−Lin−CD71−CD51+Vcam1+ CAR cells. Data are depicted as the mean and s.e.m. with n = 3 per group. Statistics were performed using an unpaired two-tailed t-test. The experiment was performed two times independently with a similar outcome.**e, Quantification of the number of growth factors (GFs) and cytokines expressed by each cell type and the fraction of total mRNA devoted to producing GFs and cytokines.**f,** Relative expression of cytokines and GFs in Adipo-CAR cells and fibroblasts. **g,** Expression of Cxcl12 (left), Kitl (middle) and summed expression of all cytokines and GFs (right) in the indicated niches from spatial transcriptomics data. The P values for differential expression of Cxcl12 and Kitl (left and middle) relative to non-vascular niches were computed using the moderated t-test used in limma/voom47 and corrected for multiple testing using the Benjamini–Hochberg method. Right, a two-sided Wilcoxon rank-sum test was applied. Arteriolar niches, n = 28; sinusoidal niches, n = 14; and non-vascular niches, n = 11. The boxplot elements are defined in Methods (‘Data visualization’ section).**h,** Expression of cytokines, chemokines and GFs in the different niches measured by LCM-seq. Only factors with significant differences between niches are included. *P < 0.001; **P < 0.01; *P < 0.05; and n.s., not significant.*

- Intracellular flow cytometric analyses confirmed that both **CAR populations are the main producers of Cxcl12**, whereas endothelial cells produced detectable but significantly lower levels.
- C**AR cell populations were also among the main producers of key cytokines required for B-cell and myeloid lineage commitment**, such as Il7 and m-CSF (Csf1;
- **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’

Differential localization of professional cytokine-producing cells to cellular scaffolds results in the establishment of specific micro-niches. In line with this, s**patial transcriptomics revealed that the five niches investigated displayed unique production patterns of sig-nalling mediators**

These observations suggest that the specific localization of **professional cytokine-producing BM-resident cells** results in the establishment of unique niches with both the **arteriolar and sinusoidal niches** as key sites for the production of HSC maintenance and differentiation factors.

### 2.5 Systems-level analysis of intercellular signalling interactions of BM cell types

applied R**NA-Magnet to soluble signalling mediators** (for example cytokines, growth factors and so on) and their **receptors** to gain a systems-level overview of potential intercellular signalling interactions.

- RNA-Magnet incorporates information on surface receptors with low messenger RNA expression, is based on a highly curated list of ligand–receptor pairs -> supervised

***Figure 8. a**, Inference of signalling interactions between cell types by RNA-Magnet. If a cell type is enriched in the expression of ligands for receptors expressed by a second cell type, a line is drawn between these cell types, with the colour indicating the ligand-producing cell type and the line width indicating the strength of the enrichment. See Methods for details and Extended Data Fig. 9e for a fully labelled version of the figure.**b,** Summary of the niche composition, as estimated from LCM-seq (see also Fig. 3f). **c,** Inference of the signalling interactions between niches and cell types. The line width indicates the strength of the enrichment for expression of the ligand–receptors pairs. The cell types are arranged as in b. Meg., megakaryocyte; ery., erythroid; eo., eosinophil; baso., basophil; prog., progenitor.

- **RNA-Magnet analyses formed two disconnected signalling clusters** consisting of either **mature immune or non-haematopoietic cells**, suggesting that immune and non-immune cells preferentially communicate within their respective group.
- HSPC populations frequently received signals from non-haematopoietic cells, suggesting that HSPCs gradually switch from a mesenchy-mal to an immune signalling niche following lineage commitment.
- **CAR cell populations were the most important source of the signal**s sensed by all myeloid and lymphoid progenitors.

## 3. Discussion

- single-cell and spatially resolved tran-scriptomics to transcriptionally map all major BM-resident cell types and spatially allocate them to distinct BM niches.
- cellular and spatial origins of key cytokines regulating BM haematopoiesis.
- Key HSC factors, such as Cxcl12 and Scf, are produced by two distinct and previously unappreciated subpopulations, which we have termed Osteo- and Adipo-CAR cells according to their gene expression profile.
- CAR cell subsets produce the high-est quantities of cytokines among all BM cell types, including the main cytokines mediating myeloid and lymphoid differentiation, in line with a recent study demonstrating that IL7 and Cxcl12 are produced by the same BM cell type.
- CAR cells act as ‘professional cytokine-producing cells’ and constitute cen-tral niche cells orchestrating different aspects of haematopoiesis.
- dipo-CAR cells localize to sinusoidal endothelia, Osteo-CAR cells localize to the non-vascular regions or cover arteriolar endothelia.
- arteriolar and sinusoidal vascular scaffolds represent key sites for the production of factors required for HSC maintenance and differ-entiation, with the Adipo- and Osteo-CAR cell subsets constituting central cellular hubs

# RNA Magent
https://static-content.springer.com/esm/art%3A10.1038%2Fs41556-019-0439-6/MediaObjects/41556_2019_439_MOESM1_ESM.pdf

- Surface receptors are frequently expressed lowly at the mRNA level
- Low efficiencies of reverse transcription (‘dropout’) therefore cause these genes to frequently be missed at the single cell level
- While their expression can be imputed using MAGIC, the absolute mRNA level of receptor genes is a poor proxy of absolute amounts of protein

Transformed MAGIC estimates of gene expression to a fuzzy logic variable to encode if the gene is ‘expressed’ or ‘not expressed’:

where xi,g is the MAGIC gene expression estimate of a gene g in cell i. Parameters were set to x0 = 0.5 and k = 10. Here, x0 specifies the threshold value as a fraction of maximal gene expression above which the gene is considered expressed, and k specifies the degree of fuzziness, with f(x) tending towards Boolean logic for k → ∞.

- Scoring of interaction strength between single cells and reference cell populations

RNA-Magnet provides scores of interaction strengths between a ‘sending’ cell population (i.e. a ligand-secreting cell population) and target ‘receiver’ cells. We therefore first computed population means of gene expression, and then applied a fuzzy logic AND operation to 6 determine if the ligand l is expressed by the sending population K, and the receptor r is expressed by the recipient cell c:


Where R is the set of all receptors and L(r) is the set of all ligands binding to receptor r. In the case of physical (receptor-receptor or receptor-ECM) interactions, we assume that if r is a receptor of l, l is also a receptor of r. S(c, K) was then normalised to sum to 1 across all populations included.

- highly similar cell types display differential localization -> considering the scores relative to an average score seen in similar cells: We computed a local background level of interaction scores for each cell c and sending population K by