Single-cell RNA sequencing reveals cellular and molecular heterogeneity in fibrocartilaginous enthesis formation

  1. Tao Zhang
  2. Liyang Wan
  3. Han Xiao
  4. Linfeng Wang
  5. Jianzhong Hu
  6. Hongbin Lu  Is a corresponding author
  1. Department of Sports Medicine, Xiangya Hospital Central South University, China
  2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, China
  3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, China
  4. Department of Spine Surgery and Orthopaedics, Xiangya Hospital, Central South University, China

Abstract

The attachment site of the rotator cuff (RC) is a classic fibrocartilaginous enthesis, which is the junction between bone and tendon with typical characteristics of a fibrocartilage transition zone. Enthesis development has historically been studied with lineage tracing of individual genes selected a priori, which does not allow for the determination of single-cell landscapes yielding mature cell types and tissues. Here, in together with open-source GSE182997 datasets (three samples) provided by Fang et al., we applied Single-cell RNA sequencing (scRNA-seq) to delineate the comprehensive postnatal RC enthesis growth and the temporal atlas from as early as postnatal day 1 up to postnatal week 8. And, we furtherly performed single-cell spatial transcriptomic sequencing on postnatal day 1 mouse enthesis, in order to deconvolute bone-tendon junction (BTJ) chondrocytes onto spatial spots. In summary, we deciphered the cellular heterogeneity and the molecular dynamics during fibrocartilage differentiation. Combined with current spatial transcriptomic data, our results provide a transcriptional resource that will support future investigations of enthesis development at the mechanistic level and may shed light on the strategies for enhanced RC healing outcomes.

Editor's evaluation

This paper represents a valuable single-cell level analysis of tendon enthesis development. The study allows further understanding of this specific process with clinical implications. The authors provided convincing evidence for the heterogeneity of postnatal enthesis growth and the molecular dynamics and signaling networks during enthesis formation.

https://doi.org/10.7554/eLife.85873.sa0

Introduction

RC and its enthesis are essential components of the shoulder, which are critical in facilitating coordinated shoulder movements and stability (Schett et al., 2017). Compared to other RC tissues, the attachment site of supraspinatus (SS) tendon is vulnerable to injury and difficult to achieve complete regeneration, due to its high heterogeneity in composition and structure (Nourissat et al., 2015; Schett et al., 2017). Histologically, the attachment site of the SS tendon is a classic fibrocartilaginous enthesis, also termed BTJ (Chen et al., 2021a; Rossetti et al., 2017). In its native state, the fibrocartilaginous enthesis exhibits gradations in tissue organization, cell phenotype, and matrix composition (Moffat et al., 2008; Rossetti et al., 2017). Fibrocartilaginous entheses manage to disperse stress and facilitate load transfer between vastly different materials like tendons and bone, with modulus ranging from 200 MPa to 20 GPa (Rossetti et al., 2017). Unfortunately, the intrinsic regenerative capacity of fibrocartilaginous enthesis is not well-understood, which limits the exploitation of the best and most rigorously proven early intervention programs (Derwin et al., 2018; Schett et al., 2017; Xiao et al., 2022). Therefore, understanding the complex process of the enthesis morphogenesis and maturation during development may inform strategies for enhanced BTJ healing.

Currently, the mechanism underlying the growthof the enthesis fibrocartilage is less understood. Details are scarce, but the fibrocartilage layer is formed by a pool of site-specific progenitor cells, and initially organizes as an unmineralized cartilaginous attachment unit (Jensen et al., 2018). Such development pattern shares an overlapping biological behavior with the growth plate, which is a process of mesenchymal stem cells differentiating into chondrogenic cells and then sequentially into fibrocartilage cells (Killian, 2022). A unique enthesis progenitor pool has been identified since the embryonic stages, as cells sandwiched between primary cartilage and tendon, expressing a mixed transcriptome of both chondrogenic and tenogenic genomic features (Scleraxis and SRY-related transcription factor) (Blitz et al., 2013; Sugimoto et al., 2013). These enthesis progenitors can differentiate into either chondrocytes or tendon fibroblasts, under the regulation of Krüppel-like (Klf) transcription factors (Kult et al., 2021). In the later embryonic and postnatal stages, cells from the enthesis progenitor pool ultimately either differentiate into or are replaced by a Hh-positive cell population marked by Gli1 and Ptch1 (Felsenthal et al., 2018; Schwartz et al., 2015). Gli1+ cells and their progenies are retained in the enthesis region throughout postnatal development and eventually populate the entire fibrocartilage region between tendon and bone, thereby contributing to enthesis growth (Felsenthal et al., 2018; Jensen et al., 2018; Schwartz et al., 2015). However, compared with our in-depth understanding of enthesis development during the embryonic stage, the cell-type composition and distribution in the enthesis at different postnatal stages, as well as biochemical markers for enthesis stem cell progenitor used in tissue engineering, remain to be not well-understood and require novel methods for elucidation.

As is proved, the normal enthesis maintains a gradient of cell phenotypes, from tendon fibroblast to chondrocyte then to mineralizing chondrocyte and to osteoblast/osteocyte (Chen et al., 2021b; Moffat et al., 2008; Rossetti et al., 2017). It is unclear how this gradient in cell phenotypes develops and how it is regulated by the local environment (e.g. extracellular matrix, muscle loading, and growth factors) (Derwin et al., 2018). scRNA-seq is a powerful method to analyze various cell types and provide insights into tendon enthesis postnatal development (Gulati et al., 2020; Kult et al., 2021). Recently, Fang Fei et al., reported an exciting single-cell work to define the enthesis cell transcriptomes at postnatal day 11, 18, and 56, revealed the clonogenicity and multipotency of enthesis Gli1+ progenitors (Fang et al., 2022). However, compared to the abundant and growing amount of single-cell resources in bone, cartilage, and tendon, more enthesis single-cell resources are needed to cover long-range timepoints of BTJ development or entheseal diseases. Here, we applied single-cell transcriptomics along with the spatial transcriptomic sequencing to analyze the cellular and molecular dynamics during postnatal tendon enthesis growth. The results provided here for deciphering postnatal tendon enthesis development may facilitate future studies of enthesis regeneration.

Results

The development of fibrocartilage in the enthesis occurs postnatally

The heterogeneity of the fibrochondrocytes in enthesis development has been an open question. Aiming to evaluate the postnatal development pattern of enthesis fibrocartilage, we first stained the shoulder sections and compared the morphological parameters of the enthesis cells at P1, P3, P7, P14, P28, and P56 (Figure 1a). At postnatal day 1, an articular cavity was formed and the supraspinatus tendon (ST) was observed attached to the humeral head. The cells at the ST attachment site were highly dense and homogeneous, and were visibly different from tendon cells and the primary cartilage cells inside the humeral head, suggesting that postnatal enthesis is formed by site-specific progenitor cells since the embryonic stage. At P1, P3, and P7, the fibrocartilage layer and subchondral bone could hardly be visualized. Fibrocartilage was not evident at the enthesis of the mouse rotator cuff until 2–3 weeks after birth. To examine the cellular morphological changes, we further measured the 2D parameters (including 2D area, roundness, frete diameter, and minimal frete diameter) of enthesis cells, statistics results demonstrated that cell size remarkably increased during postnatal development, from postnatal day 7–14 (Figure 1d). After postnatal day 14, the enthesis cells were observed as typical chondrocyte phenotype, column-like stacked alongside the direction of tendon fiber, with more prominent patterns at P28 and P56. From toluidine blue/fast green stainings, fibrochondrocytes can be observed since P14, with a larger cellular size in comparison with condensed primary chondrocytes in P7. The extracellular matrix (ECM) of fibrocartilage with mineralization can be observed after P28. Type 2 collagen is a major ECM component of cartilage, thus we stained enthesis with col2a1 antibody. The IHC results show the col2a1 protein levels are relatively low at P14, and significantly increase after P28 (Figure 1b and c).

The development of fibrocartilage in the enthesis occurs postnatally.

(a) H&E staining of P1, P3, P7, P14, P28, and P58 mouse supraspinatus tendon entheses (n=4). Scale bars, 100 μm. (b) Immunohistochemistry images of cartilage-abundant collagen II at P14, P28 and P56. Scale bars, 100 μm. (c) Comparison of the col2a1 positive area derived from IHC results. Error bars represent SEM. N=3-4. **p<0.01, ****p<0.0001. (d) Comparison of the morphological parameters (2D area, roundness, frete diameter, and minimal frete diameter) between P1, P3, P7, P14, P28, and P56. Error bars represent SEM. *p<0.05, ***p<0.001, ****p < 0.0001.

Unbiased clustering identified known cell populations in postnatal enthesis development

To determine the cellular composition of the developing enthesis, we integrated our dataset with open-source GSE182997 datasets (three samples) provided by Fang et al., (Figure 2a). After the elimination of doublets, dead, and apoptotic cells, blood cells (erythrocytes and progenitors), endothelial cells, immune cells (B cells and T cells), myeloid cells, and growth plate chondrocytes, we got high-quality transcriptomic data from 8368 single cells, including 1285 P1 cells, 4059 P7 cells, 918 P11 cells, 307 P14 cells, 329 P18 cells, 568 P28 cells, and 897 P28 cells (Figure 2—figure supplement 1). Unbiased clustering based on Uniform Manifold Approximation and Projection (UMAP) identified major cell populations. Based on the differentially expressed genes (DEGs) (Figure 2—figure supplement 2), all the cell clusters were annotated, including BTJ chondrocytes, BTJ tendons, Tenocytes, Osteocytes, Enthesoblasts, Tenoblasts, and Mesenchymal progenitors (Figure 2b). We used an entropy-based metric (ROUGE) for assessing the purity of single-cell clusters, all these cell subtypes achieved high ROGUE values of >0.9 (Liu et al., 2020), which suggested accuracy in unsupervised cell clustering (Figure 2—figure supplement 2). Generally, enthesis progenitors had relatively high expression of Ly6a, Cd44, and Pdgfrα, in agreement with progenitors found in tendons and bone marrow (Harvey et al., 2019; Tikhonova et al., 2019). Enthesoblasts were defined based on their relatively decreased stemness transcriptional signatures, as well as co-expressed tenogenic (e.g. Scx, Col1a1) and chondrogenic markers (e.g. Sox9, Acan) (Jensen et al., 2018). The correlation analysis showed the consistency between our datasets and Fang et al., reported datasets (Figure 2—figure supplement 2). The barplot diagram showed cell identity change after birth (Figure 2c).

Figure 2 with 2 supplements see all
Unbiased clustering identified Known Cell Populations in postnatal enthesis development.

(a) Schematic workflow of the study design. (b) Distributions of seven cell clusters on UMAP plot, including bone-tendon junction (BTJ) chondrocytes, BTJ tendons, tenocytes, osteocytes, enthesoblasts, tenoblasts, and mesenchymal progenitors. (c) Fractions of cell clusters in enthesis development at P1, P7, P11, P14, P18, P28, and P56. (d) The average expression of curated feature genes for previously reported enthesis marker genes and enthesis-specific extracellular matrix (ECM) genes. (e) Spatial transcriptomic spot map revealing the expression of chondrocyte marker genes in each spatial spot. (f) Representative immunofluorescence staining to validate the spatial distribution of Sox9+ and Scx+ cells in the enthesis area, at P1, P7, and P14. Scale bars, 100 μm. N=3-5. *p<0.05, ***p<0.001, ****p<0.0001.

We used the spatial transcriptome data of mice 1 day after birth to verify the cell population defined in our P1 single-cell dataset, and the enthesis chondrocytes group was consistent with the spatial anatomical position of the bone-tendon junction (Figure 2e). We checked the previously reported enthesis marker genes Sox9 and Scx, as well as enthesis-specific ECM genes (Col2a1), which were ubiquitously expressed in BTJ chondrocyte cells (Figure 2e). We then performed an immunofluorescence assay to validate the spatial distribution of enthesis-related genes, we found that Sox9+ and Scx+ cells were detected in the enthesis area, mostly in the neonatal stage and significantly decreased in postnatal weeks 2 and 4, as expected (Figure 2f).

Identifying developmental trajectories for tenocytes, chondrocytes, and osteocytes differentiation in enthesis

We next sought to investigate the trajectory and regulatory genes to govern postnatal bone-tendon junction cell development. We first predicted the differentiation state of each cell group from scRNA-seq data by using Cytotrace (Gulati et al., 2020). The Cytotrace results showed that across all cell clusters, the ‘stemness’ degree of mesenchymal progenitors is higher than other cell types (Figure 3a). Among tendon and enthesis-associated cell groups, CytoTRACE scores of the enthesis cells (BTJ chondrocytes cells) skewed toward a moderate predicted stem potential, which was slightly higher than that in tendon cells, suggestive of a higher degree of stemness for enthesis cells developing into fibrocartilage cells (Figure 3b).

Identifying developmental trajectories for tenocytes, chondrocytes, and osteocytes differentiation in enthesis.

(a) Uniform Manifold Approximation and Projection (UMAP) plot of enthesis single-cell RNA sequencing (scRNA-seq) data overlaid with CytoTRACE scores. (b) Boxplot of predicted differentiation score distributions for each cell cluster. (c) Results of RNA velocity analysis show that postnatal enthesis fibrocartilage origin from mesenchymal progenitors, instead of tendon cells. (d) Cellrank identified three differentiation terminal cell types, including bone-tendon junction (BTJ) chondrocytes, tenocytes, and osteocytes. (e) Fate probabilities uncovered putative BTJ chondrocytes lineage drivers. (f) Distribution over cluster membership for each of the cells assigned to a certain terminal state. (g) Genes that correlate positively with the BTJ chondrocyte fate correlate moderately with the tendon fate and vice versa. (h) Representative tenogenic, chondrogenic, and osteogenic gene expression dynamics along four terminal differentiation trajectories.

We next implemented RNA velocity analysis. The trajectory results showed enthesoblasts originated from mesenchymal progenitors, and enthesoblasts had a greater chance to branch into enthesis chondrocytes (Figure 3c). Next, we employed developmental potential calculated by Cytotrace to compute the terminal states among all the cell groups by Cellrank, which is a toolkit based on Markov state modeling (Lange et al., 2022), and four clusters of cells were predicted as the terminal differentiate clusters: BTJ chondrocytes, tenocytes, and two subsets of osteocytes (Figure 3d), suggesting that the differentiation profile of these cell groups was separative from each other. The fate probabilities analysis showed that only BTJ chondrocytes contributed to the fate of enthesis chondrocytes (Figure 3e and f), which was consistent with RNA velocity, suggesting that the fibrocartilage in postnatal enthesis origin from enthesis site-specific progenitors, instead of tendon cells. However, we found BTJ chondrocytes correlate slightly with tenocytes, because the genes that correlate positively with the BTJ chondrocyte fate correlate moderately with the tendon fate and vice versa (Figure 3g).

To explore gene expression dynamics along the trajectories, we measured the dynamics of genes in pseudotime along the differentiation trajectories of these three differentiation terminal cell types. We found that the expressions of known differential regulator genes were upregulated significantly higher in enthesis-associated trajectories, such as Scx, Sox9, and Tnmd, which were confirmed indispensable for enthesis formation. We also found genes related to cartilage ECM (Acan, Scrg1) were upregulated along the pseudotime, and highly expressed until the terminal differentiation state, suggesting the collagen and matrix protein synthesis were predominant in postnatal enthesis growth. Meanwhile, mineralization-related genes (Ibsp, Col11a1) were observed to increase in enthesis chondrocyte differentiation, but less than their expression levels in osteocytes (Figure 3h).

Reconstruction of the trajectory and gene dynamics in BTJ chondrocytes differentiation

As enthesis progenitors differentiating into fibrochondrocytes is the pivotal step of enthesis development, we sought to determine their gene dynamics. Mesenchymal progenitors, enthesoblasts, and enthesis chondrocytes were subsets to receive trajectory analysis. We used monocle3 to identify the unsupervised pesudotime order within the three clusters of single cells (Figure 4a). In consistence with the differentiating order predicted by Cytotrace algorithm, the mesenchymal progenitors had the highest differential potential, and the differentiate routine starting from mesenchymal progenitors to enthesoblasts and finally to enthesis chondrocytes (Figure 4b and Figure 4—figure supplement 1). To illustrate the gene ontology changes of BTJ chondrocytes differentiation between different timepoints, we next performed time-dependent Gene Ontology (GO) enrichment analysis. The biological processes alike extracellular matrix organization, collagen fibril organization, and regulation of cellular response to growth factors stimulus were significantly upregulated with enthesis developmental time increasing to postnatal day 14, 28, and 56 (Figure 4c and Figure 4—figure supplement 1c).

Figure 4 with 1 supplement see all
The trajectory and gene dynamics in bone-tendon junction (BTJ) chondrocytes differentiation.

(a) Pseudotime analysis of the three clusters of mesenchymal progenitors, enthesoblasts, and enthesis chondrocytes. (b) Cytotrace scores and predicted ordering of the three subclusters. (c) Heatmap revealing the scaled expression of differentially expressed genes (DEGs) and their annotated function for each timepoint. (d) Heatmap showing the modules of coregulated genes grouped by Louvain community analysis. (e) Uniform Manifold Approximation and Projection (UMAP) plots and functional annotation of the representative gene modules, showing the top five gene ontology (GO) annotations of indicated biological process among each timepoint. (f) Line plots showing representative gene trends in modules 15 and 22, respectively. (g) Feature plots showing top driving genes in BTJ chondrocyte growth. (h) Heatmap showing the gene expression dynamics along differentiation trajectories of BTJ chondrocytes. (i) Spatial transcriptomic spot map reveals the expression of driving genes in each spatial spot. (j) Immunofluorescence shows distribution and dense level of Tnn and Mfge8 during the different postnatal times. Scale bars, 100 μm. N=3-4. *p<0.05, **p<0.01.

To explore gene expression dynamics along the BTJ chondrocytes differentiation trajectories, we next examined gene patterns that varied with BTJ chondrocytes into 25 modules using Louvain community analysis. We found genes in Module 15 and Module 22 were expressed preferentially in BTJ chondrocytes, which were annotated for cartilage development, ossification, and biomineralization (Figure 4e). In contrast, stemness-related genes (Ly6a, Cd34, and Cd44) were downregulated along the trajectories, whose gene modules annotated for mesenchyme morphogenesis, regulation of organ formation, and regulation of mesenchymal cell proliferation (Figure 6—figure supplement 1).

In light of the gene dynamics along the postnatal BTJ chondrocytes differentiation, we measured the gene expression in pseudotime. The heatmap and the gene feature plots show the putative genes which driving BTJ chondrocyte differentiation, as well as the top driving genes predicted by cellrank (Figure 4g and h). Among the most highly significant driving genes (Klf9, Fxyd5, Klf4, Mfge8, Sox9, Clec3a, Wwp2, Tnn), we then confirmed the expression of Klf9, Klf4, Clec3a, Mfge8, and Tnn in single-cell spatial transcriptomics, except for previously reported Sox9 and Wwp2 which relative to chondrogenesis (Blitz et al., 2013). We found the expression of Mfge8 and Tnn had not been reported, and we validated the expression of Mfge8 and Tnn proteins in enthesis, as the expression of Tnn increased and reached the top at postnatal week 2 (P14), and the expression of Mfge8 increased and maintained a higher level after birth (Figure 4j).

Biological processes and regulators governing enthesis chondrocytes differentiation

To determine cellular functions across varied cell subpopulations or development stages, we then compared the GO ontology across major cell clusters. As expected, BTJ chondrocytes were enriched with genes annotated for cartilage development, extracellular matrix organization, and biomineral tissue development, and enthesoblasts were enriched with fibroblast proliferation and developmental growth involved in morphogenesis (Figure 5a and b). The time-dependent analysis of each GO annotation in BTJ chondrocytes showed that biological processes alike chondrocyte proliferation, development, and biomineral tissue development decreased with time increased to P56. While the activity of collagen synthesis maintained a steady level from after birth to postnatal days 28 and 56 (Figure 5c).

Biological processes and regulators governing enthesis chondrocytes differentiation.

(a) Heatmap shows the typical biological processes enriched in each cellular cluster. (b) Dot plots show the gene ontology (GO) clusters in bone-tendon junction (BTJ) chondrocytes in enthesis development. (c) Time-dependent analysis of GO annotations including chondrocyte proliferation, development, and biomineral tissue development decreased with time increased. (d) Heatmap shows the most significant regulatory regulons in each subcluster. (e, g) Uniform Manifold Approximation and Projection (UMAP) plots and spatial expression of Sox9 and Mef2a/Mef2c regulons and their target genes in enthesis chondrocytes. (f, h) Immunofluorescence shows distribution and dense level of Sox9 and Mef2a/Mef2c during different postnatal time. Scale bars, 100 μm.

We performed single-cell regulatory network inference and clustering (SCENIC) to investigate the gene regulatory networks that might govern enthesis growth (Figure 5d). The results showed that previously reported enthesis-related regulon Sox9 was significantly expressed both in BTJ tenocytes and chondrocytes, in comparison with other cell types. We then checked the Sox9 expression in single-cell dataset and spatial transcriptomics, the results showed that Sox9 regulon is mostly expressed in chondrocytes and some parts of tenocytes adjunct to enthesis. Target gene analysis showed the downstream targets of Sox9 regulon were mostly genes known associated with cartilage development (Acan, Scrg1, Hapln1, Chad) (Figure 5e). In-vivo validation results showed that in an early stage of postnatal growth (P1 and P7), Sox9-positive cells were widely located in the tendon enthesis and humeral head. And the expression of Sox9 decreased with time increased, partly visible in the enthesis cell (Figure 5f). In addition, we checked the expression of Mef2a/Mef2c, which had been reported relative to biomineralization and chondrocyte hypertrophy (Chen et al., 2023; Leupin et al., 2007). Immunofluorescence revealed that these Mef2a/Mef2c positive populations were found at the SST enthesis, and were more abundant at the mid-stage of enthesis differentiation (P14-P28), consistent with the emergence of fibrochondrocytes observed in histological stainings (Figure 5h).

Intercellular crosstalk signaling networks regulating enthesis postnatal growth

To seek further insights into the critical factors which may regulate the enthesis postnatal growth, we refined the CellChat and cellphoneDB input for downstream analysis (including seven clusters of BTJ chondrocytes, BTJ tendons, Tenocytes, Osteocytes, Enthesoblasts, and Mesenchymal progenitors) and performed the signaling communication analysis. Both CellChat and cellphoneDB results identified the aggregated signaling network for intercellular crosstalk. Relative active bidirectional signaling interactions among these cell subclusters revealed highly regulated cellular communications (Figure 6a). We then identified the signaling roles of each subcluster, the results showed that the cluster of BTJ chondrocytes predominately showed incoming patterns, as suggestive of signaling receivers (Figure 6b). We further identified signal components that contribute most to the incoming signaling among all these subclusters (Figure 6c and d).

Figure 6 with 1 supplement see all
Intercellular crosstalk signaling networks regulating enthesis postnatal growth.

(a) Overview of the cellular network regulating the postnatal enthesis growth predicted by CellChat and Cellphone DB. (b) The dominant senders (sources) and receivers (targets) among seven cell clusters. (c) Identify signals contributing most to the incoming signaling of bone-tendon junction (BTJ) chondrocytes. (d) Heatmap shows the most significant signaling networks among each subcluster. (e, g) Overview of FGF and BMP signalings networks in enthesis development. Hierarchy plots show the inferred signaling networks among all cell clusters. Heatmaps show the signaling roles of cell groups. Bar plots show the ligand-receptor pairs contributed significantly to BTJ chondrocytes. Feature plot shows the validation of Bmpr2 and Fgfr2 in spatial transcriptomic data. (f, h) Validations Fgfr2 and Bmpr2 protein in enthesis by immunofluorescence stainings. Scale bars, 100 μm.

FGF and BMP signalings are crucial to articular cartilage development, yet their specific roles in enthesis development need further investigation. To determine the important factors, we analyzed the intercellular signaling networks of FGF and BMP signalings. First, the expression pattern of FGF-FGFR signaling was noticed, as both autocrine and paracrine in BTJ chondrocyte. BTJ chondrocytes, tenocytes, and mesenchymal progenitors were leading senders of FGF signaling (Figure 6e). We observed Fgfr2 expressed mostly in BTJ chondrocyte cells via Fgf2-Fgfr2 or Fgf9-Fgfr2, and we validated the expression of FGFR2 protein in enthesis by immunofluorescence staining (Figure 6f, and Figure 6—figure supplement 1). The transforming growth factor-β (TGF-β) superfamily includes a family of proteins, such as TGF-βs (TGF-β1, TGF-β2, and TGF-β3) and bone morphogenetic proteins (e.g. BMP2, BMP4). In the BMP signaling network, the BTJ chondrocytes acted as critical receivers and contributors by secreting BMP ligand Bmpr2, especially (Figure 6g and Figure 6—figure supplement 1). In immunofluorescence results, BMPR2 was observed highly expressed in enthesis from P7 to P56. Ligand-receptor analysis points to Bmp2 and Bmp4 sent from mesenchymal progenitors and enthesoblasts to BTJ chondrocytes, suggesting the important role of BMP signaling in enthesis differentiation (Figure 6—figure supplement 1). In addition, the communication network of TGF-β and PTH signaling pathways was checked, as BTJ chondrocytes majorly received the stimulation by TGF-β1, TGF-β3, and Pthlh, which had been reported positive in chondrocyte differentiation and biomineralization, respectively (Bobzin et al., 2021; Xiao et al., 2022; Figure 6—figure supplement 1).

Discussion

Deciphering how a complex enthesis is formed from fetal-like into mature status may shed light on the strategies for enhanced BTJ healing (Derwin et al., 2018; Xiao et al., 2022). However, the cellular complexity and heterogeneity of developing RC enthesis are poorly understood, as previous studies could hardly resolve it at the single-cell level (Zhang et al., 2022). So far, there is only one transcriptomic study for embryonic mouse enthesis has been carried out using single-cell RNA sequencing (Kult et al., 2021). However, the development of fibrocartilage at the enthesis of mouse RC occurs no earlier than 2 weeks after birth (Galatz et al., 2007), suggesting the investigation at the postnatal stage of cellular and genomic mechanisms in enthesis development is needed. According to the works performed by Fang F et al., (Fang et al., 2022), they found enthesis progenitors (Gli1+ progenitors) and validated their stemness in-vitro and in-vivo, within the timepoints from P11 to P56. In this study, we applied single-cell transcriptome analysis to delineate the comprehensive postnatal enthesis growth with temporal atlas from as early as postnatal day 1 up to postnatal day 56. We next used the spatial transcriptome sequencing on postnatal day-1 mice enthesis to verify the anatomical position of different cell populations. This study may facilitate a better understanding of the enthesis development and add to the single-cell datasets repository of enthesis.

According to prior studies, three distinct populations appear where supraspinatus tendon attaches to the humeral head cartilage: tendon midsubstance progenitors, enthesis progenitors, and primary cartilage progenitors (Blitz et al., 2013; Dyment et al., 2015). Enthesis morphogenesis involves predominately enthesis progenitors transforming into fibrocartilage, during which process enthesis progenitors are organized as an unmineralized cartilaginous attachment unit and then mineralizes via endochondral ossification postnatally (Galatz et al., 2007). We used H&E staining to characterize the morphological changes of the cellular components of enthesis after birth. We found that fibrocartilage did not appear at the BTJ site in mice until 2–3 weeks after birth, and the statistical results showed that enthesis cells’ size increased significantly during postnatal development at 7–14 days after birth (Figure 1d). At 14 days after birth, these cells showed a typical chondrocyte phenotype. This is consistent with the findings of Galatz et al., 2007. According to the literature, the development of fibrocartilage in the enthesis occurs postnatally, which is not evident until 2 weeks after birth (Galatz et al., 2007). We consider that the postnatal 7–14 days were the critical stages for fibrocartilage cell differentiation in enthesis, based on the results of our cell morphology studies. Although there is no qualitative evidence for cells with specific gene markers, morphological changes can reflect changes in cell composition and function, which is a very important characteristic change in the study of the development of cell populations. Consistent with the timing we found, Schwartz AG et al., reported that Gli1-expressing cells significantly increased and populated at the enthesis site since postnatal day 7, well before the onset of mineralization, and persisted in the mature enthesis (Schwartz et al., 2015). Gli1+ cells and their progenies are retained in the enthesis region throughout postnatal development, contributing hugely to enthesis growth (Felsenthal et al., 2018; Jensen et al., 2018; Schwartz et al., 2015). Fei F et al., also focused on Gli+ progenitor cells, their data provided important clues to questions related to the development of the enthesis. In order to obtain sequencing data over a broader time span, we combined our own data with that of Fei F et al. And, we further mapped BTJ chondrocytes onto spatial positions by spatial transcriptome sequencing and found that they were consistent with the anatomical positions of enthesis, which verifies the accuracy of our definition of cell populations.

It is generally recognized that difficulties in restoring the mechanical properties after BTJ injury are largely due to the failure of fibrocartilage recapitulation (Shengnan et al., 2021). Moreover, in modern RC reconstruction surgeries, anchors fix the tendon to the insertion areas without access to bone marrow, which means BMSCs are not likely the main stem cell source for repair (Bi et al., 2007; Schwartz et al., 2015; Utsunomiya et al., 2013). Therefore, it is needed to investigate the native cellular origination and molecular biology of fibrocartilage formation, in order to enlighten developmental engineering strategies. To investigate the trajectory of postnatal bone-tendon junction cell development. Cellrank analysis was performed, and results showed that the directionality of differentiation between BTJ chondrocytes, tenocytes, and two subsets of osteocytes were independent of each other. The fate probabilities analysis, consistent with RNA velocity, showed that only BTJ chondrocytes contributed to the fate of enthesis chondrocytes (Figure 3e and f), suggesting that the fibrocartilage in postnatal enthesis origin from enthesis site-specific progenitors, instead of tendon cells. We also noticed that Tnn was one of the driving genes in chondrogenic fibroblasts differentiating into fibrochondrocytes, and we confirmed the existence of the tenascin N protein (also named tenascin W) in the developing fibrocartilage layer. We still know very little about the basic biology of tenascin W, which has been reported to be expressed in developing and mature bone, specifically in a subset of stem cell niches (Meloty-Kapella et al., 2006). Details are scarce, but the stimulating effects of tenascin-W on osteoblastic progenitors’ differentiation and migration have been reported (Meloty-Kapella et al., 2008; Morgan et al., 2011), suggesting its potential role in facilitating enthesis progenitors differentiating into fibrochondrocytes.

To investigate the key factors that regulate enthesis development, CellChat analysis was performed. We found that enthesis cells mainly received signals from other cell types, instead of sending signaling factors. And, we focused on the growth factors signaling pathways that were involved in the enthesis cells network, mainly including the previously reported FGF family (Bobzin et al., 2021; Roberts et al., 2019), BMP family (Blitz et al., 2013), TGF-β family (Tan et al., 2020; Xiao et al., 2022), and PTH family (Felsenthal et al., 2018; Schwartz et al., 2015; Figure 6 and Figure 6—figure supplement 1). Among them, we first found that FGF signaling was widely expressed in enthesis cells. As previously mentioned, the development of fibrocartilage in enthesis is similar to a growth plate with an endochondral-like zone. According to the literature, pre-hypertrophic cells in the growth plate express high levels of Fgfr3, and hypertrophic chondrocytes express high levels of Fgfr1 (Ornitz and Itoh, 2015). Yet we found enthesis cells mainly expressed Fgfr2, suggesting that the regulation of FGFs in enthesis progenitors differentiating into fibrocartilage cells was different from growth plate cartilage. Recent work has confirmed that enthesis development in the mouse mandible was regulated by FGF signaling via FGFR2-FGF2 signaling (Roberts et al., 2019). These findings suggest a potential role of FGFR2 in enthesis cells differentiating into fibrocartilage during enthesis development. We next noticed the BMP signaling (specifically Bmp2, Bmp4) was expressed in enthesis cells, as one key feature of the Sox9 and Scx positive progenitors is their dependence on Bmp2 and Bmp4 for specification and differentiation (Blitz et al., 2013; Bobzin et al., 2021). Blitz et al., found that BMP4 derived from the tendon tip induces enthesis progenitors differentiating into chondrocytes, and conditional inactivation of Bmp4 using Scx-Cre blocks formation of the cartilage anlage prefiguring the bone eminence (Blitz et al., 2009). Interestingly, it has been indicated that BMP2 could upregulate tenascin-N expression through a p38-dependent signaling pathway (Scherberich et al., 2005), and we found Tnn was among the top driving genes in fibrochondrocyte differentiation. We also found the TGF-β signaling in enthesis cells, as TGF-β was important due to its crucial role in cartilage and tendon development (Killian, 2022). Canonical TGF-β ligands may be diffuse into the near tendon and enthesis, positive in recruiting chondrogenic cells, and the secretion of TGF-β1 has been confirmed mechanically mediated (Subramanian et al., 2018; Xiao et al., 2022).

There are some certain limitations in the current study. First, we could not remove all the humeral head cells because the enthesis tissues only contain 4–5 layers of cells and are located adjunct to bone marrow and growth plate. Future use of high-precision microdissection approaches to isolate region-specific cells will address the limitation. Second, it is undeniable that spatial transcriptomics are better reliable to address possible dissociation artifacts and gain spatial information. However, utilization of spatial transcriptomic sequencing on enthesis is limited owing to the difficulty in sectioning without decalcification, which restricted our attempt to acquire spatial transcriptomics on mature enthesis with tough bony tissue. Finally, this study was performed predominately on single-cell RNA and spatial transcriptomics datasets, despite we verified the molecules inferred by the analysis algorithm, the whole study was designed as descriptive research, Future studies will label the markers found in this study on transgenic mice and investigate their in-vivo function in enthesis development.

In summary, our study deciphered the cellular complexity and heterogeneity of postnatal enthesis growth by providing descriptive single-cell transcriptomic and spatial datasets. We then revealed the molecular dynamics during fibrocartilage differentiation, providing a valuable resource for further investigation of tendon enthesis development at the mechanistic level, which may facilitate a better understand of the enthesis development and add to the single-cell datasets repository of the enthesis.

Materials and methods

Collection of cells from the supraspinatus tendon enthesis

Request a detailed protocol

All animal experimental protocols were approved by the Animal Ethics Committee of Central South University (No. 2022020058). The humeral head- supraspinatus tendon samples were dissected from the left shoulders of C57/BL6 mice at postnatal day-1, day-7, day-14, and day-28. In general, samples were harvested from pooled sibling limbs of two litters (five to six limbs per pool). Following dissection, the humeral heads and tendons were trimmed to retain the enthesis part, and all the samples were minced immediately and digested in type I collagenase (1 mg/ml, Gibico) and type II collagenase (1 mg/ml, Gibico) diluted in low-glucose DMEM (Gibco) solution at 37 °C for 30–40 min. Freshly isolated cells were resuspended into FACS buffer containing 2% FBS (Gbico) in PBS. Cell suspensions were stained with antibodies including Ter119-Alexa700 and Cd45- Alexa700 (Biolengend) to remove blood cells. DAPI (BD) stain was used to exclude dead cells. Flow cytometry was performed on BD FACS Aria II, single cells were gated using doublet-discrimination parameters and collected in FACS buffer.

Single-cell spatial transcriptomic sequencing by stereo-seq

Request a detailed protocol

Single-cell spatial transcriptomic sequencing was performed on the BGI stereo-seq platform (Chen et al., 2022). Briefly, the tissue section of the postnatal day-1 left shoulder of C57/bl6 mice was placed on the Stereo-seq chip (1 cm * 1 cm), then incubated and stained with a mixture of nucleic acid reagent (Invitrogen, Q10212). Section images were captured using a Zeiss Axio Scan Z1 microscope (at EGFP wavelength, 10 ms exposure). Tissue sections were then permeated to release RNAs from the permeated tissue and captured by a Stereo-seq chip. RNAs were then reverse transcribed and the cDNA-containing chips were then amplified with Hot Start DNA Polymerase (QIAGEN). In the library preparing procedure, samples were tagmented with Tn5 transposases (Vazyme) and amplified. After amplification, the PCR products were used for DNB (DNA Nano Ball) generation. Finally, the DNBs were sequenced on the DNBSEQTM T10 sequencing platform (MGI, Shenzhen, China).

Spatial mapping of cell states with cell2location

Request a detailed protocol

Cell2location was used to deconvolute and map single-cell clusters onto spatial transcriptomics spots. In brief, we first estimated reference signatures of cell states using scRNA-seq data of each region and a negative binomial regression model provided in the cell2ocation package (Kleshchevnikov et al., 2022). The regression model for the single-cell data was initialized with default settings. The model was then trained using a maximum of 30,000 epochs. The inferred reference cell type signatures were used for cell2location cell-type mapping for corresponding regions that estimate the abundance of each cell state in each spot.

Droplet-based scRNA-seq

Request a detailed protocol

8000–10,000 cells were loaded for each age group by Chromium instrument and its chemistry kit V3 (10 X Genomics) according to the manufacturer’s guidance. Each cell was encapsulated with a barcoded Gel Bead in a single partition, then amplified to generate single-cell cDNA libraries and sequenced on an Illumina NovaSeq 6000 platform at a sequencing depth of ~500 million reads. The Cellranger pipeline (version 6.1.1) was used to align the raw reads to the mouse reference genome GRCm38 and to generate feature-barcode matrices. All the low-quality reads were filtered with default parameters.

Single-cell data processing, quality control, and integration

Request a detailed protocol

All the feature-barcode matrices were loaded by the Seurat package (v4.1.0) (Hao et al., 2021), doublets or cells with poor quality were removed (less than 200 genes and greater than 2 Median absolute deviations above the median, or more than 5% genes mapping to the mitochondrial genome). After quality control, all the feature data were scaled with the sctransform algorithm, to avoid unwanted variation including percentages of mitochondrial reads, number of detected genes, and predicted cell cycle phase effect. All the datasets were integrated and batch-corrected by using SCVI with default parameters. Furthermore, this integrated data was analyzed and subclustered to exclude uninterested clusters (including immune cells, red blood cells, endothelial cells, smooth muscle cells, and neural cells).

Dimensionality reduction, clustering, and DEGs analysis

Request a detailed protocol

We used the Uniform Manifold Approximation and Projection (UMAP) and Potential of Heat diffusion for Affinity-based Transition Embedding (PHATE) (Moon et al., 2019) method to visualize the dataset in low dimensions. Furthermore, the K-nearest neighbor (KNN) method and the Louvain algorithm were applied to cluster the cells, with 50 PCs selected and resolution set to 2.4, resulting in nine major cell clusters for subsequent analyses. For second-round chondrocyte sub-clustering, we reconstructed the SNN graphs for BTJ clusters, and three subclusters were determined resolution set as 0.6 for each fibrochondrocyte cluster. The FindAllMarkers function in Seurat was used to calculate DEGs among different clusters, the ‘test.use’ function was set to a statistical framework called MAST (Finak et al., 2015). Genes met the criteria that (1) expressing a minimum fraction of 10% in either of the two tested populations; (2) at least a 0.1-fold difference (log-scale) between the two tested populations; (3) adjusted p values less than 0.01, were considered as signature genes. Clusters were annotated according to the expression of those highly variable genes reported in the literature.

Time-dependent gene signature clustering

Request a detailed protocol

DEGs between different timepoints were acquired FindAllMarkers function in Seurat, temporal pattern analysis, and visualization were conducted by using R package Tcseq according to a standard pipeline. Through the above algorithm, the time-dependent DEGs were divided.

Trajectory analysis and cell state analysis

Request a detailed protocol

Before trajectory analysis, the S4 Seurat object was transformed into an anndata object using the Seuratdisk package and loaded by Scanpy (Wolf et al., 2018). Then, all the bam files were processed with Velocyto (La Manno et al., 2018) to quantify the spliced and unspliced mRNA counts. Subsequently, the Velocyto outputs were loaded into scVelo (Bergen et al., 2020) and merged with the anndata object from Scanpy to compute RNA velocity vectors. Low abundance genes (less than 30 total counts) were filtered from the merged dataset. After RNA velocity analysis, we used Cellrank (Lange et al., 2022) package to compute infer the terminal cell state and cluster absorption probabilities using nearest-neighbour relationships and RNA velocity with equal weight in CellRank’s Markov chain model.

Cell-cell interaction analysis

Request a detailed protocol

Cell-cell interaction analysis was performed using CellChat (Jin et al., 2021) package, according to a standard pipeline.

Gene regulatory network analyses

Request a detailed protocol

We applied Single-Cell Regulatory Network Inference and Clustering (SCENIC) (Aibar et al., 2017) to identify the cluster-specific gene regulatory networks. The pySCENIC grn method was performed for building the initial co-expression gene regulatory networks (GRN). The regulon data was then analyzed using the RcisTarget package to create TF motifs referring to the mm9-tss-centered-10kb-7 database. The regulon activity scores were calculated using the Area Under the Curve (AUC) method. Besides, we used Cellcall package to analyze the cluster-specific TF enrichment and intercellular communication by combining the expression of ligands/receptors and downstream TF activities for certain L-R pairs. Genes that were expressed in less than 10% of the cells of a certain cell type were excluded.

GO enrichment analysis

Request a detailed protocol

GO enrichment of cluster differentially expressed genes was performed using the R package clusterProfiler (Wu et al., 2021), with a Benjamini–Hochberg (BH) multiple testing adjustment and a false-discovery rate (FDR) cutoff of 0.1. The Gene Ontology Resource database (http://geneontology.org) was used for GO pathway analysis. Module scores for each gene set were calculated using the AddModuleScore function implemented in Seurat. Gene sets used for scoring (Chondrocyte proliferation, Proteoglycans synthesis, Cartilage homeostasis, Collagen synthesis, Regulation of bone development, Biomineralization, Negative regulation of bone mineralization) were selected from the Gene Ontology Browser of MGI Database (C5: biological process gene sets). Visualization was performed using the R package ggplot2.

Sample harvest and histological observation

Request a detailed protocol

The left shoulder of the C57/BL6 mice was harvested on postnatal days 1, 3, 7, 14, 28, and 56. Specimens were obtained and fixed with 10% formalin buffer for 24 hr and rinsed by dual evaporated water then gradually dehydrated by sequential immersion in 70%, 80%, 90%, and 100% alcohol (each for 2 hr), finally dried in the air before use. The samples were embedded in paraffin and then sectioned for histological studies with H-E and Toluidine blue/Fast green staining. Histologic sections were observed using light microscopy (CX31, Olympus, Germany).

Immunofluorescence staining

Request a detailed protocol

The left shoulder of the C57/BL6 mice was harvested on postnatal days 1, 3, 7, 14, 28, and 56. Then fixed with 4% neutral buffered formalin for 24 hr. After decalcifying, dehydrating, and embedded in OCT, specimens were longitudinally sectioned with 10 µm. For immunofluorescence staining, the sections were washed with PBS, permeabilized with 0.1% TritonX-100, and then blocked with 5% bovine serum albumin (BSA; Sigma-Aldrich). Sections were incubated with primary antibodies anti-Sox9 (Abcam, ab185996), anti-Scx (Santa Cruz, sc518082), anti-Tnn (Thermo Fisher, ab_2900654), anti-Mef2a/Mef2c (Abclonal, A2710), and anti-Mfge8 (Abclonal, A12322) at 4 °C overnight, then incubated with Alexa-Fluor 488 conjugated secondary antibody (Abcam, ab150129) and Alexa-Fluor 594(Abcam, ab150120) conjugated secondary antibody at room temperature for 1 hr and counterstained with DAPI (Invitrogen, USA). All the images were observed and captured using a Zeiss AxioImager.M2 microscope (Zeiss, Germany) equipped with an Apotome.2 System. Densities of Sox9, Scx, Tnn, Mef2a/Mef2c, or Mfge8 positive cells of each captured image were measured using 200 x magnification graphs for each slide by the Image J software (National Institutes of Health, Bethesda, MD). The antibodies used in this study were listed in Supplementary file 1.

Statistical analysis

Request a detailed protocol

One-way ANOVA and Student’s t-test were performed to assess whether there were statistically significant differences in the results between time groups. Values of p<0.05 were considered to be significantly different. Data were analyzed using Prism 7 software (GraphPad).

Data availability

All single-cell datasets created during this study are publicly available at the Gene Expression Omnibus (GSE223751). Any additional information required to re-analyze the data in the paper is available from the corresponding author upon request.

The following data sets were generated
    1. Zhang T
    2. Wan LY
    (2022) NCBI Gene Expression Omnibus
    ID GSE223751. Single-cell RNA sequencing reveals cellular and molecular heterogeneity in fibrocartilaginous enthesis formation.
The following previously published data sets were used
    1. Fang F
    2. Thomopoulos S
    (2022) NCBI Gene Expression Omnibus
    ID GSE182997. Single-cell RNA-seq of the tendon enthesis cells.

References

    1. Chen C
    2. Wu X
    3. Han T
    4. Chen J
    5. Bian H
    6. Hei R
    7. Tang S
    8. Li Z
    9. Lu Y
    10. Gu J
    11. Qiao L
    12. Zheng Q
    (2023)
    Mef2a is a positive regulator of Col10a1 gene expression during chondrocyte maturation
    American Journal of Translational Research 15:4020–4032.
    1. Ornitz DM
    2. Itoh N
    (2015) The fibroblast growth factor signaling pathway
    Wiley Interdisciplinary Reviews. Developmental Biology 4:215–266.
    https://doi.org/10.1002/wdev.176

Decision letter

  1. Murim Choi
    Senior and Reviewing Editor; Seoul National University, Republic of Korea
  2. Xiao Chen
    Reviewer; Zhejiang University School of Medicine, China

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Single-cell RNA sequencing reveals cellular and molecular heterogeneity in fibrocartilaginous enthesis formation" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Murim Choi as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Xiao Chen (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

Although the editors and reviewers agreed that the work is valuable, it needs further improvements to warrant publication in eLife. Major issues are:

1) The paper should be discussed relative to the results in Fang et al., Cell Stem Cell, 2022.

2) The importance of P7 as a critical differentiation timepoint is not well supported.

3) Validation is needed (e.g., IHC and/or FISH) for many of the scRNAseq results (e.g., CellChat pathways).

Reviewer #1 (Recommendations for the authors):

1. The results and their novelty should be discussed in comparison to the recent Cell Stem Cell study describing enthesis development using scRNAseq and lineage tracing approaches (https://doi.org/10.1016/j.stem.2022.11.007).

2. Figure 1d: The PCA for histomorphologic parameters (which typically have high variations) does not show any meaningful separations or groupings. I do not agree with the authors' conclusion that this analysis reveals P7 as a critical developmental timepoint. In general, PCA may not be appropriate for this data set.

3. According to the methods, it appears that the entire humeral head-supraspinatus tendon was used for cell isolation for scRNAseq. This results in the inclusion of cells from bone, growth plate, enthesis and tendon. As such, only a very small percentage of cells came from the enthesis, as is clear from the cell clusters in Figure 2b and 2c. This is a flaw in the approach; inclusion of such a wide range of cells makes interpretation of "enthesis" cells difficult, as described in more detail in the comments below.

4. The differentiation/pseudotime analysis described in Figure 3 is difficult to follow. I do not think it is useful to combine cell transcriptomes from vastly different tissues and then define a velocity map. There is too much varied information for the algorithm to create valid connections, as the there will be many many branches/paths from mesenchymal stem cell to osteoblast, tenocyte, chondrocyte, etc. Presumably, embedded in these maps are trajectories for osteoblast differentiation, chondrocyte differentiation, tenocyte differentiation, etc. There are too many layers of overlapping information to deduce anything meaningful for the small number of cells associated with the enthesis.

5. The authors uses the term "function" throughout the paper (e.g., "functional definition of fibrocartilage subpopulations"). However, this is a descriptive study, and "function" (or mechanism) can only be theoretically inferred from the various algorithms used to analyze the data. A role for any of the pathways or processes can only be defined with gain- and/or loss-of-function studies.

6. "C2 highly expressed biomineralization-related genes (Clec3a, Tnn, Acan)". The three example genes are not related to biomineralization.

7. The functional characterization of the three enthesis cell clusters is not convincing. For example, activation of metabolism-related processes is a vague result than can mean a lot of things (including changes in differentiation), yet the authors interpret it very specifically as " role in postnatal fibrochondrocyte formation and growth".

8. The pseudotime analysis of the enthesis cell clusters is not convincing. The three clusters are quite close and overlapping on the UMAP. Furthermore, the authors focus on Tnn as a novel and unique gene, yet the expression pattern shown in Figure 5g implies even expression of this gene across all three clusters.

9. The TC1 markers (Ly6a, Dlk3, Clec3b) imply a non-tendon-specific cell population. Perhaps a tendon progenitor pool or an endothelial cell phenotype is more appropriate.

10. Pseudotime analyses assume that your data set includes cells from progenitor through mature cell populations. It is unclear that the timepoints studied here included cells from early progenitor states.

11. The CellChat analysis is not useful, as the authors included 18 cell types. The number of possible interactions among so many cell types is enormous, and deducing valid connections between any two cell types in this case is questionable.

12. The authors should validate their key scRNAseq results with in situ hybridization. Only a single gene, Tpp, was validated on tissue sections. This validation is particularly important for this study because the authors included a wide variety of tissues/cells in their isolation and analyses.

13. The authors should demonstrate functional necessity of at least one gene/pathway identified by the scRNAseq analyses (e.g., through gene knockout).

Reviewer #2 (Recommendations for the authors):

1. As known, Fei Fang et al. have established single-cell transcriptomes of mouse supraspinatus tendon enthesis cells (Cell Stem Cell, 2022). It is suggested that the authors introduced Fei Fang et al.'s work in Introduction and emphasize the significant novelty compared with Fei Fang et al.'s work.

2. In Figure 1, the authors highlighted P7 was a critical stage for enthesis differentiation. But this section was less associated with the following content. The authors should link these results with the scRNASeq data. Is there any time-dependent change/signaling in scRNASeq data at this critical time point?

3. In the H and E staining of Figure 1A, the tendon structure was separated and random. It is suggested that the authors provide high-quality staining figures.

4. Figure 2 showed that the Scx+ or Sox9+ cells was decreased in enthesis over time. At least it should be co-staining to show the distribution and frequency of double positive and single positive cell populations. However, a previous study has demonstrated this finding (PLOS ONE, 2020). It is suggested to verify some new findings by IF or IHC staining.

5. There are some conflicts about trajectory analysis. In Figure 3C, RNA velocity showed that the arrow flowed from BTJ to MTJ and CTFb. However, in Fig3d, PAGA plot indicated that BTJ cells is independent of other cells. Furthermore, in supplementary figure S3, RNA velocity showed that the trajectory flowed from TC to BTJ. These figures were inconsistent with the described results. Please provide detailed explanation to avoid misleading readers.

6. Figure 5 showed that C1 was the original cluster, and whether C1 cluster expressed canonical progenic/stem cell markers.

7. The authors performed cell-cell interaction based on cellchat analysis. But the cell-cell interaction was not actively examined.

Reviewer #3 (Recommendations for the authors):

1. Fang et al. (A mineralizing pool of Gli1-expressing progenitors builds the tendon enthesis and demonstrates therapeutic potential. Cell stem cell. 2022) defined enthesis cell transcriptomes and differentiation trajectories, and identified Gli1+ progenitor population for enthesis. Please further clarify the innovation of your research, and in depth introduction or discussion is needed to compare and contrast the results between the two research.

2. In Figure 1, more evidence are needed to prove that Neonatal to postnatal day 7 (P7) is the critical stage for enthesis fibrocartilage cell differentiation, for example, immunofluorescence staining or qPCR for enthesis fibrocartilage cell makers, instead of relying on H and E only.

3. Line123. Figure 2e showed that the expression of Clec3a and Col2a1 were low in c4," which were ubiquitously expressed in bone-tendon junction cell (c4)" seems to be an inexact expression.

4. Line 117, which cell clusters belong to "fibroblast-associated cells"?

5. Line 125, it is better to co-staining the scx and Sox9 to validate the existence of BTJ cells. Scx and Sox9 are known markers of BTJ, do you have find new makers for BTJ by scRNA-seq?

6. Line 148, "stemness" degree? Are there other evidence, such as stem cell maker expression, to show that "growth plate cells and fibroblasts associated clusters are higher than other cell types". The expression of "stemness" seems exaggerated.

7. There is no description of figure 4b in the results.

8. In figure 5, 2-3 makers identified by scRNA-seq for fibrocartilage formation are suggested to be validated by immunofluorescence stainning or other methods, instead of only proving the Tnn expression in postnatal BTJ growth.

9. There are no verification of the signaling network for the enthesis postnatal growth which were revealed by Cellchat. It is suggested to validate the key signaling, such as Bmpr2 signaling.

https://doi.org/10.7554/eLife.85873.sa1

Author response

Essential revisions:

Although the editors and reviewers agreed that the work is valuable, it needs further improvements to warrant publication in eLife. Major issues are:

1) The paper should be discussed relative to the results in Fang et al., Cell Stem Cell, 2022.

2) The importance of P7 as a critical differentiation timepoint is not well supported.

3) Validation is needed (e.g., IHC and/or FISH) for many of the scRNAseq results (e.g., CellChat pathways).

Thank you for your kind suggestions. According to your comments, we carefully revised the paper, and the issues mentioned in prior version has been addressed:

(1) We discussed the relation and consistency between our dataset and Fang et al. in writing. Meanwhile, we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al. Consequently, we are now better able to discuss the enthesis development from P1 to P58, including 7 timepoints to ideally cover the whole postnatal enthesis growth.

(2) We removed the statement that P7 was a critical differentiation timepoint in revised paper. And we furtherly added toluidine blue staining to observe the extracellular matrix of fibrocartilage. At the same time, we performed immunochemistry to test the protein level of collagen II at enthesis.

(3) We performed immunofluorescence staining to validate the results from cell trajectory inference (Tnn and Mfge8 expression, Figure 4), transcription factor calculation (Mef2a/2c and Sox9 expression, Figure 5), and communication results (Fgf2r and Bmp2r expression, Figure 6). And we performed spatial transcriptome sequencing on enthesis slide at postnatal day 1 to exam the location of these results.

Reviewer #1 (Recommendations for the authors):

1. The results and their novelty should be discussed in comparison to the recent Cell Stem Cell study describing enthesis development using scRNAseq and lineage tracing approaches (https://doi.org/10.1016/j.stem.2022.11.007).

Thank you for your suggestions. We discussed the relation and consistency between our dataset and Fang et al. (revised in Figure 2, Figure 2—figure supplement 2 and lines 75-79). Meanwhile, we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang F et al. Consequently, we are now better able to discuss the enthesis development from P1 to P58, including 7 timepoints to ideally cover the whole postnatal enthesis growth.

2. Figure 1d: The PCA for histomorphologic parameters (which typically have high variations) does not show any meaningful separations or groupings. I do not agree with the authors' conclusion that this analysis reveals P7 as a critical developmental timepoint. In general, PCA may not be appropriate for this data set.

Thank you for your suggestions. We removed the PCA results in revised paper. And we realized the improper statement about P7 as a critical developmental timepoint. We furtherly added toluidine blue staining to observe the extracellular matrix of fibrocartilage. At the same time, we performed immunochemistry to test the protein level of collagen II at enthesis.

We found fibrocartilage was not evident at the enthesis of mouse rotator cuff until 2-3 weeks after birth, as evidenced by toluidine blue stained cartilage ECM and IHC validated Col2a1 expression. Yet we found the cell sizes of enthesis cells remarkably increased during postnatal day 7-14, and the enthesis cells were observed as typical chondrocyte phenotype, column-like stacked alongside the direction of tendon fiber. These findings led us the conclusion that the development of fibrocartilage in the enthesis occurs postnatally. (revised in Figure 1 and lines 94-104)

3. According to the methods, it appears that the entire humeral head-supraspinatus tendon was used for cell isolation for scRNAseq. This results in the inclusion of cells from bone, growth plate, enthesis and tendon. As such, only a very small percentage of cells came from the enthesis, as is clear from the cell clusters in Figure 2b and 2c. This is a flaw in the approach; inclusion of such a wide range of cells makes interpretation of "enthesis" cells difficult, as described in more detail in the comments below.

Thank you for your suggestions, and we are sorry for the misinterpretation about sample preparation in prior manuscription.

In our practice, the humeral head- supraspinatus tendon samples were dissected from the left shoulders of C57/BL6 mice at postnatal day-1, day-7, day-14, and day-28. In general, samples were harvested from pooled sibling limbs of two litters (five to six limbs per pool). Following dissection, the humeral heads and tendons were trimmed to retain the enthesis part, and all the samples were minced immediately and digested. (revised in Lines 456-458)

Although, we can’t deny the fact that we could not remove all the humeral head cells under microscope because the enthesis tissues only contained 4-5 layers of cells and located adjunct to bone marrow and growth plate. The enthesis datasets inevitably contains some cells came from bone marrow and growth plate in humeral head, as shown by the work reported by Fang F et al. (Figure 1a, Cell Stem Cell. 2022) (revised in Lines 406-409)

In order to better understand the preciseness of our cluster and annotation results, we furtherly performed single cell spatial transcriptomic sequencing on postnatal day 1 mice enthesis ice section, we then deconvoluted and mapped P1 single-cell clusters (BTJ chondrocytes) onto spatial transcriptomics spots, and the cell2location results showed that BTJ chondrocytes we had annotated in single cell dataset could be properly mapped onto enthesis site (Figure 2e). (revised in Lines 134-136)

However, utilization of spatial transcriptomic sequencing on enthesis are limited owing to the difficulty in sectioning without decalcification, which restricted our attempt to acquire spatial transcriptomics on mature enthesis with tough bony tissue. (revised in Lines 409-413)

4. The differentiation/pseudotime analysis described in Figure 3 is difficult to follow. I do not think it is useful to combine cell transcriptomes from vastly different tissues and then define a velocity map. There is too much varied information for the algorithm to create valid connections, as the there will be many many branches/paths from mesenchymal stem cell to osteoblast, tenocyte, chondrocyte, etc. Presumably, embedded in these maps are trajectories for osteoblast differentiation, chondrocyte differentiation, tenocyte differentiation, etc. There are too many layers of overlapping information to deduce anything meaningful for the small number of cells associated with the enthesis.

Thank you for your suggestions. We realized the bad effects caused by too many layers of overlapping information in our prior results, and we re-analyzed the differentiation pseudotime across cell types.

Because we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al., thus we re-clustered and re-annotated the cell clusters. We first performed Cytotrace and Cellrank analysis across all cell clusters to predict cellular differential potential and terminal differentiate states. And we identified differentiation trajectories of these 3 terminal cell types (revised in Figure 3 and lines 155-184).

Our next step was performed by refining the input for pseudotime analysis (including mesenchymal progenitors, enthesoblasts, and enthesis chondrocytes), and we focused on BTJ chondrocyte differentiation trajectory to reconstruct the gene dynamics (revised in Figure 4 and lines 197-225).

5. The authors uses the term "function" throughout the paper (e.g., "functional definition of fibrocartilage subpopulations"). However, this is a descriptive study, and "function" (or mechanism) can only be theoretically inferred from the various algorithms used to analyze the data. A role for any of the pathways or processes can only be defined with gain- and/or loss-of-function studies.

Thank you for your suggestions. We realized the improper statement of functionality in the prior version of paper. We scrutinized the paper and corrected the statement. (revised in Figure 5 and lines 241-249)

6. "C2 highly expressed biomineralization-related genes (Clec3a, Tnn, Acan)". The three example genes are not related to biomineralization.

Thank you for your suggestions. In combination with GSE182997 datasets provided by Fang et al., we re-analyzed the clustering outcomes of BTJ chondrocytes and we realized it is inappropriate to subcluster the current BTJ chondrocytes. And we removed the results related to the three enthesis cell subclusters in revised version. At the same time, we recalculated the biomineralization-related genes (Runx2, Ibsp, Spp1, Col11a1) in the revised manuscription. (revised in Figure 3h and 4f)

7. The functional characterization of the three enthesis cell clusters is not convincing. For example, activation of metabolism-related processes is a vague result than can mean a lot of things (including changes in differentiation), yet the authors interpret it very specifically as " role in postnatal fibrochondrocyte formation and growth".

Thank you for your suggestions. Relating to question 6, we re-analyzed the clustering outcomes of BTJ chondrocytes and we believed it is inappropriate to subcluster the current BTJ chondrocytes. And we removed the results related to the three enthesis cell subclusters in revised version.

8. The pseudotime analysis of the enthesis cell clusters is not convincing. The three clusters are quite close and overlapping on the UMAP. Furthermore, the authors focus on Tnn as a novel and unique gene, yet the expression pattern shown in Figure 5g implies even expression of this gene across all three clusters.

Thank you for your suggestions. We re-performed differentiation pseudotime analysis (including mesenchymal progenitors, enthesoblasts, and enthesis chondrocytes), and we focused on BTJ chondrocyte differentiation to reconstruct the gene dynamics. (revised in Figure 4 and lines 197-225)

We found the expression of BTJ chondrocytes differentiate driving genes Mfge8 and Tnn had not been reported, and we furtherly validated the expression of Mfge8 and Tnn in enthesis by immunofluorescence staining. (revised in Figure 4)

9. The TC1 markers (Ly6a, Dlk3, Clec3b) imply a non-tendon-specific cell population. Perhaps a tendon progenitor pool or an endothelial cell phenotype is more appropriate.

Thank for your suggestions. We decided not to discuss the tendon fate in this study, in consideration of the tissue digest and cell isolation methods was designed for extracting enthesis chondrocytes, instead of tendon cells. And we have removed the results and discussion in association with tendon subpopulations.

10. Pseudotime analyses assume that your data set includes cells from progenitor through mature cell populations. It is unclear that the timepoints studied here included cells from early progenitor states.

Thank you for your suggestions. According to the works performed by Fang F et al. (Cell Stem Cell. 2022), they found enthesis progenitors (Gli1+ progenitors) and validated their stemness in-vitro and in-vivo, within the timepoints from P11 to P56 (revised in lines 325-328).

From enthesis histological results, enthesis specific fibrochondrocytes are observable no early than postnatal 14 days, suggested that the enthesis chondrocyte formation is a postnatal progress. We integrated our dataset with GSE182997 datasets (3 sample) provided by Fang F et al. Consequently, we discussed the enthesis development from P1 to P58, including 7 timepoints to ideally cover the whole postnatal enthesis growth (revised in lines 324-331).

11. The CellChat analysis is not useful, as the authors included 18 cell types. The number of possible interactions among so many cell types is enormous, and deducing valid connections between any two cell types in this case is questionable.

Thank you very much for your suggestions. We refined the cellchat input for downstream analysis, including 7 clusters of enthesis-related cells (BTJ chondrocytes, BTJ tendons, Tenocytes, Osteocytes, Enthesoblasts, and Mesenchymal progenitors). Then we checked the expression and location of CellChat inferred target genes in spatial transcriptomic data. Meanwhile we performed immunofluorescence staining to check the Fgfr2 and Bmpr2 expression inferred by cellchat. (revised in Figure 6 and lines 286-303)

12. The authors should validate their key scRNAseq results with in situ hybridization. Only a single gene, Tpp, was validated on tissue sections. This validation is particularly important for this study because the authors included a wide variety of tissues/cells in their isolation and analyses.

Thank you for your suggestions. validate the results from cell trajectory inference (Tnn and Mfge8 expression, Figure 4), transcription factor calculation (Mef2a/2c and Sox9 expression, Figure 5), and communication results (Fgf2r and Bmp2r expression, Figure 6). And we performed spatial transcriptome sequencing on enthesis slide at postnatal day 1 to exam the location of these results in enthesis.

13. The authors should demonstrate functional necessity of at least one gene/pathway identified by the scRNAseq analyses (e.g., through gene knockout).

Thank you very much for your suggestions. We aimed to provide a transcriptional resource for further investigation of fibrocartilage development.

We found expression of Tnn, Clec3a, and Mfge8 as significant driving genes in enthesis chondrocytes differentiation, suggesting their possibilities as new makers for BTJ chondrocytes, and we checked their expression in spatial transcriptomic data. Furtherly we stained the protein expression of Tnn and Mfge8 to validate their location.

The whole study was designed as descriptive research, Future studies will label the markers found in this study on transgenic mice and investigate their in-vivo function in enthesis development. (revised in Figure 4 and lines 412-415)

Reviewer #2 (Recommendations for the authors):

1. As known, Fei Fang et al. have established single-cell transcriptomes of mouse supraspinatus tendon enthesis cells (Cell Stem Cell, 2022). It is suggested that the authors introduced Fei Fang et al.'s work in Introduction and emphasize the significant novelty compared with Fei Fang et al.'s work.

Thank for your suggestions, we discussed the relation and consistency between our dataset and Fang et al. Meanwhile, we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al. Consequently, we hope to be better able to discuss the enthesis development from P1 to P58, including 7 timepoints to ideally cover the whole postnatal enthesis growth.

2. In Figure 1, the authors highlighted P7 was a critical stage for enthesis differentiation. But this section was less associated with the following content. The authors should link these results with the scRNASeq data. Is there any time-dependent change/signaling in scRNASeq data at this critical time point?

Thank for your suggestions, we removed the PCA results in revised paper. And we furtherly added toluidine blue staining to observe the extracellular matrix of fibrocartilage. At the same time, we performed immunochemistry to test the protein level of collagen II at enthesis.

We found fibrocartilage was not evident at the enthesis of mouse rotator cuff until 2-3 weeks after birth, as evidenced by toluidine blue stained cartilage ECM and IHC validated col2a1 expression. Yet we found the cell sizes of enthesis cells remarkably increased during postnatal day 7-14, and the enthesis cells were observed as typical chondrocyte phenotype, column-like stacked alongside the direction of tendon fiber. These findings led us the conclusion that the development of fibrocartilage in the enthesis occurs postnatally. (revised in Figure 1 and lines 90-104)

And we summarized the time-dependent gene dynamic and bioprocess change in BTJ chondroctyes differentiation. (revised in Figure 4 and lines 197-225)

3. In the H and E staining of Figure 1A, the tendon structure was separated and random. It is suggested that the authors provide high-quality staining figures.

Thank for your suggestions, and we restained the sections with H-E and toluidine blue to better show the structure change in enthesis growth. At the same time, we performed immunochemistry to test the protein level of collagen II at enthesis. (revised in Figure 1a)

4. Figure 2 showed that the Scx+ or Sox9+ cells was decreased in enthesis over time. At least it should be co-staining to show the distribution and frequency of double positive and single positive cell populations. However, a previous study has demonstrated this finding (PLOS ONE, 2020). It is suggested to verify some new findings by IF or IHC staining.

Thank for your suggestions. We co-stained Scx+ and Sox9+ and recalculated Scx and Sox9 positive cells in enthesis. We also checked Scx and Sox9 expression in our P1 spatial transcriptomic data. (revised in Figure 2e, 2f and lines 136-141)

In differentiation pseudotime analysis, we identified the most highly significant driving genes (Klf9, Fxyd5, Klf4, Mfge8, Sox9, Clec3a, Wwp2, Tnn), we then confirmed the expression of Klf9, Klf4, Clec3a, Mfge8, and Tnn in single cell spatial transcriptomics, except for previously reported Sox9 and Wwp2 which relative to chondrogenesis (Blitz et al., 2013). We found the expression of Mfge8 and Tnn had not been reported, and we validated the expression of Mfge8 and Tnn proteins in enthesis. (revised in Figure 4 and lines 216-225)

5. There are some conflicts about trajectory analysis. In Figure 3C, RNA velocity showed that the arrow flowed from BTJ to MTJ and CTFb. However, in Fig3d, PAGA plot indicated that BTJ cells is independent of other cells. Furthermore, in supplementary figure S3, RNA velocity showed that the trajectory flowed from TC to BTJ. These figures were inconsistent with the described results. Please provide detailed explanation to avoid misleading readers.

Thank you for your suggestions. Because we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al., thus we re-clustered and re-annotated the cell clusters. Therefore, we re-analyzed the differentiation pseudotime across different cell types.

We first performed Cytotrace and Cellrank analysis across all cell clusters to predict cellular differential potential and terminal differentiate states. And we identified differentiation trajectories of these 3 terminal cell types (revised in Figure 3 and lines 163-175).

Our next step was performed by refining the input for pseudotime analysis (including mesenchymal progenitors, enthesoblasts, and enthesis chondrocytes), and we focused on BTJ chondrocyte differentiation to reconstruct the gene dynamics (revised in Figure 4 and lines 209-225).

6. Figure 5 showed that C1 was the original cluster, and whether C1 cluster expressed canonical progenic/stem cell markers.

Thank you for your suggestions. In combination with GSE182997 datasets (3 sample) provided by Fang et al., we re-analyzed the clustering outcomes of BTJ chondrocytes and we believed it is inappropriate to subcluster the current BTJ chondrocytes. And we removed the results related to the three enthesis cell clusters in revised version.

We compared canonical progenic stem cell markers (Ly6a, Cd34, and Cd44) between BTJ chondrocytes, enthesoblasts and mesenchymal progenitors (revised in Figure 4—figure supplement 1).

7. The authors performed cell-cell interaction based on cellchat analysis. But the cell-cell interaction was not actively examined.

Thank for your suggestions, we refined the cellchat input for downstream analysis, including 7 clusters of enthesis-related cells (BTJ chondrocytes, BTJ tendons, Tenocytes, Osteocytes, Enthesoblasts, and Mesenchymal progenitors). Then we checked the expression and location of CellChat inferred target genes in spatial transcriptomic data. Meanwhile we performed immunofluorescence staining to check the Fgfr2 and Bmpr2 expression inferred by cellchat. (revised in Figure 6 and lines 286-303)

Reviewer #3 (Recommendations for the authors):

1. Fang et al. (A mineralizing pool of Gli1-expressing progenitors builds the tendon enthesis and demonstrates therapeutic potential. Cell stem cell. 2022) defined enthesis cell transcriptomes and differentiation trajectories, and identified Gli1+ progenitor population for enthesis. Please further clarify the innovation of your research, and in depth introduction or discussion is needed to compare and contrast the results between the two research.

Thank you for your suggestions. We discussed the relation and consistency between our dataset and Fang et al. Meanwhile, we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al. Consequently, we are now better able to discuss the enthesis development from P1 to P58, including 7 timepoints to ideally cover the whole postnatal enthesis growth. (revised in Figure 2 and lines 115-133)

2. In Figure 1, more evidence are needed to prove that Neonatal to postnatal day 7 (P7) is the critical stage for enthesis fibrocartilage cell differentiation, for example, immunofluorescence staining or qPCR for enthesis fibrocartilage cell makers, instead of relying on H and E only.

Thank for your suggestions, we removed the PCA results in revised paper. And we realized the improper statement about P7 as a critical developmental timepoint. We furtherly added toluidine blue staining to observe the extracellular matrix of fibrocartilage. At the same time, we performed immunochemistry to test the protein level of collagen II at enthesis. (revised in Figure 1 and lines 86-104)

3. Line123. Figure 2e showed that the expression of Clec3a and Col2a1 were low in c4," which were ubiquitously expressed in bone-tendon junction cell (c4)" seems to be an inexact expression.

Thank you for your suggestions. In combination with GSE182997 datasets (3 sample) provided by Fang et al., we re-analyzed the clustering outcomes of BTJ chondrocytes and we believed it is inappropriate to subcluster the current BTJ chondrocytes. And we removed the results related to the three enthesis cell clusters in revised version.

4. Line 117, which cell clusters belong to "fibroblast-associated cells"?

Thank you for your suggestions. We are sorry for the misinterpretation of cell annotations. Because we integrated our dataset with open source GSE182997 datasets (3 sample) provided by Fang et al., thus we re-clustered and re-annotated the cell clusters. And the clustering descriptions had been revised in the current version paper. (revised in Figure 2 and lines 115-141)

5. Line 125, it is better to co-staining the scx and Sox9 to validate the existence of BTJ cells. Scx and Sox9 are known markers of BTJ, do you have find new makers for BTJ by scRNA-seq?

Thank you for your suggestions. We co-stained Scx+ and Sox9+ and recalculated Scx and Sox9 positive cells in enthesis. We also checked Scx and Sox9 expression in our P1 spatial transcriptomic data. (revised in Figure 2e, 2f and lines 134-141)

Meantime, we found expression of Tnn, Clec3a, and Mfge8 was significantly upregulated in enthesis chondrocytes differentiation, suggesting their possibilities as new makers for BTJ chondrocytes, and we checked their expression in spatial transcriptomic data. Furtherly we stained the protein expression of Tnn and Mfge8 to validate their location (revised in Figure 4 and lines 216-225).

In future study, we will label these markers on transgenic mice and investigate their in-vivo function in enthesis development. (revised in lines 412-416).

6. Line 148, "stemness" degree? Are there other evidence, such as stem cell maker expression, to show that "growth plate cells and fibroblasts associated clusters are higher than other cell types". The expression of "stemness" seems exaggerated.

Thank you for your suggestions. We re-analyzed the datasets and we filtered out doublets, dead, and apoptotic cells, blood cells (erythrocytes and progenitors), endothelial cells, immune cells (B cells and T cells), myeloid cells, and growth plate chondrocytes.

The downstream analysis in revised version was focused on the developmental trajectories for tenocytes, chondrocytes, and osteocytes differentiation in enthesis. (revised in Figure 2, 3 and lines 155-225)

7. There is no description of figure 4b in the results.

Thank you for your suggestions. We corrected the mistake in revised manuscription. (revised in Figure 4 and lines 200-203)

8. In figure 5, 2-3 makers identified by scRNA-seq for fibrocartilage formation are suggested to be validated by immunofluorescence stainning or other methods, instead of only proving the Tnn expression in postnatal BTJ growth.

Many thanks for your suggestions. We performed immunofluorescence staining to validate the results from cell trajectory inference (Tnn and Mfge8 expression, figure 4), transcription factor calculation (Mef2a/2c and Sox9 expression, figure 5), and communication results (Fgf2r and Bmp2r expression, figure 6). And we also performed spatial transcriptome sequencing on enthesis slide at postnatal day 1 to exam the location of these results in enthesis. (revised in Figure 4 and lines 216-225)

9. There are no verification of the signaling network for the enthesis postnatal growth which were revealed by Cellchat. It is suggested to validate the key signaling, such as Bmpr2 signaling.

Thank you for your suggestions, we refined the cellchat input for downstream analysis, including 7 clusters of enthesis-related cells (BTJ chondrocytes, BTJ tendons, tenocytes, osteocytes, enthesoblasts, and mesenchymal progenitors). Then we checked the expression and location of CellChat inferred target genes in spatial transcriptomic data. Meanwhile we performed immunofluorescence staining to check the Fgfr2 and Bmpr2 expression inferred by cellchat. (revised in Figure 6 and lines 286-303)

https://doi.org/10.7554/eLife.85873.sa2

Article and author information

Author details

  1. Tao Zhang

    1. Department of Sports Medicine, Xiangya Hospital Central South University, Changsha, China
    2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft
    Contributed equally with
    Liyang Wan
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6746-8834
  2. Liyang Wan

    1. Department of Sports Medicine, Xiangya Hospital Central South University, Changsha, China
    2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Visualization, Methodology, Writing – original draft
    Contributed equally with
    Tao Zhang
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2194-1080
  3. Han Xiao

    1. Department of Sports Medicine, Xiangya Hospital Central South University, Changsha, China
    2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  4. Linfeng Wang

    1. Department of Sports Medicine, Xiangya Hospital Central South University, Changsha, China
    2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Data curation, Validation
    Competing interests
    No competing interests declared
  5. Jianzhong Hu

    1. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    2. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    3. Department of Spine Surgery and Orthopaedics, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Conceptualization, Supervision, Project administration
    Competing interests
    No competing interests declared
  6. Hongbin Lu

    1. Department of Sports Medicine, Xiangya Hospital Central South University, Changsha, China
    2. Key Laboratory of Organ Injury, Aging and Regenerative Medicine of Hunan Province, Changsha, China
    3. National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
    Contribution
    Conceptualization, Supervision, Funding acquisition, Project administration
    For correspondence
    hongbinlu@hotmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7749-3593

Funding

National Natural Science Foundation of China (82230085)

  • Hongbin Lu

National Natural Science Foundation of China (82272572)

  • Hongbin Lu

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors would like to thank Professor Hui Xie, Xiang-Hang Luo, and other staff from the Movement System Injury and Repair Research Center, Xiangya Hospital, Central South University, Changsha, China, for their kind assistance during the experiments. Funding this study was supported by the National Natural Science Foundation of China (NO. 82230085, and 82272572).

Ethics

All animal experimental protocols were approved by the Animal Ethics Committee of Central South University (No. 2022020058).

Senior and Reviewing Editor

  1. Murim Choi, Seoul National University, Republic of Korea

Reviewer

  1. Xiao Chen, Zhejiang University School of Medicine, China

Version history

  1. Received: December 30, 2022
  2. Preprint posted: February 3, 2023 (view preprint)
  3. Accepted: September 10, 2023
  4. Accepted Manuscript published: September 12, 2023 (version 1)
  5. Version of Record published: September 21, 2023 (version 2)

Copyright

© 2023, Zhang, Wan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 774
    Page views
  • 130
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Tao Zhang
  2. Liyang Wan
  3. Han Xiao
  4. Linfeng Wang
  5. Jianzhong Hu
  6. Hongbin Lu
(2023)
Single-cell RNA sequencing reveals cellular and molecular heterogeneity in fibrocartilaginous enthesis formation
eLife 12:e85873.
https://doi.org/10.7554/eLife.85873

Share this article

https://doi.org/10.7554/eLife.85873

Further reading

    1. Developmental Biology
    Nicolas G Brukman, Clari Valansi, Benjamin Podbilewicz
    Research Article Updated

    The fusion of mammalian gametes requires the interaction between IZUMO1 on the sperm and JUNO on the oocyte. We have recently shown that ectopic expression of mouse IZUMO1 induces cell-cell fusion and that sperm can fuse to fibroblasts expressing JUNO. Here, we found that the incubation of mouse sperm with hamster fibroblasts or human epithelial cells in culture induces the fusion between these somatic cells and the formation of syncytia, a pattern previously observed with some animal viruses. This sperm-induced cell-cell fusion requires a species-matching JUNO on both fusing cells, can be blocked by an antibody against IZUMO1, and does not rely on the synthesis of new proteins. The fusion is dependent on the sperm’s fusogenic capacity, making this a reliable, fast, and simple method for predicting sperm function during the diagnosis of male infertility.

    1. Biochemistry and Chemical Biology
    2. Developmental Biology
    Sima Stroganov, Talia Harris ... Michal Neeman
    Research Article Updated

    Background:

    Fetal growth restriction (FGR) is a pregnancy complication in which a newborn fails to achieve its growth potential, increasing the risk of perinatal morbidity and mortality. Chronic maternal gestational hypoxia, as well as placental insufficiency are associated with increased FGR incidence; however, the molecular mechanisms underlying FGR remain unknown.

    Methods:

    Pregnant mice were subjected to acute or chronic hypoxia (12.5% O2) resulting in reduced fetal weight. Placenta oxygen transport was assessed by blood oxygenation level dependent (BOLD) contrast magnetic resonance imaging (MRI). The placentae were analyzed via immunohistochemistry and in situ hybridization. Human placentae were selected from FGR and matched controls and analyzed by immunohistochemistry (IHC). Maternal and cord sera were analyzed by mass spectrometry.

    Results:

    We show that murine acute and chronic gestational hypoxia recapitulates FGR phenotype and affects placental structure and morphology. Gestational hypoxia decreased labyrinth area, increased the incidence of red blood cells (RBCs) in the labyrinth while expanding the placental spiral arteries (SpA) diameter. Hypoxic placentae exhibited higher hemoglobin-oxygen affinity compared to the control. Placental abundance of Bisphosphoglycerate mutase (BPGM) was upregulated in the syncytiotrophoblast and spiral artery trophoblast cells (SpA TGCs) in the murine gestational hypoxia groups compared to the control. Hif1α levels were higher in the acute hypoxia group compared to the control. In contrast, human FGR placentae exhibited reduced BPGM levels in the syncytiotrophoblast layer compared to placentae from healthy uncomplicated pregnancies. Levels of 2,3 BPG, the product of BPGM, were lower in cord serum of human FGR placentae compared to control. Polar expression of BPGM was found in both human and mouse placentae syncytiotrophoblast, with higher expression facing the maternal circulation. Moreover, in the murine SpA TGCs expression of BPGM was concentrated exclusively in the apical cell side, in direct proximity to the maternal circulation.

    Conclusions:

    This study suggests a possible involvement of placental BPGM in maternal-fetal oxygen transfer, and in the pathophysiology of FGR.

    Funding:

    This work was supported by the Weizmann Krenter Foundation and the Weizmann – Ichilov (Tel Aviv Sourasky Medical Center) Collaborative Grant in Biomedical Research, by the Minerva Foundation, by the ISF KillCorona grant 3777/19.