Mapping Distinct Bone Marrow Niche Populations and Their Differentiation Paths

Mice

For profiling of freshly isolated tissue and of cultured samples, bone marrow was harvested from 8–16 week old adult C57BL/6J mice (Jackson Laboratories). Male and female mice were used for all experiments. All experiments were approved by the Institutional Animal Care and Use Committees of Beth Israel Deaconess Medical Center and the University of Alabama at Birmingham.

Single-cell RNA-sequencing

For scRNA-seq, we used inDrops (Klein et al., 2015) following a previously described protocol (Zilionis et al., 2017) with the following modifications: the sequence of the primer on the hydrogel beads was 5′-CGATGACGTAATACGACTCACTATAGGGTGTCGGGTGCAG[bc1,8nt]GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG[bc2,8nt]NNNNNNTTTTTTTTTTTTTTTTTTTV-3′; the sequence of the PE2-N6 primer (step 151 in (Zilionis et al., 2017) was 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGNNNNNN-3′; and the sequences of the PCR primers (steps 157 and 160 in (Zilionis et al., 2017)) were 5′-AATGATACGGCGACCACCGAGATCTACACxrefXXTCGTCGGCAGCGTC-3′ and 5′-CAAGCAGAAGACGGCATACGAGATGGGTGTCGGGTGCAG-3′. Following droplet barcoding reverse transcription, emulsions were split into aliquots of approximately 1,000 (in vivo samples) or 3,000 (cultured samples) single-cell transcriptomes and frozen at −80°C. For the in vivo samples, two libraries (n = 1,533 cells total) were prepared for mouse 1 and three libraries (n = 3,574 cells total) were prepared for mouse 2. For the three cultured samples, one library per sample was prepared (n = 2,837, n = 2,164, and n = 2,520 total cells, respectively). These cell numbers correspond to the final number of transcriptomes detected upon removal of background barcodes and stressed or dying cells (see section below).

Sequencing and read mapping

All libraries were sequenced on two runs of a NextSeq 500 (one for the in vivo samples and one for cultured samples). Raw sequencing data (FASTQ files) were processed using the previously described (Zilionis et al., 2017) inDrops.py bioinformatics pipeline (available at https://github.com/indrops/indrops), with a few modifications: Bowtie v.1.1.1 was used with parameter -e 80; all ambiguously mapped reads were excluded from analysis; and reads were aligned to the Ensembl release 85 Mus muscuius GRCm38 reference.

Cell filtering and normalization

Each library was initially filtered to include only abundant barcodes (> 500 total counts for all in vivo libraries; > 800 counts for cultured sample 1; > 1,000 counts for cultured samples 2 and 3), based on visual inspection of the histograms of total transcript counts per cell barcode. Next, we excluded putatively stressed or dying cells with > 30% (in vivo samples) or > 20% (cultured samples) of their transcripts coming from mitochondrial genes.

The gene expression counts of each cell were then normalized using a variant of total-count normalization that avoids distortion from very highly expressed genes. Specifically, we calculated ${\widehat{x}}_{i,j}$, the normalized transcript counts for gene j in cell i, from the raw counts x_i,j as follows: ${\widehat{x}}_{ij}=x_{ij}\overline{X}/X_i$, where $x_i=\sum _j{x}_{ij}$ and $\overline{X}$ the average of X_i over all cells. To prevent very highly expressed genes from correspondingly decreasing the relative expression of other genes, we excluded genes comprising > 5% of the total counts of any cell when calculating $\overline{X}$ and X_i.

In order to focus on heterogeneity within the stromal cell population, we first clustered the data and excluded hematopoietic and endothelial cell clusters based on their expression of previously known marker genes, as well as putative cell doublets. Specifically, we identified genes that were highly variable (top 25% by v-score (Klein et al., 2015), a measure of above-Poisson noise) and expressed at reasonably high levels (at least 3 counts in at least 5 cells). The counts for these genes were z-score normalized and used to perform principal components analysis (PCA), keeping the top 35 dimensions. After PCA, a k-nearest-neighbor (kNN) graph (k = 4) was constructed by connecting each cell to its four nearest neighbors, using Euclidean distance in the principal component space. Finally, we applied spectral clustering (scikit-learn SpectralClustering function with assign_labels = ‘discretize’) to the kNN graph and visualized the clustering by projecting the graph into two dimensions using a force-directed graph layout (SPRING(Weinreb et al., 2018a)).

We then identified enriched genes in each cluster and assigned cell type labels based on well-characterized cell type-specific marker genes (Figure S1D). Using this approach, we excluded putative endothelial cells, granulocytes, lymphoid progenitors, megakaryocytes, and erythroid progenitors.

For the in vivo samples, we also used Scrublet(Wolock et al., 2019) to identify two clusters of cell doublets that co-expressed marker genes of distinct cell types. 142 putative doublets were excluded.

Clustering and visualization of stromal cells

We repeated cell clustering and visualization using only the non-hematopoietic, non-endothelial clusters. Gene filtering, PCA, and kNN graph construction were performed as above, except only the top 25 principal components were used, and only seven spectral clusters were generated.

Permutation test for gene enrichment

To find significantly enriched genes in each cell cluster, we used a parameter-free permutation-based test to calculate p values, with the difference in means as the test statistic (Engblom et al., 2017). We accounted for multiple hypotheses testing with a false discovery rate of 5% using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995). To be considered for differential gene expression analysis, genes had to be expressed by at least 5 of the cells in the cluster of interest.

Gene set enrichment analysis (GSEA)

We used the online GSEA tool (http://software.broadinstitute.org/gsea/login.jsp)(Mootha et al., 2003; Subramanian et al., 2005) to find terms enriched in cluster-specific genes. As input, we used significantly enriched genes with > 2-fold higher average expression (adding a pseudocount of 0.1 transcript counts) in the cluster of interest compared to the remaining cells. The following gene set collections were tested: H (hallmark), C2 (curated), and C5 (Gene Ontology).

Population balance analysis (PBA)

The PBA algorithm calculates a scalar “potential” for each cell that is analogous to a distance, or pseudotime, from an undifferentiated source, and a vector of fate probabilities that indicate the distance to fate branch points. These fate probabilities and temporal ordering were computed using the Python implementation of PBA (available online https://github.com/AllonKleinLab/PBA), as described(Weinreb et al., 2018b).

The inputs to the PBA scripts are a set of comma-separated value (.csv) files encoding: the edge list of a kNN graph (k = 50) of the cell transcriptomes (A.csv); a vector assigning a net source/sink rate to each graph node (R.csv); and a lineage-specific binary matrix identifying the subset of graph nodes that reside at the tips of branches (S.csv). These files are provided online at http://kleintools.hms.harvard.edu/paper_websites/bone_marrow_stroma/. PBA is then run according to the following steps:

(1) Apply the script ‘compute_Linv.py-e A.csv’, here inputting edges (flag ‘-e’) from the SPRING kNN graph (see above). This step outputs the random-walk graph Laplacian, Linv.npy.

(2) Apply the script ‘compute_potential.py -L Linv.npy -R R.csv’, here inputting the inverse graph Laplacian (flag ‘-L’) computed in step (1) and the net source/sink rate to each graph node (flag ‘-R’). This step outputs a potential vector (V.npy) that is used for temporal ordering (cells ordered from high to low potential).

(3) Apply the script ‘compute_fate_probabilities.py -S S.csv -V V.npy -e A.csv -D 1 ‘, here inputting the lineage-specific exit rate matrix (flag ‘-S’), the potential (flag ‘-V’) computed in step (2), the same edges (flag ‘-e’) used in step (1) and a diffusion constant (flat ‘-D’) of 1. This step outputs fate probabilities for each cell.

Estimation of net source/sink rate vector R

A complete definition of the vector R in terms of biophysical quantities has been published previously(Weinreb et al., 2018b). We assigned negative values to R for the five cells with the highest expression of marker genes for each of the three terminal lineages. Specifically, for each lineage, we identified genes enriched in the most mature cell cluster (cluster 1 for adipocytes, cluster 6 for osteoblasts, and cluster 7 for chondrocytes), keeping genes expressed in > 25% of cells with an average expression level of > 0.5 transcript counts and > 2-fold higher average expression within the cluster than in the rest of the cells. We then identified the five cells with the highest average z-score normalized expression of these marker genes. We used the same procedure to identify ten starting cells (cells with highest score of cluster 3 [MSC] genes). We assigned different exit rates to each of the three lineages using a fitting procedure that ensured that cells identified as the putative starting MSCs would have a uniform probability to become each fate. We assigned a single positive value to all remaining cells, with the value chosen to enforce the steady-state condition $\sum _i{R}_i=0$. In the fitting procedure, all exit were initially set to one and iteratively incremented or decremented until the average fate probabilities of the putative starting MSCs were within 1 % of uniform. The separate lineage exit rates were then used to form the lineage-specific exit rate matrix S.

Extracting and ordering ceils for each lineage

To isolate the differentiation trajectory for each lineage (adipocyte, osteoblast, and chondrocyte), we ordered cells on the basis of their graph distance from the earliest predicted MSC progenitors, keeping only cells for which the probability of the given fate increased or remained constant with graph distance. Graph distance was measured by PBA potential, and starting with the cell closest to the MSC origin, we added the cell with next highest potential to the trajectory if the PBA-predicted lineage probability for cell i was at least 99.5% of the average lineage probability of the cell(s) already in the trajectory.

More formally, the procedure is as follows: order all N cells in the experiment from highest to lowest PBA potential V, with decreasing potential corresponding to increasing distance from MSCs. Let E_i be an indicator variable for the membership of ordered cell i in the erythroid trajectory (E_i = 1 if cell i is in the trajectory; otherwise, E_i = 0). If P_i is the PBA-predicted lineage probability for ordered cell i, then E_i = 1 if

Cells on a given lineage’s trajectory were then ordered by decreasing potential. Defining t_j as the index of the jth cell on a given trajectory,

Throughout this paper, we report this cell order (akin to the “pseudotime” in other publications) as a percentage of ordered cells, with the first, least differentiated cell at 0% and the most mature cell at 100%.

Significant dynamic genes

To find genes with significant changes in expression across each lineage’s cell ordering, we used a modified version of permutation test described above (see “Permutation test for gene enrichment”). Specifically, we applied a sliding window (n = 50 cells) to the cell ordering and used the difference in means between the windows with the highest and lowest expression as the test statistic, comparing the observed difference to the differences obtained after permuting the cell ordering. To be considered for analysis, genes were required to have a mean expression of at least 0.01 transcript counts in the input cells.

Dynamic gene clustering

To find groups of genes with similar expression patterns along each lineage’s differentiation ordering, we clustered the smoothed expression traces for all significantly variable genes (see previous section) with at least two-fold change between the windows with minimum and maximum expression. In detail, we smoothed the gene expression traces using a Gaussian kernel (σ = 10% of cell ordering), z-score normalized the smoothed traces, and clustered the traces using k-means clustering.

Mapping cultured cell transcriptomes to freshly isolated cell data

For Figures 5A and S5B, cells from the cultured samples were projected into the same principal component space as the in vivo data, then mapped to their most similar in vivo neighbors. In detail, counts were first converted to TPM for all samples. Then, using only the in vivo cells, the top 25% most variable genes (measured by v-score) with at least three transcript counts in at least five cells were z-score normalized and used to find the top 35 principal components. Next, the cultured cells were z-score normalized using the gene expression means and s.d. from the in vivo data and transformed into the in vivo principal component space. Lastly, each cultured cell was mapped to its closest in vivo neighbor in principal component space (Euclidean distance). In the visualization in Figure 5A, the number of cultured cells mapping to each in vivo cell was smoothed over the kNN graph (see section “Smoothing over the kNN graph”). For Figure S5B, we compared cells in the in vivo MSC cluster to the cultured cells mapping to them.

Smoothing over the kNN graph

We smoothed data over the kNN graph for visualization of the density of cultured cells mapping to in vivo cells (Figure 5A). Smoothing was performed by diffusing the number of mapped cells over the graph, as described. In brief, if G is the kNN graph, then the smoothing operator S is S = expm(- βL), where L is the Laplacian matrix of G, β is the strength of smoothing (β = 2), and expm is the matrix exponential. Then the smoothed vector X* of a vector of raw values X (number of mapped cells) is X* = SX.

RNA velocity

In order to generate the input for Velocyto (v0.17.13) (La Manno et al., 2018), which requires annotation of exons and introns for read alignments, the raw reads were reprocessed using dropEst (vO.8.5) (Petukhov et al., 2018). We first ran droptag with the default parameters, then aligned reads to the mouse genome (mm10) using STAR (v2.7.0a) (Dobin et al., 2013), allowing unique alignments only (‘-outSAMmultNmax 1’). Then dropEst was run with default settings, aside from the following: ‘-m -V -b -F -L eiEIBA’. Cell barcodes were error-corrected using the Velocyto ‘dropest-bc-correct’ command, followed by generation of Velocyto loom files using ‘run-dropest’.

Velocyto.py was run following an example notebook (https://github.com/velocyto-team/velocyto-notebooks/blob/master/python/DentateGyrus.ipynb). Briefly, the loom files generated by dropEst were merged and then filtered to include cell barcodes used in this paper’s other analyses (see ‘Cell filtering and normalization’ section). Following gene filtering (3000 most variable genes with a minimum of 3 counts and detection in 3 cells), spliced and unspliced counts were normalized separately based on total counts per cell, with a target size of the mean total counts across cells. PCA was run using 33 components, followed by KNN imputation with 66 neighbors. Gamma fitting, RNA velocity calculations, and Markov process simulations were conducted as in the Dentate Gyrus example (code for the full analysis is available at kleintools.hms.harvard.edu/paper_websites/bone_marrow_stroma).

Monocle

Monocle (v2.10.1) (Qiu et al., 2017) was run on our normalized counts matrix. Using the same gene filter as in the other analyses (see ‘Cell filtering and normalization’ section), we generated an embedding using the reduceDimension() function (‘max_components = 2, method = “DDRTree” norm_method = “none,” pseudo_expr = 0, relative_expr = FALSE, scaling = TRUE’). After generating an initial ordering with the orderCells() function, we identified the state corresponding to MSCs and re-ran orderCells() using this state as the root state. To compare the Monocle osteoblast cell ordering to that of PBA, we selected cells in the MSC and osteoblast states of the Monocle embedding and then ordered cells by Monocle pseudotime.

Stromal cell differentiation

The CD45/TER119⁻ CD31⁻ CD51⁺ Sca1^+/− or tdTOM⁺ cells were sorted and cultured for 7–14 days. MSCs were induced toward osteogenic and adipogenic lineages after plating to 70%−80% confluency. The cells were then cultured in differentiation induction media for osteoblast, adipocyte and chondrocyte differentiation. The media were changed twice a week for 21 days. For osteoblast differentiation, the cells were cultured in complete osteogenic medium from STEMCELL Technologies (05504). At day 21 after induction, the cells were fixed with ice-cold methanol and stained with 1% Alizarin Red S for 5–10 min. The excess dye was removed, the wells were dehydrated, and imaged. For adipogenesis, the cells were cultured in STEMCELL Technologies (05503). The cells were fixed with 10% formalin and stained with 0.5% Oil Red O in isopropanol. Chondrocyte differentiation was induced in monolayer culture, and the cells were cultured in 12-well plates in GIBCO (A1007101). Chondrocyte differentiation was assessed by staining with Toluidine blue, followed by removal of the excess dye and three washes with distilled water. All images were taken using a bright field microscope.