The cellular landscape of the endochondral bone during the transition to extrauterine life

  1. Alejandro Díaz Rueda
  2. Irepan Salvador-Martínez
  3. Ismael Sospedra-Arrufat
  4. Ana Alcaina-Caro
  5. Ana Fernández-Miñán
  6. Ana M. Burgos-Ruiz
  7. Ildefonso Cases
  8. Alberto Mohedano
  9. Juan J. Tena
  10. Holger Heyn
  11. Javier Lopez-Rios
  12. Gretel Nusspaumer

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

The cellular complexity of the endochondral bone underlies its essential and pleiotropic roles during organismal life. While the adult bone has received significant attention, we still lack a deep understanding of the perinatal bone cellulome. Here, we have profiled the full composition of the murine endochondral bone at the single-cell level during the transition from fetal to newborn life and in comparison to the adult tissue, with particular emphasis on the mesenchymal compartment. The perinatal bone contains different fibroblastic clusters with blastema-like characteristics in organizing and supporting skeletogenesis, angiogenesis, and hematopoiesis. Our data also suggests dynamic inter- and intra-compartment interactions as well as a bone marrow milieu that seems prone to anti-inflammation, which we hypothesize is necessary to ensure the proper program of lymphopoiesis and the establishment of central and peripheral tolerance in early life. Our study provides an integrative roadmap for the future design of genetic and cellular functional assays to validate cellular interactions and lineage relationships within the perinatal bone.

Article activity feed

  1. Review Commons

    Note: This response was posted by the corresponding author to Review Commons. The content has not been altered except for formatting.

    Learn more at Review Commons

    Reply to the reviewers

    1. General Statements

    We would like to thank the reviewers for their critical input on the manuscript and we are glad that, overall, they recognize that the extensive analysis of the endochondral perinatal bone we describe in this manuscript can constitute a useful resource both for the bone development and hematopoietic fields. Their input has allowed us to revise the manuscript such that it is much improved in our opinion. In this section, we wish to comment on the main common aspects raised by the reviewers, while specific point-by-point responses are provided below.

    Fist, we are aware of the lack of functional assays mentioned by the reviewers, a limitation we explicitly mentioned in our original manuscript. While this is certainly a direction we will take in the future, we consider that such experiments are out of the scope and intentions of our study, given the magnitude of the resources and time they require (e.g. generation of new mouse alleles for cell fate tracking or selective ablation of specific populations, cell transplants into immunocompromised newborns, etc.). As stated in our original manuscript, this study is meant to be a resource that provides new findings and hypotheses that might be relevant for more specialized groups to functionally evaluate (e.g. teams working on thymus seeding progenitors, on adipogenesis or on immune tolerance in newborns, to name a few). As such, we believe our work has an intrinsic value. In fact, this is the first study with single cell resolution that not only compares bone populations before and after birth and with the adult tissue, but also one of the few in which all cell compartments (mesenchymal, endothelial and hematopoietic) are considered. Our manuscript hence brings a new layer of analysis not available in more directed studies, such as those based on flow cytometry (FC), in which not all populations are detected, either by lineage fraction discrimination or due to the lack of surface markers with validated antibodies for FC. This is relevant as our study identifies several new cluster-specific genetic markers and reveals their dynamic/changing expression (perinatal vs adult), or identifies that loci previously targeted for lineage tracing studies are not cluster-specific, which in our view will be useful for the interpretation of previous reports.

    The other major point brought up in the reviewers’ reports is that our analysis would be nicely complemented by the spatial localization in the perinatal bone of the various populations we describe in our study. We also agree with the reviewers on this point, which we had considered, but for which we found severe technical limitations. Spatial transcriptomics with cellular resolution would be the ideal method to address this aspect, and we tested two different methods on our samples and under several fixation and permeabilization conditions. Unfortunately, and in contrast to brain tissue used as control, these attempts have been unsuccessful in consistently detecting even ubiquitous transcripts in perinatal bone samples. As spatial transcriptomics is a technology in constant development and several new platforms and approaches are becoming available, we expect that one or several of these various methods, at the moment mostly optimized for soft tissues, will be eventually set-up for mRNA detection with true cellular resolution in perinatal and adult bone samples.

    Finally, immunofluorescence (IF) against specific markers is not a suitable approach in this case to unequivocally localize related cell populations such as the ones we describe (e.g. fibroblastic clusters). While flow cytometry has the unique advantage of performing lineage exclusion using cocktails of antibodies conjugated to the same fluorophore to label populations of cells which are not the aim of the study (e.g. hematopoietic and endothelial cells can be excluded by the use of TER119 plus CD45 and CD31, respectively), IF would require the availability of multiple specific antibodies, each conjugated to a different fluorophore, which are not available. In this regard, we would also like to point out that several studies that report the localization of specific cell populations in the bone have done so by taking advantage of genetic reporters (e.g. knock-in alleles encoding intracellular GFP or RFP). As previously mentioned, we consider that the generation of such new genetic tools is out of the scope of this manuscript.

    1. Point-by-point description of the revisions

    Reviewer #1 (Evidence, reproducibility and clarity (Required)):

    In this new manuscript, Rueda and colleagues present an extensive bioinformatics analysis of single cell transcriptomic data obtained for mouse endochondral bone cell populations before and after birth. They describe gene markers of mesenchymal and hematopoietic cells pointing to differences with adult bone populations, and they use gene ontology and trajectory analyses to infer possible roles of these cells in the developing bone. The data could provide a valuable resource for further understanding endochondral bone development and the changes driving this process in peri-natal stages. However, they are also significant weaknesses. *

    A major weakness is that the scRNA-seq data lack validation through other techniques and functional assays. Namely, in situ data are missing to locate the various cell populations in the developing bones, especially the different types of fibroblastic cells identified by the authors. Such data would go a long way to understand the possible functions of the cell populations. Although the authors tried to complement their data with a review of the literature, most of the conclusions remain purely speculative and not sufficiently supported by scientific and statistic rigor. This makes the Results section more like a discussion than a description of the results. For instance, the authors proposed important regulatory functions for the fibroblastic clusters, but there is no data supporting this other than broad GO terms associated with genes expressed in these cells. Related to this point, the title of the manuscript does not accurately reflect the content of the study**. *

    We thank the reviewer for the critical evaluation of our study and for considering it of potential interest for the field, and we have revised the manuscript to take into account his/her comments. We agree with the reviewer that including data to localize in situ the different cell populations would be highly informative. In fact, we had already attempted to perform these experiments using one of the most validated approaches, in situ sequencing (ISS). Despite assaying several fixation and permeabilization conditions, we could not reliably detect even ubiquitously expressed genes in all cells in PN1 bone sections. After discussing with a number of providers that have recently launched instruments capable of performing spatial transcriptomics technology, they all agreed that bone tissue is generally difficult to use for spatial transcriptomics technology. In summary, this data suggests that further optimization of ISS or of alternative spatial transcriptomics approaches will be needed in the future to robustly detect transcripts in bone sections with cellular resolution so as to localize in situ the various cell populations we describe in our study.

    Finally, and given our attempt to interpret our analysis of the scRNAseq data in the context of the vast literature that considers both the mesenchymal and the hematopoietic compartments, we agree with the reviewer on the speculative nature of some our conclusions that he/she mentions at the end of the paragraph, an aspect also brought up by the other reviewers. Hence, starting with the title (being now “The cellular landscape of the endochondral bone during the transition to extrauterine life”), we have systematically modified such statements throughout the text to accurately make this distinction.

    *Other points:

    1. The authors missed to report in the Results section which skeletal elements they used for their analyses and which skeletal elements were used for the adult dataset that they compared their data with. Differences in skeletal elements and in the ways whereby these samples were collected and processed could explain differences detected in the two types of datasets. Also, the sex and age of the samples for the adult dataset should be reported.*

    We now state also in the Results section that we collected forelimb long bones (excluding the handplate) for perinatal stages. In addition, we also indicate that the benchmark study by Baccin et al. used adult bone samples of mixed origin (femurs, tibiae, hips and spines) from 8-12 weeks old females. We agree with the reviewer that both this difference, as well as those related to the extraction protocols, might contribute to some of the variability we report. We now mention both these possibilities in the Discussion.

    • It is unclear whether PN1 is the day that mice are born (classically referred to as P0) or the next day.*

    As the reviewer indicates, P0 is the day of birth, and PN1 is the following day, which is the stage we chose for analysis. We have now indicated this clearly in the Materials and Methods section.

    • It is unclear whether the cells obtained for each biological replicate were pooled for the scRNA-seq assays or were treated individually. It is thus unclear how reproducible the data are.*

    In order to capture biological variability, each sample represented pooled littermates (5 fetuses for E18.5 and 4 pups at PN1), and processed as a single scRNA-seq library per stage to minimize technical variation. As our samples contained individuals from both sexes, already indicated in the original manuscript, we have now deconvoluted our datasets and computed male/female cell clustering so as to capture biological variability in duplicates (except for the sex, which is not considered as highly relevant at these stages). We assigned a “female” or “male” sex to a cell if it had at least one transcript read from a female or male specific transcript, respectively. If cells had at least one transcript read from both male and female specific genes, the cell was tagged as “undetermined”. Cells without any sex-specific transcript reads were tagged as “NA”. For the E18.5 sample we identified 21% female cells, 42% male cells, 3.7% undetermined and 33.3% NA cells. For the PN1 sample we identified 42.3% female cells, 28.1% male cells, 4.4% undetermined and 25.2% NA cells. This analysis, now shown in the new Fig. S2 and explained in Materials and Methods, reveals that all mesenchymal, hematopoietic and endothelial clusters are detected in both biological replicates. Finally, the changes we highlighted in the manuscript in the mesenchymal compartment between E18.5 and PN1 (TC, SPF and AFP) are maintained independently if the cells are processed as a single pool per stage or separated according to sex.

    • It is not clear in the gating strategy chosen for the flow cytometry as shown in Fig. 1A why the green gate containing cells expressing high levels of CD9, CD140 and CD31 has been extended in between the purple and orange gates containing CD140 and CD31 negative cells.*

    While the option mentioned by the reviewer is certainly plausible, this would have diluted the number of hematopoietic cells with intermediate CD9 levels present in our datasets. As our aim was to make sure even less abundant populations from all compartments would be captured in the scRNAseq libraries, we selected the sorting strategy depicted in Fig. 1A.

    • Are cells from all the sequenced samples homogenously distributed in the scRNA-seq clusters? Authors should provide this information and add statistic when they describe changes in the amount of cells per cluster between E18.5 and PN1 stages.*

    As mentioned in comment 3, we have now deconvoluted the datasets according to sex, which shows all clusters are represented in both biological duplicates and that overall follow similar trends in the E18.5 and PN1 samples.

    • On the basis of what markers the AFP population has been called adipogenic? Authors present Ptch2 and Notch3 as markers of this cluster, but not adipogenic progenitor genes.*

    Fig. 1C represents differentially expressed markers between clusters, which is why we chose these two representative markers for the AFP population. AFP cells also express adipogenic genes such as Pparg, Lpl or Gas6, although not exclusively. Cluster annotation is based on their molecular signature per se, GO and SCENIC analysis, which identified adipogenic regulons as active in the AFP cluster (see Fig. S10).

    • Authors claim that there is a good correlation between OsC and osteo-CAR clusters. However, OsC cells do not express Cxcl12 and other typical CAR cell markers.*

    We thank the reviewer for raising this very important point, as both our study and other recent ones (Liu et al., 2022, doi.org/10.1038/s41467-022-28775-x; Kara et al., 2023, doi:10.1016/j.devcel.2023.02.003), show that the most representative genes that historically define CARs (e.g. Cxcl12-high and LepR) are still not expressed at these stages, which indicates that these cells are not yet present at perinatal stages. Accordingly, we did not annotate any perinatal cluster as CAR cells. However, we did observe that other genes such as Runx2, Sp7, Spp1 or Alpl define populations belonging to the OsC cluster that map to the same integrated coordinates as the adult osteo-CAR cluster defined by Baccin et al. (Fig. 2C, bottom panel and Fig. S3, bottom panels). These observations stress the importance of performing ontogenic analysis for each marker defining specific populations, and that data obtained from adult tissue cannot be extrapolated to perinatal stages. We have also corrected the figure legend, which was certainly confusing in this respect.

    • In Figure 6 expression of PaS cell markers should be shown for both adult and perinatal populations. Additionally, have the authors tested that the sorted cells in panel C have the same progenitor properties as the PaS cells?*

    As requested by the reviewer, we have added the expression of PaS cell markers in adults to Fig. 6 (new panels in Fig. 6B). We are certainly considering exploring in the future the progenitor properties of the sorted cells in comparison to PaS, but these in vivo experiments will require extensive experimentation such as kidney subcapsular transplants in newborns in an immunocompromised background. We consider that these complex in vivo experiments are out of the scope of this manuscript, conceived as a resource paper.

    Reviewer #1 (Significance (Required)):

    *As indicated in the comments for the authors, the new scRNA-seq data could become a useful resource for subsequent studies, but they are at present insufficient to represent a significant scientific advancement. The main concern is that new cell populations appear to have been identified by the authors, but a number of questions were not answered such as regarding their actual location in the skeletal elements, their origins, their fates and their functions. Generating such data would require a major amount of effort and require substantial revision of the manuscript.

    Our study uncovers, in an unbiased and unsupervised manner, the heterogeneity of the entire perinatal bone with cellular resolution. As the reviewer points out, addressing the origin, fates and functions of the various cell clusters we describe would require a major financial effort and years to be completed. We consider that those aims are well beyond the aims of our manuscript, which is intended as a resource for the large scientific community in the field.

    *Reviewer #2 (Evidence, reproducibility and clarity (Required)): *

    **In the study by Rueda et al. the authors use single cell RNA-sequencing to investigate differences in cellular composition bones/bone marrow between late gestational stage mouse embryos and their perinatal counterparts. The authors describe specific differences in the relative abundance of putative cell types and use established bioinformatic tools to infer interactions as well as molecular mechanisms determining specific functions. The employed methods are well described and the results are presented in a very clear and understandable manner. Despite that, the findings do not provide any substantial knowledge advance but are rather confirming work of published literature while supplementing available single cell RNA-sequencing datasets of mouse bones at adult ages. As such, this work provides an interesting resource but does not report novel biology.

    Major comment: The authors explore the interesting transition of embryonic to perinatal bone/bone marrow using single RNA-sequencing. This fills a gap for the field of bone and hematopoietic researchers. There is little to criticize about the presented data. However, while it provides a nice resource, the knowledge gained is incremental. As acknowledged by the authors themselves, their study lacks functional validation of any findings made or conclusions drawn from bioinformatic tools in this manuscript. They use published work to validate their findings but do not go beyond that to confirm putative new biology. Some examples are listed in the minor comments.*

    We thank the reviewer for the overall positive comments on our manuscript as a resource study and his/her critical input that we have taken into account when preparing the revised version of the manuscript. Despite the lack of functional validation (already discussed in the General Statements section), we feel that our molecular analysis does provide new valuable insight into the biology of the perinatal bone. For instance, this is the first report that categorizes the heterogeneity of all perinatal bone populations with single cell resolution, and the first that explores the cellular changes in the bone that accompany birth. It also provides an important resource for the generation of more specific genetic models for cell fate tracking or for the interpretation of previous results. Finally, while it is only an inference, our interactome analysis predicts interactions between specific mesenchymal and hematopoietic populations, opening new possibilities for specialists in the specific fields to functionally address in a directed manner (e.g. interactions between the mesenchymal compartment and the Eo/Bas or the ICL-TSP2 subpopulations, which, to the best of our knowledge, have not been previously postulated).

    *Minor comments: *

    Remark: 10X Chromium does not provide whole transcriptomic coverage but rather captures the most highly abundant transcripts without for example being able to distinguish alternatively spliced gene variants. Based on that, interpretation of gene expression, or the absence of a gene in the dataset, should be interpreted carefully.*

    The reviewer is correct, as the method used only captures the 3’UTR of each transcript. We have therefore added a sentence in Materials and Methods to address the limitation of the method. Still, our approach is widely used in the field, as it allows capturing several thousand cells and one facilitating the direct comparison with other datasets, as we ourselves did when integrating the adult dataset from Baccin et al.

    • The fact that bone marrow adipose tissue begins to accumulate after birth is well known. It is therefore not surprising that adipogenic progenitor populations start to accumulate perinatally (established by studies cited by the authors). Thus, these results only confirm the validity of the dataset. This represents an example on how the majority of findings have been presented here.*

    We fully agree that some of our results confirm, at the single cell level, knowledge previously stablished with other methods. However, and continuing with the case of adipose tissue mentioned by the reviewer, the analysis of our datasets with unbiased tools allowed the identification of fibroblastic populations, such as AFP or GFP, which shown by GO terms and, most importantly, by highly-relevant regulons identified by SCENIC, to be potentially associated with thermogenesis and brown fat differentiation. As far as we know, the specific transcriptional regulators involved in brown fat differentiation in the bone are poorly defined. In addition, adipogenesis is not the only aspect we highlight, being other novel association the putative interaction between fibroblastic mesenchymal populations and Eo/Bas and ILC-ISLP2 hematopoietic cells. These are just two examples of relevant aspects uncovered by that our holistic analysis of all bone population, and that might be further explored by specialized groups in the respective fields.

    • Given that the authors do not provide functional validation of putative new molecular interactions (by CellPhoneDB) their conclusions should be presented in a more tempered manner and acknowledged as inference rather than fact.*

    We agree with the reviewer and accordingly, we have tuned-down several statements throughout the text (see also response to Reviewer #1).

    • Similarly, the authors claim "...we identified Ptx4 as a novel tenogenic-specific gene...". This is too strong a conclusion as this has not been functionally validated. It* should at least be tested by immuno-(co)-staining.

    Being Ptx4 a secreted molecule, it would be very difficult to reliably assign the signal to a specific population by co-immunolocalization with bona fide tenogenic markers such as Scx or Tnmd. Besides, when pointing out Ptx4 specific expression in the tenogenic branch of the TC cluster, we intended to suggest the potential use of this locus for the generation of novel genetic tools. We have reformulated this sentence to clearly indicate this and avoid claiming that Ptx4 is a novel tenogenic marker.

    • The authors identify "uncommitted clusters" as mesenchymal progenitor populations without actual showing that they are even related by lineage. This is a general pitfall in analyzing single cell RNA-sequencing data and making trajectory/pseudotime inferences. It is now well-established that the mesenchymal compartment is highly heterogeneous and composed of multiple distinct cellular lineages. Trajectory inference tools such as PHATE do not distinguish those different mesenchymal lineages. As such, the presented results cannot be considered valid unless there is proper functional validation.*

    We agree with the reviewer on the limitation of this type of analysis, such as its inability to resolve phenotypic convergence (e.g. the case for osteoblasts generated from reprogrammed hypertrophic chondrocytes or from perichondrial cells). We have therefore removed the PHATE data from the manuscript.

    • The description of PaS being mainly associated with compact bone is neither correct nor supported by cited studies. The authors show potential additional markers to target Pas in mice, but fail to validate their point that these markers could be used in human tissue as well.*

    We thank the reviewer for pointing this out and we apologize for our incorrect wording. What we intended to mean is that PaS cells can only be efficiently extracted by enzymatic treatment of the bone fraction after bone marrow aspiration in adults. We have now corrected these instances in the revised manuscript.

    Concerning the validation in human samples of the proposed additional markers for the PaS population, we agree that this is an important point, but one that would require the processing of fetal/newborn human bone tissue for FC, which is beyond our capacities and the scope of the current manuscript

    Figure 1c: the legend for dotsize is off scale.

    We thank the reviewer for spotting this mistake, which inadvertedly happened during figure assembly and is now corrected.

    Reviewer #2 (Significance (Required)):

    Strength:

    • thorough analysis of single cell RNA-sequencing datasets including integration of published work
    • good writing and figure presentation
    • dataset fills gap for the field as the presented ages have not been published

    Limitations:

    • lacking functional validation
    • lack of new biology - mostly confirmation of known facts

    Advance:

    • knowledge gain is incremental
    • good resource

    Audience:

    • fills a gap in the bone and hematopoietic research field as a resource

    My expertise: Skeletal stem cell lineage biology, single cell RNA-sequencing of bone cell populations.*

    __*Reviewer #3 (Evidence, reproducibility and clarity (Required)):

    *__

    Summary The authors present data on a very interesting model, mouse bone just before and after birth. In this timeframe, the organism has to adapt from a buoyant, nurtured environment to stronger gravitational forces acting upon the skeletal structure, changed oxygen uptake, changed demands to the immune system and its development, and an overall changed metabolism. The authors introduce these changes and their importance in a clear, easy-to read introduction, and this clear structure and language continue throughout the manuscript. Comparing scRNA-seq of bone E18.5 and adult stage, comparable findings by Liu Y. et al., (https://doi.org/10.1038/s41467-022-28775-x*) have been previously shown. However, this manuscript showed additional postnatal day 1 (P1) data. *All computational analyses are well done, for the most part well described, and, notably, the integration of previously published data allows us to put the results of this study into context and compare them to the adult situation. Data sharing is not optimal, but it is already very good. The only downside is that most of the computational analyses are done at a very limited level of depth and merely provide initial insights and an overview of the data presented.

    We thank the reviewer for the overall positive view of our manuscript and his/her critical comments which we have tried to address in our revisions. We wish to apologize for our omission in citing the Liu et al. study, which is now corrected in the revised text. In this respect, we would like to point out that the Liu et al. study is mostly centered in the endothelial compartment, whereas our work is more focused on the mesenchymal populations. Hence, both studies are complementary. Of note, Liu et al. were able to detect Wnt2 expression in E18.5 endothelial cells using targeted single-cell RNA-seq for a panel of specific genes, while in our data, more focused on the mesenchymal compartment, Wnt2 expression maps mostly to the SFP fibroblastic cluster, with low expression in few endothelial cells. In our view, this apparent discrepancy is not such, but the result of different strategies of sorting and enrichment, and illustrates the need of having complementary studies and datasets (e.g the SFP populations may also be an additional source of Wnt2 to promote hematopoietic stem and progenitor cell proliferation, as reported by Liu et al.).

    As for the limited depth on our analysis mentioned by the reviewer, we would like to point out that we made a major effort to put our observations in the context of the vast literature on both the mesenchymal and hematopoietic compartments, which forced us to synthetize in the main text. When possible, we added additional data as part of the supplementary information (e.g. full CellPhoneDB inferred interactions as Excel tables).

    *Further comments will be given in bullet-point form, split by their impact on the overall message of the manuscript. *

    Major • The authors provide FACS results for the cell clusters and types they defined from the scRNA-Seq data but do not provide any results on where and in which cellular contexts these cells are found in the bone and whether their spatial proximity and proteome (via staining, for example) make it likely or unlikely to see the cell-cell communication suggested in their CellPhoneDB analysis. The authors should either provide such results or adjust sentences as follows to not overstate their results: "These analyses also unveiled the complex MC-HC connectome, in particular the abundant interactions of fibroblastic SFP, AFP, CLFP, and GFP populations with HPC, and quite outstandingly, with the ILC-TSP2 and Eo/Bas clusters."*

    We agree with the reviewer but, as previously commented in the General Statements and in our response to Reviewer #1, spatial localization of all populations is technically challenging in the bone and the methods we have tested fall short for the precise and reliable localization of specific bone cell populations with cellular resolution. Following his/her suggestion, we have systematically edited the text so as to not overstate any message stemming from our expression analysis.

    • To address the issue of lineage commitment, the authors could offer some functional assessments between E18.5 and P1 or Adult bone BMSCs stromal cells subset that were sorted using FACS (Fig. 6).*

    As commented in our response to Reviewer #2, and given the complex lineage relations in the bone, addressing this point would require extensive in vivo experimentation through transplant surgery in immunocompromised newborns or genetic analysis using novel mouse alleles, both of which we consider out of the scope of our study, conceived as a resource. Of note, in Fig. 6 we did not analyze by FC any adult population, but in the revised Fig. 6A, we now provide the expression of all markers in perinatal and adult datasets. In relation with lineage relations, we have also removed the PHATE analysis from the revised manuscript, as suggested by Reviewer #2.

    Minor • The visualization of UMAP embeddings is very inconsistent across the manuscript and misleading or irritating in some cases. For example, in Figure 1b, the separation of the background grid is not clearly visible between E18 and PN1. In other figures in the same manuscript, borders around the figures solve this issue. Additionally, axes are either missing or unlabeled, whereas for UMAP embeddings, irrelevant axis tick labels and grid lines are present in most figures. It would benefit the overall flow and visualization of the manuscript if UMAP figures were more consistent.

    We apologize for these inconsistencies and we thank the reviewer for pointing them out. We have now separated both panels in Fig. 1B and added borders so as to separate both UMAP plots. We have also added the missing labeling of axes throughout the manuscript so as to make all figures more consistent. We have chosen to keep both grid lines and tick labels as they help in the comparison of Harmony-integrated datasets.

    • For Figure 1b specifically, it might also make sense to outline the main cell populations in both UMAPs, as in Figure 2a.*

    Agreed and done.

    • Average gene expression cannot take on values below 0, as that is the lower bound for expression counts. Figure 1c seems to show the colorbar dropping below 0 though. This might just be a problem of confusing color bar label placement, but it should be addressed. It should also be assured that, indeed, there has not been a mix-up and expression values are limited to >=0.*

    We should have explained this better. These dotplots in Fig.1 and the heatmap in Fig. S1 use the normalized and scaled expression value (mean=0; standard deviation=1), which means that it might be negative expression values. These instances are interpreted as genes in which the expression levels are lower than the mean expression level in the dataset and facilitate the visualization of differential gene expression in the different clusters. We have now indicated this clearly in the figure legends.

    • For the PHATE analysis, was there any batch correction applied to address potential batch effects between the E18 and PN1 datasets?*

    We have removed from the manuscript all PHATE analysis. Still, as we use Harmony integration as a batch-correction tool, we now describe it now in detail in the Materials and Methods section.

    • Figure 4a: From the text, it is clear that CellPhoneDB was used to calculate significant interactions between cell types. However, it is not clear which threshold (even if default) was used to determine what constitutes a significant interaction.*

    We apologize for this omission. We have indicated both in the revised figure legend and in Materials and Methods the threshold (p-value ≤0.05; as calculated by CellPhoneDB) that was used to represent all significant interactions and shown in Fig. 4A.

    • Figure 4b: It is unclear why a collection of chord representations was chosen here, as chord diagrams of this kind generally do not provide any useful additional information apart from an interaction being found to be significant (by a certain threshold) between two cell types. Lacking are generally more interesting parameters, such as the interaction score of such interactions or the expression of the involved ligands and receptors, in comparison to other cell types, where the respective interaction was not predicted to be significant. In this particular case, it is also unclear what is encoded by the width of the respective arrows. This should be made clear. Additionally, a suggestion could be to either present this information in an array of two DotPlots, one for ligands and receptors, respectively, or to encode additional information in, for example, the arrow or connector width, with the connector encoding the mean ligand expression and the arrow head encoding the mean receptor expression in the chord diagram.*

    We initially chose to use chord plots as we thought it would be a visual way to represent significant interactions but, as the reviewer points out, they do not provide any additional information. In the line with the reviewer’s suggestion, we have substituted all Fig. 4B chord representations for bubble plots in which are both encoded the mean scaled expression of the ligand/receptor pair (the output of the CellPhoneDB tool) and the mean percentage of cells in the clusters expressing the corresponding molecules. We believe that this modification makes this figure more informative and visually easier to interpret.

    • The authors do not mention how many genes were used as marker genes for GFP, SFP, etc. for the GO term enrichment analysis. This number (if low), the significance cut-offs, and the method used to determine DEGs could potentially have an impact on the GO enrichment results. The authors should therefore, already in the main manuscript text, mention the number of genes used for each of these cell subtypes and the method used to determine them. The text mentions cellranger, but the underlying methodology is not mentioned.*

    In the revised manuscript, we have included how differential gene expression between clusters was calculated (DEGs were obtained using the FindAllMarkers() function in Seurat, using the default parameters -by default Seurat uses the Wilcoxon Rank Sum test for statistical testing) and the genes used for GO analysis (DEGs were filtered to include genes with an adjusted p-value ≤0.005; gene lists provided as new Supplementary Table 1). The resulting number of genes used for GO analysis at E18.5/PN1 was 218/280 (AFP), 455/480 (CLFP), 185/234 (GFP) and 436/305 (SFP). Retrieved GO terms were filtered by a ratio fold of enriched/expected ≥ 2 and manually curated.

    Reviewer #3 (Significance (Required)):

    This study and single cell RNA-sequencing data further analyze the distinctions between the neonatal and adult stages of hematopoietic cells and bone stromal cells. This study also demonstrated the cellular heterogeneity of hematopoietic and bone stromal cells, as well as how cellular cross-talk supports osteogenic and hematopoietic cells. This sequencing data will be useful in the future to comprehend how the bone and marrow adapt to a stronger gravitational force operating on the skeletal structure, as well as to changed oxygen consumption, requirements for the development of the immune system, and an overall altered metabolism.*

  2. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #3

    Evidence, reproducibility and clarity

    Summary

    The authors present data on a very interesting model, mouse bone just before and after birth. In this timeframe, the organism has to adapt from a buoyant, nurtured environment to stronger gravitational forces acting upon the skeletal structure, changed oxygen uptake, changed demands to the immune system and its development, and an overall changed metabolism. The authors introduce these changes and their importance in a clear, easy-to read introduction, and this clear structure and language continue throughout the manuscript. Comparing scRNA-seq of bone E18.5 and adult stage, comparable findings by Liu Y. et al., (https://doi.org/10.1038/s41467-022-28775-x) have been previously shown. However, this manuscript showed additional postnatal day 1 (P1) data.

    All computational analyses are well done, for the most part well described, and, notably, the integration of previously published data allows us to put the results of this study into context and compare them to the adult situation. Data sharing is not optimal, but it is already very good. The only downside is that most of the computational analyses are done at a very limited level of depth and merely provide initial insights and an overview of the data presented. Further comments will be given in bullet-point form, split by their impact on the overall message of the manuscript.

    Major

    • The authors provide FACS results for the cell clusters and types they defined from the scRNA-Seq data but do not provide any results on where and in which cellular contexts these cells are found in the bone and whether their spatial proximity and proteome (via staining, for example) make it likely or unlikely to see the cell-cell communication suggested in their CellPhoneDB analysis. The authors should either provide such results or adjust sentences as follows to not overstate their results:
      • "These analyses also unveiled the complex MC-HC12 connectome, in particular the abundant interactions of fibroblastic SFP, AFP, CLFP13, and GFP populations with HPC, and quite outstandingly, with the ILC-TSP2 and 14 Eo/Bas clusters."
    • To address the issue of lineage commitment, the authors could offer some functional assessments between E18.5 and P1 or Adult bone BMSCs stromal cells subset that were sorted using FACS (Fig. 6).

    Minor

    • The visualization of UMAP embeddings is very inconsistent across the manuscript and misleading or irritating in some cases. For example, in Figure 1b, the separation of the background grid is not clearly visible between E18 and PN1. In other figures in the same manuscript, borders around the figures solve this issue. Additionally, axes are either missing or unlabeled, whereas for UMAP embeddings, irrelevant axis tick labels and grid lines are present in most figures. It would benefit the overall flow and visualization of the manuscript if UMAP figures were more consistent.
    • For Figure 1b specifically, it might also make sense to outline the main cell populations in both UMAPs, as in Figure 2a.
    • Average gene expression cannot take on values below 0, as that is the lower bound for expression counts. Figure 1c seems to show the colorbar dropping below 0 though. This might just be a problem of confusing color bar label placement, but it should be addressed. It should also be assured that, indeed, there has not been a mix-up and expression values are limited to >=0.
    • For the PHATE analysis, was there any batch correction applied to address potential batch effects between the E18 and PN1 datasets?
    • Figure 4a: From the text, it is clear that CellPhoneDB was used to calculate significant interactions between cell types. However, it is not clear which threshold (even if default) was used to determine what constitutes a significant interaction.
    • Figure 4b: It is unclear why a collection of chord representations was chosen here, as chord diagrams of this kind generally do not provide any useful additional information apart from an interaction being found to be significant (by a certain threshold) between two cell types. Lacking are generally more interesting parameters, such as the interaction score of such interactions or the expression of the involved ligands and receptors, in comparison to other cell types, where the respective interaction was not predicted to be significant. In this particular case, it is also unclear what is encoded by the width of the respective arrows. This should be made clear. Additionally, a suggestion could be to either present this information in an array of two DotPlots, one for ligands and receptors, respectively, or to encode additional information in, for example, the arrow or connector width, with the connector encoding the mean ligand expression and the arrow head encoding the mean receptor expression in the chord diagram.
    • The authors do not mention how many genes were used as marker genes for GFP, SFP, etc. for the GO term enrichment analysis. This number (if low), the significance cut-offs, and the method used to determine DEGs could potentially have an impact on the GO enrichment results. The authors should therefore, already in the main manuscript text, mention the number of genes used for each of these cell subtypes and the method used to determine them. The text mentions cellranger, but the underlying methodology is not mentioned.

    Significance

    This study and single cell RNA-sequencing data further analyze the distinctions between the neonatal and adult stages of hematopoietic cells and bone stromal cells. This study also demonstrated the cellular heterogeneity of hematopoietic and bone stromal cells, as well as how cellular cross-talk supports osteogenic and hematopoietic cells. This sequencing data will be useful in the future to comprehend how the bone and marrow adapt to a stronger gravitational force operating on the skeletal structure, as well as to changed oxygen consumption, requirements for the development of the immune system, and an overall altered metabolism.

  3. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    In the study by Rueda et al. the authors use single cell RNA-sequencing to investigate differences in cellular composition bones/bone marrow between late gestational stage mouse embryos and their perinatal counterparts. The authors describe specific differences in the relative abundance of putative cell types and use established bioinformatic tools to infer interactions as well as molecular mechanisms determining specific functions. The employed methods are well described and the results are presented in a very clear and understandable manner. Despite that, the findings do not provide any substantial knowledge advance but are rather confirming work of published literature while supplementing available single cell RNA-sequencing datasets of mouse bones at adult ages. As such, this work provides an interesting resource but does not report novel biology.

    Major comment: The authors explore the interesting transition of embryonic to perinatal bone/bone marrow using single RNA-sequencing. This fills a gap for the field of bone and hematopoietic researchers. There is little to criticize about the presented data. However, while it provides a nice resource, the knowledge gained is incremental. As acknowledged by the authors themselves, their study lacks functional validation of any findings made or conclusions drawn from bioinformatic tools in this manuscript. They use published work to validate their findings but do not go beyond that to confirm putative new biology. Some examples are listed in the minor comments.

    Minor comments:

    1. Remark: 10X Chromium does not provide whole transcriptomic coverage but rather captures the most highly abundant transcripts without for example being able to distinguish alternatively spliced gene variants. Based on that, interpretation of gene expression, or the absence of a gene in the dataset, should be interpreted carefully.
    2. The fact that bone marrow adipose tissue begins to accumulate after birth is well known. It is therefore not surprising that adipogenic progenitor populations start to accumulate perinatally (established by studies cited by the authors). Thus, these results only confirm the validity of the dataset. This represents an example on how the majority of findings have been presented here.
    3. Given that the authors do not provide functional validation of putative new molecular interactions (by CellPhoneDB) their conclusions should be presented in a more tempered manner and acknowledged as inference rather than fact.
    4. Similarly, the authors claim "...we identified Ptx4 as a novel tenogenic-specific gene...". This is too strong a conclusion as this has not been functionally validated. It should at least be tested by immuno-(co)-staining.
    5. The authors identify "uncommitted clusters" as mesenchymal progenitor populations without actual showing that they are even related by lineage. This is a general pitfall in analyzing single cell RNA-sequencing data and making trajectory/pseudotime inferences. It is now well-established that the mesenchymal compartment is highly heterogeneous and composed of multiple distinct cellular lineages. Trajectory inference tools such as PHATE do not distinguish those different mesenchymal lineages. As such, the presented results cannot be considered valid unless there is proper functional validation.
    6. The description of Pas being mainly associated with compact bone is neither correct nor supported by cited studies. The authors show potential additional markers to target Pas in mice, but fail to validate their point that these markers could be used in human tissue as well.
    7. Figure 1c: the legend for dotsize is off scale.

    Significance

    Strength:

    • thorough analysis of single cell RNA-sequencing datasets including integration of published work
    • good writing and figure presentation
    • dataset fills gap for the field as the presented ages have not been published

    Limitations:

    • lacking functional validation
    • lack of new biology - mostly confirmation of known facts

    Advance:

    • knowledge gain is incremental
    • good resource

    Audience:

    • fills a gap in the bone and hematopoietic research field as a resource

    My expertise: Skeletal stem cell lineage biology, single cell RNA-sequencing of bone cell populations

  4. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    In this new manuscript, Rueda and colleagues present an extensive bioinformatics analysis of single cell transcriptomic data obtained for mouse endochondral bone cell populations before and after birth. They describe gene markers of mesenchymal and hematopoietic cells pointing to differences with adult bone populations, and they use gene ontology and trajectory analyses to infer possible roles of these cells in the developing bone. The data could provide a valuable resource for further understanding endochondral bone development and the changes driving this process in peri-natal stages. However, they are also significant weaknesses.

    A major weakness is that the scRNA-seq data lack validation through other techniques and functional assays. Namely, in situ data are missing to locate the various cell populations in the developing bones, especially the different types of fibroblastic cells identified by the authors. Such data would go a long way to understand the possible functions of the cell populations. Although the authors tried to complement their data with a review of the literature, most of the conclusions remain purely speculative and not sufficiently supported by scientific and statistic rigor. This makes the Results section more like a discussion than a description of the results. For instance, the authors proposed important regulatory functions for the fibroblastic clusters, but there is no data supporting this other than broad GO terms associated with genes expressed in these cells. Related to this point, the title of the manuscript does not accurately reflect the content of the study.

    Other points:

    1. The authors missed to report in the Results section which skeletal elements they used for their analyses and which skeletal elements were used for the adult dataset that they compared their data with. Differences in skeletal elements and in the ways whereby these samples were collected and processed could explain differences detected in the two types of datasets. Also, the sex and age of the samples for the adult dataset should be reported.
    2. It is unclear whether PN1 is the day that mice are born (classically referred to as P0) or the next day.
    3. It is unclear whether the cells obtained for each biological replicate were pooled for the scRNA-seq assays or were treated individually. It is thus unclear how reproducible the data are.
    4. It is not clear in the gating strategy chosen for the flow cytometry as shown in Fig. 1A why the green gate containing cells expressing high levels of CD9, CD140 and CD31 has been extended in between the purple and orange gates containing CD140 and CD31 negative cells.
    5. Are cells from all the sequenced samples homogenously distributed in the scRNA-seq clusters? Authors should provide this information and add statistic when they describe changes in the amount of cells per cluster between E18.5 and PN1 stages.
    6. On the basis of what markers the AFP population has been called adipogenic? Authors present Ptch2 and Notch3 as markers of this cluster, but not adipogenic progenitor genes.
    7. Authors claim that there is a good correlation between OsC and osteo-CAR clusters. However, OsC cells do not express Cxcl12 and other typical CAR cell markers.
    8. In Figure 6 expression of PaS cell markers should be shown for both adult and perinatal populations. Additionally, have the authors tested that the sorted cells in panel C have the same progenitor properties as the PaS cells?

    Significance

    As indicated in the comments for the authors, the new scRNA-seq data could become a useful resource for subsequent studies, but they are at present insufficient to represent a significant scientific advancement. The main concern is that new cell populations appear to have been identified by the authors, but a number of questions were not answered such as regarding their actual location in the skeletal elements, their origins, their fates and their functions. Generating such data would require a major amount of effort and require substantial revision of the manuscript.

  5. Version published to 10.1101/2023.07.18.549529v3 on bioRxiv
  6. Version published to 10.1101/2023.07.18.549529v2 on bioRxiv
  7. Version published to 10.1101/2023.07.18.549529v1 on bioRxiv