{'fraction_non_alphanumeric': 0.03459216289807032, 'fraction_numerical': 0.020030780158636203, 'mean_word_length': 4.937860255869535, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Live-animal imaging of native hematopoietic stem and progenitor cells

The biology of hematopoietic stem cells (HSCs) has predominantly been studied under transplantation conditions1,2. Particularly challenging has been the study of dynamic HSC behaviors given that live animal HSC visualization in the native niche still represents an elusive goal in the field. Here, we describe a dual genetic strategy in mice that restricts reporter labeling to a subset of the most quiescent long-term HSCs (LT-HSCs) and that is compatible with current intravital imaging approaches in the calvarial bone marrow (BM)3–5. We find that this subset of LT-HSCs resides in close proximity to both sinusoidal blood vessels and the endosteal surface. In contrast, multipotent progenitor cells (MPPs) display a broader distance distribution from the endosteum and are more likely to be associated with transition zone vessels. LT-HSCs are not found in BM niches with the deepest hypoxia and instead are found in similar hypoxic environments as MPPs. In vivo time-lapse imaging reveals that LT-HSCs display limited motility at steady-state. Following activation, LT-HSCs display heterogenous responses, with some cells becoming highly motile and a fraction of HSCs expanding clonally within spatially restricted domains. These domains have defined characteristics, as HSC expansion is found almost exclusively in a subset of BM cavities exhibiting bone-remodeling activities. In contrast, cavities with low bone-resorbing activities do not harbor expanding HSCs. These findings point to a new degree of heterogeneity within the BM microenvironment, imposed by the stages of bone turnover. Overall, our approach enables direct visualization of HSC behaviors and dissection of heterogeneity in HSC niches.

behaviors given that live animal HSC visualization in the native niche still represents an elusive goal in the field. Here, we describe a dual genetic strategy in mice that restricts reporter labeling to a subset of the most quiescent long-term HSCs (LT-HSCs) and that is compatible with current intravital imaging approaches in the calvarial bone marrow (BM) [3][4][5] . We find that this subset of LT-HSCs resides in close proximity to both sinusoidal blood vessels and the endosteal surface. In contrast, multipotent progenitor cells (MPPs) display a broader distance distribution from the endosteum and are more likely to be associated with transition zone vessels. LT-HSCs are not found in BM niches with the deepest hypoxia and instead are found in similar hypoxic environments as MPPs. In vivo time-lapse imaging reveals that LT-HSCs display limited motility at steady-state. Following activation, LT-HSCs display heterogenous responses, with some cells becoming highly motile and a fraction of HSCs expanding clonally within spatially restricted domains. These domains have defined characteristics, as HSC expansion is found almost exclusively in a subset of BM cavities exhibiting bone-remodeling activities. In contrast, cavities with low bone-resorbing activities do not harbor expanding HSCs. These findings point to a new degree of heterogeneity within the BM microenvironment, imposed by the stages of bone turnover. Overall, our approach enables direct visualization of HSC behaviors and dissection of heterogeneity in HSC niches. Current live animal HSC tracking studies require transplantation of the HSCs that are imaged, typically in the calvarium of an irradiated recipient whose BM microenvironment is severely altered 4,6 . Therefore, while engraftment biology can be studied in these models, stem cell and progenitor's behavior is likely different than in a fully unperturbed state 1,2,4 . The recent description of HSC-reporter lines in mice has facilitated the identification of these cells in bone sections and after tissue clearing; nevertheless, these reporters are still not fully HSC specific and require the use of additional markers for their identification [7][8][9] . Despite these advances there is still considerable uncertainty about the exact localization of HSC and progenitor cells. Even less is known about the nature of distinct niches that support HSC proliferation or maintain HSC quiescence 7 .

Characterization of an HSC-specific reporter mouse line
Our previous work demonstrated that the expression of the Myelodysplastic syndrome 1 (Mds1) gene is highly enriched in LT-HSCs 10 . Mds1 is transcribed from its own promoter in the MECOM locus, which also produces the well-known Evi1 gene product and the Mds1-Evi1 gene fusion product 11 . We targeted an EGFP expression cassette to the first transcriptional start site of Mds1 (Extended data figure 1a). The resulting allele is predicted to be a hypomorph for Mds1 and Mds1-Evi1 but have no effect on the expression of Evi1. Mds1 GFP/+ mice displayed normal hematopoietic parameters, HSC frequencies and cell cycle properties, and response to myelosuppression (Extended data figure 1b-f). Flow cytometric characterization of these mice confirmed the complete absence of GFP expression in any mature lineage positive (Lin+) hematopoietic cells (Extended data figure 2a and Supplementary File 2). GFP expression was predominantly restricted to a small fraction of cKit+ Sca1+ cells (Figure 1a). Using standard phenotypical parameters 12,13 , we found that 28.83 ± 11.99% of BM cells gated solely on GFP could be categorized as LT-HSCs, 26.61 ± 9.86% as short-term HSCs (ST-HSCs) and 30.12 ± 12.8% as MPPs ( Figure  1a and Supplementary File 1). Within the phenotypic LT-HSC compartment, approximately 60% of cells expressed GFP, compared to 6% of MPPs ( Figure 1b). We refer to this mixed population as hematopoietic stem and progenitor cells (HSPCs). Importantly, the reporter is not expressed in non-hematopoietic compartments of the BM (Extended data figure 2b, c and Supplementary File 3).
With the aim of eliminating the labeling of MPPs in the Mds1 GFP/+ model, we reasoned that the additional expression of a gene associated with early differentiation could facilitate exclusive LT-HSC identification. We noticed increased brightness of the reporter in phenotyical LT-HSCs, which was inversely correlated with the expression of Flt3, a gene whose expression has been associated with loss of long-term self-renewal 14,15 (Extended data figure 2d). Taking advantage of the fact that the GFP coding sequences in the Mds1 GFP/+ allele are flanked by loxP sites (Extended data figure 1a), we introduced a Flt3-Cre allele into our model (Extended data figure 3a). This allele drives Cre-recombination in cells beginning in the ST-HSC compartment 14,15 (Extended data figure 3b). Characterization of Mds1 GFP/+ Flt3-Cre mice revealed an extremely rare GFP+ population (to be referred as MFG) that corresponds to only 0.022 ± 0.013% of the lineage negative BM (Figure 1a To further validate our combined Mds1 GFP/+ Flt3-Cre model, we performed single cell RNAseq in cells isolated exclusively on the basis of GFP expression 17 . The resulting transcriptomes were then extrapolated to a published single cell transcriptional map of LT-HSCs, ST-HSCs, and multiple MPP populations (MPP2/3/4) 18 . Strikingly, virtually all MFG + transcriptomes map to the most unprimed cluster of cells where phenotypic LTHSCs also reside ( Figure 1d and Extended data figure 4a). A small fraction of MFG+ cells display megakaryocyte lineage priming (Figure 1d and Extended data figure 3c), a feature recently described in multipotent HSCs 18,19 . This analysis also highlights the efficiency of our approach in restricting GFP expression to Mds1 + Flt3cells (Extended data figure 4b). MFG cells expressed transcripts also enriched in dormant HSCs 16 (Extended data figure 4c, d). Additionally, single cell FLUIDIGM analysis for a 280-gene hematopoietic gene panel 20 , demonstrates clustering of MFG cells with LT-HSCs but no other progenitor signatures (Extended data figure 4e). Finally, we performed long-term reconstitution assays to assess the potency of MFG+ cells in comparison to cells isolated based on traditional SLAM parameters. Limiting dilution transplants using 3 to 25 cells suggest that MFG HSCs are at least as enriched as SLAM cells in transplantation capacity (Figure 1e, Extended data 4f). MFG cells also repopulated secondary recipients (Extended data figure 4g). In addition, within the SLAM compartment, long-term repopulating activity is enriched in cells expressing GFP (Extended data figure 4h). Thus, our MFG animal model allows the isolation of a highly quiescent sub-population of LT-HSCs with potent repopulation potential.

Localization of MFG-HSCs in the calvaria
Using these two reporter models, we performed live animal imaging of GFP+ cells in the calvarium 3,4 . As expected, Mds1-only HSPCs are more prevalent than MFG-HSCs ( Figure  1a, b and 2a). Both cell types were found to be in a peri-vascular location of average distance less than ~10 μm from the closest vessel ( Figure 2b, Videos 1,2). MFG-HSCs were also found in a similarly close distance to the endosteum (Figure 2c), pointing to a possible dual endosteal-vascular niche as suggested previously 3,21,22 . However, we found the MFG-HSCs to be almost exclusively associated with sinusoids rather than arterioles (Figure 2d). While HSPCs also predominantly localized close to sinusoids, a significant fraction of these were also nearby transition zone vessels ( Figure 2d) and displayed a broader distance distribution from the endosteum (Figure 2c), suggesting different micro-niches for MFG-HSC and downstream HSPCs.
Given the known developmental and structural differences between flat and long bones 23 , we also imaged femurs using a recently described quantitative deep-imaging protocol 24 . We identified a very small number of GFP+ c-Kit+ HSCs in 250 μm thick, whole-bone femoral sections of MFG mice (Extended data figure 5a-e). Approximately 70% of MFG HSCs are located within 5 μm from sinusoidal CD105-expressing cells, but are not statistically significant in comparison to random dots (Extended data figure 5f). In addition, MFG-HSCs did not significantly differ from random spots in their distance to the endosteum, (~12% were within 10 μm while >50% were more than 50 μm away, Extended data figure 5b, d, f), underscoring the difference between the calvarium and the long bone, particularly the diaphysis.

MFG-HSCs are found in similar hypoxic environments as MPPs
Low oxygen tension (hypoxia) has been historically thought to be a shared niche characteristic that is critical for maintaining stem cell quiescence 25 . However, support for the existence of a hypoxic niche has largely come from indirect evidence and measurements lacking spatial resolution 26 . Using a novel oxygen sensor and two-photon phosphorescence lifetime microscopy 27 , we measured the local pO 2 surrounding individual HSPCs and MFG-HSCs in their native microenvironment (Extended data figure 6a-f). First, we confirmed the overall hypoxic status of the calvarial BM 27 , with intravascular pO 2 in the range of 15-30 mmHg (mean ~23 mmHg, about 3% O 2 ) and extravascular pO 2 in the range of 10-25 mmHg (mean ~17 mmHg, about 2% O 2 ) (Figure 2e). We then measured pO 2 around individual HSPCs and MFG-HSCs, and found, surprisingly, similar oxygen levels (~18 and ~19 mmHg) as the average extravascular pO 2 in the BM (Figure 2e). Thus, while we detected regions in the BM with pO 2 as low as 10 mmHg, the HSPCs and MFG-HSCs were not found in these regions, suggesting that localization to regions of deepest hypoxia is not a prerequisite for maintenance of MFG-HSC quiescence.

Heterogeneous response following activation of LT-HSCs
We next examined the dynamic behaviors of HSPCs and LT-HSCs in their native niche. In vivo time-lapse imaging of the calvarium revealed that MFG-HSCs displayed low baseline motility while HSPCs showed enhanced motility (Figure 3a, b and Videos 3-6). To assess whether HSC behaviors would be affected in the context of activation, we used a cyclophosphamide/G-CSF protocol that leads to LT-HSCs expansion and subsequent mobilization 28 (Extended data figure 7a). FACS and imaging analysis demonstrated a 10fold increase in the number of MFG+ cells post-treatment ( Figure 3 c-e, video 7, and Extended data figure 7b, c). Importantly, MFG+ cells are equally enriched in the phenotypical LT-HSC fraction in this activated state (Extended data figure 7b) and display a concomitant increase in cell cycle activity (Extended data figure 7d). Time-lapse imaging of treated animals for ~6 hrs after the third G-CSF dose (Figure 3c and Videos 8) showed that, on average, MFG-HSCs became significantly more motile, with displacement measurements even higher than steady-state HSPCs (Extended data figure 7e and Video 8) but the response was heterogeneous, ranging from cells displaying limited displacement to a small number of cells that fully exited the BM into the blood stream (Video 8 and Extended data table 1). Similarly, following treatment with the myeloablative agent 5-Fluorouracil (5FU), we observed larger numbers of MFG+ cells (Extended data 8a and Videos 9, 11) and a subset exhibiting higher motility, particulary at day 20 following (Extended data figure 8b and Videos 10, 12), suggesting that enhanced motility is a common feature of the HSC response to injury, although we cannot rule out that the response is a result of indirect action on the niche by Cy/G-CSF or 5-FU.
Activated MFG-HSCs were found, on average, to be further away from the endosteum compared to native MFG-HSCs (Extended data figure 7f). Interestingly, they are found even closer to the vasculature, with an average distance of ~1 μm (Extended data figure 7g), and importantly, maintain their sinusoidal proximity (Figure 3f). By assessing the distribution of activated MFG HSCs in the entire calvarial region, unique patterns of HSC proliferation emerge. First, native MFG-HSCs are found as rare single cells within the BM, while activated MFG-HSCs appear as clusters (Figure 3c, d, e), suggesting that that MFG-HSC proliferation occurs within spatially restricted domains. Second, a subset of MFG+ cells remained as single cells while others form clusters in both the Cy/GCSF and 5FU models ( Figure 3e and Extended data figure 8c), suggesting that the proliferative response is heterogeneous among HSCs. To assess if these clusters are clonal, we generated Mds1 CreER/+ Rosa26 Confetti/+ mice 29 (Extended data figure 8d, e, f). Imaging of untreated mice revealed that labeled cells are usually found as rare single cells of different colors dispersed throughout the BM (Extended data figure 8g). Following Cy/GCSF treatment, we observed labeled cell clusters made up predominantly of cells of a single color (Extended data figure 8g, h). Quantitative analysis confirmed the increased likelihood of labeled nearby cells to be of the same color (Extended data figure 8h), providing support for clonal HSPC proliferation within confined physical domains.

Marrow cavities with bone remodeling activity correlate with HSC expansion
The observation of heterogeneous HSC proliferation in restricted physical domains prompted us to re-examine the characteristics of the microenvironment that either support clonal expansion or maintain cell quiescence. Recognizing that the bone is constantly undergoing remodeling, we hypothesized that the stages of bone turnover impose an additional degree of heterogeneity in the BM microenvironment not captured by the prevailing view centered on the endosteal vs. perivascular duality. To visualize the stages of bone turnover, we administered two (spectrally distinct) calcium binding dyes 30 48 hours apart, and imaged the calvarium immediately after the second dye injection. The two dyes mark the positions of the old and new bone fronts, respectively, and reveal where the old bone front has been eroded ( Figure 4a). We quantify the ratio of the two dyes and classify the cavities as D-type (undergoing predominantly bone deposition), R-type (predominantly bone resorption), and M-type (mixed). We confirmed that osteoblasts are biased toward Dtype cavities while osteoclasts are biased toward R-type cavities.

Discussion
Our work here describes the generation and characterization of a novel animal model where a single-color reporter can be used for the identification of LT-HSCs and their live imaging in the native niche without transplantation (Extended data table 2). We found new evidence of heterogeneity in both the HSC response to injury and the BM microenvironment, coupled to the stages of bone remodeling, that has not been recognized previously 3,27,31 . Notably, distinct cavity types are also observed in the metaphysis of lone bones (Extended data figures 11 a-f and Video 15). The existence of distinct BM cavity types implies that the traditional way of characterizing the HSC niche as endosteal or perivascular is inadequate, as the microenvironment, including the perivascular niches, contained within these cavities are likely to be different, subject to local calcium gradient 32 and downstream effects of osteoclast degradation 33 . Development of molecular profiling technology 34 capable of spatially mapping distinct bone cavities will be needed to fully characterize the regulatory factors that govern HSC quiescence versus proliferation.

Mice and genotyping
The targeting construct of Mds1 GFP/+ mice was generated by cloning and homologous recombination of the linearized targeting vector via electroporation into v6.5 ES cells. Upon selection using neomycin and clonal screening using PCR, correctly targeted ES cell clones were injected into C57Bl/6 blastocysts. Derived chimeras were initially bred to C57Bl/6 to obtain germline transmission, followed by crossing with FLPe mice 35 to remove the Frt-Neo-Frt cassette that was part of the original targeting vector. Derived mice were backcrossed onto a C57Bl/6 background for >6 generations and mice were analyzed via PCR to identify their genotype ('5-AGAGTGAAAGACCGAGTGTGTG-3', '5-GTACAGGGTAGGCTGCTCAACT-3', '5-CTCCCTCCCAGCTTTTTGCT-3'). Some of the displayed data comes from animals that still carried the Frt-Neo-Frt cassette which showed slightly lower mean fluorescence intensity in BM cells. A similar strategy was used to generate Mds1 CreER . For all experiments 2-12 months adult mice of both sexes were used while wild type littermates were used as controls. Flt3 Cre mice 14,15 , Rosa26-CAG-loxp-stoploxp-tdTomato reporter mice 36 and Rosa26-CAG-loxp-stop-loxp-Confetti reporter mice 29 were previously described. For identification of Flt3 Cre and Mds1 CreER , Cre primers were used ('5-TTACTGACCGTACACCAAAATTTGCC-3', '5-CCTGGCAGCGATCGCTATTTTCCATG-3'). Mice were bred and housed according to NIH guidelines in our AAALAC-accredited, specific-pathogen-free animal care facilities at Boston Children's Hospital or Massachusetts General Hospital. 2.3Col1-GFP mice 37 were generously provided by Dr. Jayaraj Rajagopal (Massachusetts General Hospital). All animal protocols were approved by the Animal Resources at Children's Hospital Boston, Boston Children's Hospital Institutional Animal Care and Use Committee, and Massachusetts General Hospital Institutional Animal Care and Use Committee.

HSC isolation, flow cytometry and cell sorting
BM cells were isolated by crushing of the bones using a mortar and pestle in Ca 2+ and Mg 2+ free phosphate-buffered saline (D-PBS) supplemented with 2% fetal bovine serum (FBS) and 1X penicillin/streptomycin (Pen/Strp) (Invitrogen). Viable cell number was calculated by manually counting with a hemocytometer or using a TC20 Automated Cell Counter (Bio-Rad). The cell suspension was filtered through a 70 μm strainer. For HSC identification via flow cytometry, the cells were stained for c-kit (eBioscience), Sca-1 (eBioscience, BioLegend), CD48 (BD Pharmingen) and CD150 (BioLegend) as well as a lineage marker cocktail consisted of B220, Ter-119, Gr-1, CD4 and CD8α (eBioscience). For experiments requiring lineage depletion, antibody staining for B220, Ter119, Gr-1, CD4, CD8α biotin conjugated antibodies was first performed followed by anti-biotin beads (Miltenyi Biotec) and depletion using magnetic separation columns (Miltenyi Biotec). For megakaryocyte progenitor staining the cells were stained for Lineage marker cocktail, c-kit, Sca-1, CD150, CD41 (BioLegend) and FcγR (eBioscience). For CMP/GMP/MEP identification cells were stained for Lineage marker cocktail, c-kit, Sca-1, CD34 (eBioscience) and FcγR. For mesenchymal stem cells, cell suspension was stained for the lineage cocktail, CD45, PDGFRα and Integrin-αV (eBioscience). For endothelial cells, lineage cocktail, CD45, CD31 and VE-Cadherin (eBioscience) were used. For pre/pro B, immature and mature B cell identification B220 and IgM (eBioscience) were used. For erythroid cells Ter-119 (eBioscience) was the primary marker, for monocytes and neutrophils Mac-1 (eBioscience) and Ly6-G (BD Pharmingen) were used while for T cells CD4 and CD8 were used. Antibody staining of cell suspension was always performed on ice for 45min. 4',6diamidino-2-phenylindole (DAPI, 10 μg / ml in PBS, Invitrogen) was used for exclusion of dead cells during flow cytometry. Relevant flow cytometry gating strategies for the identification of different mature cell populations, LT-HSC / MPP1 / MPP2 / MPP3/4 and endothelial cells are available in the supplementary file. Transplanted cells were double sorted to increase purity. For low cell number transplantations, all secondary sorts were performed in a plate using the automated plate reader sorting function. For FACS analysis BD LSRII Flow Cytometer was used while cell sorting was performed using a BD FACSAria II sorter (BD Biosciences). Flow cytometry data were analyzed with FlowJo (Tree Star).

Cell cycle analysis
As each animal only contains an average of 600-700 GFP+ sorted cells, the cell cycle analysis demonstrated in Figure 1c and extended data figure 3g represent GFP cells isolated from 7 Mds1 GFP/+ Flt3 Cre mice, thus the data represent an average from 7 mice that were pooled together. Upon identification and sorting purification of corresponding cellular populations as described above, cells were fixed in ice cold 70% ethanol. Cells were then washed and stained with Ki67 (BioLegend) for distinction of G0/G1 phase for 30 min on ice. DAPI was finally used for staining and analysis of G1/G2-M/S phase. Cell cycle analysis was performed using BD FACSAria II sorter.

Competitive reconstitution assays in irradiated mice and peripheral blood analysis
8-12-week-old B6.SJL-Ptprc a Pepc b /BoyJ (CD45.1) were used as BM transplant recipients. Prior transplantation mice were lethally irradiated using a gamma irradiator with a split dose of 11 Gy with 3 hr interval time between the two doses. Cells were transplanted via retroorbital injection in anesthetized mice. 100,000 whole BM CD45.1 cells were used as competitors unless stated otherwise. For limiting dilution studies, HSC frequency was calculated using Extreme Limiting Dilution Analysis software (http://bioinf.wehi.edu.au/ software/elda/ ) 38 with data at 4 months post transplant. The lower stem cell frequency reported here might be due to the incomplete backcrossing of MFG mice (6 generations), presence of constitutive CRE in hematopoietic cells and/or technical reasons. For secondary transplants, 2 million whole BM cells from primary recipients were transplanted in lethallyirradiated secondary recipients. For blood analysis of transplanted recipients, blood was collected at 4-week intervals at least for 16 weeks after transplantation. Peripheral blood was first treated with red blood cell lysis buffer to remove red blood cells followed by antibody staining for B cells (CD19, eBioscience), T cells (CD4, CD8α) and granulocytes (Ly6-G). The percentage chimerism was estimated using CD45.1 (BioLegend) and CD45.2 (eBioscience) antibody staining.

Blood cell counts, 5-fluorouracil, cyclophosphamide / GCSF and tamoxifen treatments
Blood samples were collected via the retro-orbital vein in EDTA coated tubes. Blood cell counts were performed using HEMAVET®950 (Drew Scientific) cell blood counter. For blood cell kinetic analysis upon 5-fluorouracil treatment, cell counts were performed on day 0, 3, 7, 10, 13 and 17. 5-fluorouracil was delivered via retro-orbital injection as single dose of 150 mg/kg immediately after day 0 bleeding sample was collected while control mice were injected with PBS. In addition, BM of treated mice was analyzed using flow cytometry or imaging was performed at the indicated time points. For Cyclophosphamide / GCSF experiments cyclophosphamide was delivered via intraperitoneal injection as single dose of 200 mg/kg on day 1, followed by GCSF subcutaneous injection on day 2, 3, 4 at 250 μg/kg dose per day followed by BM flow cytometry analysis or live animal imaging of the calvarial BM. For Mds1 CreER/+ Rosa26 Confetti/+ mouse experiments cyclophosphamide or PBS were administrated on day 1 followed by GCSF and tamoxifen injection on day 2 and further GCSF or PBS administration on day 3 and 4 according to each experiment. Tamoxifen was administrated at 2 mg dose via intraperitoneal injection to target labeling to the LT-HSC compartment. BM analysis of femurs and calvaria was performed 24 hrs after final GCSF dose administration according to each experiment.

Single Cell inDrops RNA sequencing
GFP cells were sorted from Mds1 GFP/+ Flt3 Cre mice and single cells were encapsulated using droplet microfluidic technology as previously described 17 . While 1200 GFP+ cells were sorted only ~400 cells were encapsulated; loss of more than half of the population is quite standard for low cell number populations in the inDrop platform. Upon library preparation of barcoded single cells, RNA sequencing was performed. To process the data, we used a previously published workflow and code available on https://github.com/ AllonKleinLab/SPRING. Ensembl release 81 mouse mm10 cDNA plus the sequence of loxP-IRES-GFP-polyA-loxP was used as reference. SPRING plots were generated using the next four-step process. Initially, the cells with few mRNA counts (< 1000 UMIs) and stressed cells (mitochondrial gene-set Z-score > 1) were filtered out. The remaining highquality cells were total-counts normalized. We next filtered genes, keeping those that were well detected (mean expression > 0.05) and highly variable (CV > 2). Finally, the data were normalized by Z-scoring each gene and applying principal components analysis (PCA), retaining the top 50 PCs. Following filtration and bioinformatics analysis, only 50 GFP+ cells passed quality control and were used for plotting. These acquired data from the GFP cells were then plotted into previously published data for LSK cells 18 upon transformation in the PCA space of the previously published data. Briefly, the two datasets were integrated using the library sklearn.decomposition.PCA (python 2.7). The fit function was used to calculate the first 50 principal components for the single cell LSK dataset 18 . Then, the normalized and filtered count matrices of the GFP and LSK were vertically combined. Before combining the two matrices, the GFP matrix was scaled in order to have a comparable amount of normalized counts in correlation to the LSK matrix. The resulting zscore combined matrix was used as input for the transform function to project the combined dataset on the original LSK dataset. The output generated by the transform function and the corresponding distance matrix that was obtained using the SPRING function get_distance_matrix, were used to generate the final SPRING plot. All data were visualized as previously described 18 . This reduced any batch effects between the two experiments. The coordinates generated by the SPRING plots were plotted using R.

Data availability
The GEO accession number for GFP cells is: GSE115908. The GEO accession number for LSK cells used for overlay was previously published (GSE90742) 18 .

Single cell fluidigm analysis
GFP cells were primarily sorted from Mds1 GFP/+ Flt3 Cre mice followed by secondary single cell sort directly in 96-well plates containing PCR buffer. Sorted plates were frozen on dry ice followed by reverse transcription, pre-amplification and high-throughput microfluidic real-time PCR for 180 transcription factors as previously described 20 . Data analysis and hierarchical clustering was performed using MultiExperiment Viewer (MeV) program. Previously published data for CMP, GMP, MEP, CLP, MPP and LT-HSC 20 using the same 180 real-time PCR platform were overlaid for comparison to the GFP cells.

Synthesis and characterization of phosphorescent Oxyphor PtG4 probe
Oxyphor PtG4 has an almost identical structure to the previously published Oxyphor PdG4 probe 39 . The synthesis of the core porphyrin as well as the synthesis of a dendritic probe similar to Oxyphor PtG4, has been published previously 40 . All synthesis steps were identical to those developed for the synthesis of PdG4 40 . MALDI-TOF was used to confirm the identity of the intermediate products and of the target probe molecule. Calibration, assessment of phosphorescence oxygen quenching and absorption spectra measurements for Oxyphor PtG4 probe were performed as previously described 39,41 .

In Vivo and Ex Vivo Imaging
All in vivo imaging experiments were performed according to procedures approved by the Institutional Animal Care and Use Committee at Massachusetts General Hospital. Briefly, mice were anesthetized with either an induction dose of 3-4% isoflurane (96% O 2 ) and maintenance dose of 1.25-2% or were given an intraperitoneal bolus injection of Ketamine (100 mg/kg) and Xylazine (15 mg/kg). Animals were deemed anesthetized by the toe pinch method. To minimize pain, mice were treated with buprenorphine (0.05-0.1 mg/kg). The hair on the calvarium was removed with scissors or mechanical trimmer and then the skin was wiped with alcohol. Next, a calvarial skin flap was created by a U-shaped incision to reveal the underlying calvarial bone as previously described 27 . Mice were injected with imaging agents (e.g. vascular labels) retro-orbitally, mounted in a custom-designed heated mouse holder, and secured to the stage of a home-built multiphoton/confocal laser-scanning video-rate microscope (for z-stack or time-lapse imaging) or an Olympus FVMPE-RS multiphoton imaging platform (for oxygenation measurements) 27 . A drop of 0.9% saline was applied to the skull to act as the immersion fluid, and a Zeiss 63× 1.15 numerical aperture water-dipping objective, an Olympus 60× 1.0 numerical aperture water-dipping objective, or an Olympus 25× 1.05 numerical aperture water-dipping objective were used for all imaging. For endpoint imaging, mice were sacrificed while under anesthesia using approved procedures. For survival imaging, the skin flap was closed with 6-0 vinyl sutures (Ethicon). Triple antibiotic ointment (Bacitracin, Neomycin, and Polymyxin-B sulfate) was applied to the top of the surgical site to minimize the chance of infection. Mice were put in a heated cage and monitored until fully awake. For 6 hr imaging sessions, mice were given an intraperitoneal injection of ~100 μL of 0.90% saline solution every hour to ensure proper hydration.
For steady-state in vivo imaging, 15-60 frames from the live video mode were averaged to acquire single 500 × 500 pixel images. Z-stacks were acquired with 1-2 μm steps and timelapse images were acquired in 30 second intervals for 20 min or longer. For Cy/GCSF in vivo imaging, z-stacks were acquired with 2 μm steps every 20 minutes for ~6 hours. Calvarial cell location maps in Figures 3d,e were created in Matlab using custom code based on the x,y,z coordinates of each cell. Data from each mouse was aligned and then overlayed based on the location of the coronal and sagittal sutures.
For pO 2 measurements, ~75 μL of 1.7 mM Pt-G4 suspended in 0.9% PBS (1X PBS, Invitrogen) was intravenously injected prior to imaging. In each pO 2 measurement location, multiple pulse excitation/emisision cycles were used to record the phophorescence lifetime decay. For each cycle, the probe was excited for ~10-20 μs followed by ~150 μs for timeresolved photon collection. Quantitative pO 2 values were obtained by fitting the phosphorescence decay with a single exponential to get an average lifetime of phosphorescence. This lifetime value was converted to pO 2 using an in vitro calibration curve for the same batch of Pt-G4.
For 5FU imaging experiments, 5FU was delivered via retro-orbital injection as single dose of 150 mg/kg as described above. In vivo BM two-photon imaging was then performed on Day 4 or Day 20 after the 5FU injection. The calvarial cell location map in Extended data figure 8c was made in similar way to the Cy/GCSF cell maps described above.
For Cyclophosphamide/GCSF ex vivo imaging, freshly fixed (4% paraformaldehyde) and excised mouse calvaria were affixed to a plastic dish, immersed in 1X PBS, and immediately imaged for 4-5 hours. Tiled z-stacks were acquired with 3 μm steps and 10% overlap between fields of view. Images were stitched together in three dimensions using Olympus FluoView software or ImageJ scripts.

Author Manuscript
Author Manuscript Author Manuscript

Author Manuscript
For examining bone remodeling activities in calvaria (in vivo) and metaphysis (ex vivo), calcium binding dyes were administered 48-hour apart via retro-orbital injection. Calvarial in vivo imaging were performed as described. Mouse tibia was freshly harvested, thinned, and imaged from the bone surface. Tiled z-stacks were acquired with 3 μm steps and stitched using ImageJ.

Image quantification
For distance measurements, the distance of each cell to blood vessels or to the nearest bone surface (i.e. endosteum identified based on SHG) was computed by hand as described previously using the Pythagorean theorem 27,42 . The bone contains an abundant amount of collagen which enables us to use SHG imaging to identify the inner bone surface (endosteum). We and others have successfully used this technique in many previous manuscripts of live BM imaging 3,27,43,44 .
The identity of blood vessels (whether arterioles, sinusoids, or transition vessels) within the calvaria were determined by a combination of morphology, location within the vessel network, location within the BM, and blood flow. Briefly, with our blood pool agent the arterioles appear as narrow (~5-10 μm diameter) and generally straight vessels with a smooth surface upstream of sinusoidal vessels, which are larger (~20 μm diameter or greater) with irregular surfaces. This definition was based on our previous work 27,44 , which confirmed that these small diameter vessels are arterioles with a faster flow speed (~2 mm/s or higher), higher pO2, and increased barrier function in comparison to sinusoidal vessels. They also stain positive for Sca1. At the transition point between arterioles and sinusoids, the vessel diameter increases. It is from this point of increase until the next vessel branching point downstream that we define as transitional vessels.
Distance measurements were performed in ImageJ v.1.51p. For display purposes, the brightness and contrast of images in the figures were adjusted, but all image analysis was performed on raw data. For motility measurements, frame-to-frame drift was corrected in three dimensions using the Template Matching plugin in ImageJ. Next, the centroid of the cell was determined for the first and last image of a 20-minute sequence, and the twodimensional displacement was calculated using the distance formula.
For cell clustering analysis, the individual tiled z-stack images were reconstructed into a single z-stack for the whole calvaria using ImageJ. Next, each cell was designated as one of three tags (red, green, or blue) based on the color of the cell during imaging, and the x, y, z coordinate was recorded by hand. Using a custom Matlab script similar to ClusterQuant 45 , we analyzed the spatial clustering (cluster size = 3) of like-colored cells in this model compared to 10,000 randomized samples to determine the statistical likelihood of the color clustering in our samples.
Graphs and statistical analyses were performed using Graph Pad Prism version 6 or higher.
The contrast and/or brightness of figure images and videos were adjusted for display purposes only.

Classification for types of BM cavities
A bone marrow cavity is defined as a 3D inclusion inside bone with a single concave endosteum (Figure 4e, extended data Figures 11d-f), while deeper down all cavities are interconnected. Once a cavity is defined using the bone SHG signal, we then proceed to classify types of BM cavities based on sequential staining of two calcium binding reagents. The first calcium binding dye (Dye1, Tetracycline or Calcein blue, Sigma, at the dose of 35mg/Kg and 30mg/Kg, respectively) was administered 48 hours before imaging to track bone resorption activities based on erosion of Dye1, and the second calcium binding dye (Dye2, Alizarin Red, 40mg/Kg) was administered 30 min before imaging to label high calcium regions, the bone fronts. The 48-hour interval was chosen based on the estimated lifespan of murine active osteoclasts 46,47 . Therefore, the double staining approach delineated approximately one bone erosion cycle occurred in the BM. As the lack of the Dye1 indicates existence of resorption while strong double staining of both dyes indicates ongoing bone deposition, the dye1/dye2 ratio contained within a single concave endosteum depicts the status of bone remodeling during the 48-hour period. For each cavity, the acquired depth covered the Dye1 and Dye2 labeled regions, typically between 80 to 120 μm beneath endosteum. For quantification of osteoblast or osteoclasts coverage (2.3Col1GFP or Cathepsin K pixels) along endosteum, z-stack double staining, Col1, and Cathepsin K images of 2 μm z-steps were rendered in 2D using maximum intensity projection in ImageJ and then analyzed (Extended data figure 9). As Cathepsin K is also expressed significantly by endothelial cells, a vascular map (Rhodamine B dextran) was acquired simultaneously and subtracted from the Cathepsin K map before retrieving the total pixel counts. For quantifying fractions of cavity types, 3D maps of calvaria were acquired and rendered in 2D using maximum intensity projection then analyzed (Figure 4f). For quantification of cavity types for figures 4g-h, total pixels of Dye1 and Dye2 were retrieved directly from the 3D stacks. Segmentation of Dye1 and Dye2 in each stack was obtained using ImageJ macros combining multiple built-in plugins. Specifically, contrast enhancement was applied consistently (0.1% saturation) for each stack. The images were smoothed using 3D image J suite plugins 48 (3D mean filter, kernel size = 1) followed by background subtraction using the rolling ball algorithm with the radius size of 100 and 250 for Dye1 and Dye2, respectively. This background subtraction step removed diffuse signals from bone autofluorescence that is more prominent in blue and green channels (Dye1) but still distinct from structured patterns of bone front staining. Segmentation was then performed using ImageJ built-in global or local thresholding algorithms to render matching binary results compared to raw stacks. The total pixel number of Dye1 and Dye2 were then obtained from the binary images to calculate Dye1 to Dye2 ratio.
For classification purposes, we defined bone cavities as (i) deposition type (D-type), Dye1/ Dye2 > 75% (ii) resorption type (R-type), Dye1/Dye2 < 25%, and (iii) mixed type (M-type), Dye1/Dye2 between 25 to 75%. Such definitions emphasize on functional perspectives of bone remodeling along endosteum (dominated by bone deposition or resorption), instead of the presence of osteoblasts or osteoclasts at the time of imaging. This is especially important given that osteoclasts went through apoptosis after each resorption cycle and may not be present at the time of imaging. Of note, bone-lining cells have been reported to occupy the bone fronts of inactive regions that lack both mature osteoblasts and calcium staining 49 . In the calvarium, as neither MDS-HSPCs nor MFG-HSCs were found in fully inactive regions, we only characterized cell distribution in D, M, and R cavities where small patches of inactive areas may be present but does not alter the distribution of the three cavity types. To quantify the number of cells per BM cavity (Figures 4g,h), Mds1 GFP/+ Flt3Cre and Mds1 GFP/+ cells were manually counted. A cell is considered belonging to a cavity if it is underneath a concave dome; there is no restriction on its distance to the endosteal surface.
Graphs and statistical analyses were performed using Graph Pad Prism version 6.

Extended Data
Christodoulou et al.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.    a, The double calcium staining strategy that identifies D, M, and R-type cavities. Dye 1, delivered 48 hours before imaging, shows the old bone front that has been eroded to varying extent, while Dye 2 delivered before imaging shows the new bone front. b-d, zoomed in regions, showing distinct cavity types defined by the Dye1 to Dye2 pixel ratios. e, A sagittal section of bone marrow cavities containing Mds1 GFP/+ cells. f, Fractions of D, M, R-type cavities in the calvaria of non-treated or treated mice (two-tailed t-test, n=155 cavities from non-treated and 80, 73 bone marrow cavities from treated animals (N=3 mice), black line represents mean, error bars represent SD). g, Quantification of MFG-HSCs in D-, M-, or Rtype cavities at steady-state and after Cy/G-CSF activation. N=4 mice per group, plotted as different symbols. Black line represents mean, error bars represent SD. h, Quantification of HSPCs in D-, M-, or R-type cavities at steady-state and after Cy/G-CSF activation. N=4 mice per group, plotted as different symbols. Two-sided Mann-Whitney test was used in all graphs unless otherwise specified, ****p<0.0001. Black line represents mean, error bars represent SD.