Single-Cell Integration Analysis of Heterotopic Ossification and Fibrocartilage Developmental Lineage: Endoplasmic Reticulum Stress Effector Xbp1 Transcriptionally Regulates the Notch Signaling Pathway to Mediate Fibrocartilage Differentiation

Introduction Regeneration of fibrochondrocytes is essential for the healing of the tendon-bone interface (TBI), which is similar to the formation of neurogenic heterotopic ossification (HO). Through single-cell integrative analysis, this study explored the homogeneity of HO cells and fibrochondrocytes. Methods This study integrated six datasets, namely, GSE94683, GSE144306, GSE168153, GSE138515, GSE102929, and GSE110993. The differentiation trajectory and key transcription factors (TFs) for HO occurrence were systematically analyzed by integrating single-cell RNA (scRNA) sequencing, bulk RNA sequencing, and assay of transposase accessible chromatin seq. The differential expression and enrichment pathways of TFs in heterotopically ossified tissues were identified. Results HO that mimicked pathological cells was classified into HO1 and HO2 cell subsets. Results of the pseudo-temporal sequence analysis suggested that HO2 is a differentiated precursor cell of HO1. The analysis of integrated scRNA data revealed that ectopically ossified cells have similar transcriptional characteristics to cells in the fibrocartilaginous zone of tendons. The modified SCENIC method was used to identify specific transcriptional regulators associated with ectopic ossification. Xbp1 was defined as a common key transcriptional regulator of ectopically ossified tissues and the fibrocartilaginous zone of tendons. Subsequently, the CellPhoneDB database was completed for the cellular ligand-receptor analysis. With further pathway screening, this study is the first to propose that Xbp1 may upregulate the Notch signaling pathway through Jag1 transcription. Twenty-four microRNAs were screened and were found to be potentially associated with upregulation of XBP1 expression after acute ischemic stroke. Conclusion A systematic analysis of the differentiation landscape and cellular homogeneity facilitated a molecular understanding of the phenotypic similarities between cells in the fibrocartilaginous region of tendon and HO cells. Furthermore, by identifying Xbp1 as a hub regulator and by conducting a ligand–receptor analysis, we propose a potential Xbp1/Jag1/Notch signaling pathway.


Introduction
Tendon to bone healing is a major concern in sports medicine. Regardless of the advances in rotator cuff repair and anterior, healing in the tendon-bone interface (TBI) is the chief problem that has attracted significant attention [1][2][3][4].
The TBI contains four coherent layers with gradual composition, structure, and mechanical properties. Histologically, the TBI is composed of bone, calcified fibrocartilage, noncalcified fibrocartilage, and tendon, which ensures the smooth transmission of force from the muscle to bone [5][6][7][8]. However, following rotator cuff repair or anterior cruciate ligament reconstruction, large amounts of fibrovascular tissue accumulate in the TBI, hindering healing and weakening the biomechanical properties of the joint [6][7][8][9][10][11][12]. Accumulating evidence has demonstrated that fibrochondrocytes localized within the TBI are essential for tendon to bone healing [13]. Rui et al. found that the tendon and bone healed well in rabbits with more fibrocartilage cells, and the promotion of fibrocartilage cell proliferation could enhance tendon-bone healing [8]. Therefore, the study of the origin and characteristics of fibrochondrocytes could help improve tendon to bone healing [13].
Neurogenic heterotopic ossification is a common complication of acute ischemic stroke; however, its exact mechanism remains unknown. Heterotopic ossification (HO) is the appearance of osteoblasts in the soft tissues and the formation of bone tissue [14,15]. Interestingly, HO demonstrates several characteristics similar to the healing process of the TBI. First, the tissue in HO is histologically similar to the TBI [14][15][16]. As an important feature of the tissue in HO and an essential process for fibrocartilage development, mineralization is observed in both HO and TBI formation [6,16,17]. Second, mechanical stimulation is pivotal for the development of both TBI and HO [16,18]. Third, bone marrowderived mesenchymal stem cells (BM-MSCs) and tendonderived stem cells, the cellular sources of fibrocartilage and HO, have been reported to possess a multidifferentiation potential [19][20][21][22][23][24][25]. Fourth, numerous studies have demonstrated that the Notch and Hedgehog (Hh) signaling pathways could drive ectopic ossification and tendon to bone healing by mediating associated transcription factors, although the exact underlying mechanisms remain unclear [20,[26][27][28][29][30]. Fifth, some genes, such as COL2A1 and MKX, were found to be associated with the upregulation of the Notch and Hh signaling pathways and have the potential to promote tendon to bone healing [27,[30][31][32][33]. Moreover, the transcript levels of these genes (including COL2A1 and MKX) were found to be upregulated in ectopic ossified tissues [27,34,35]. Therefore, considering the homogeneity between HO and the TBI, we hypothesized that research on HO might serve as the reference for a better understanding of tendon to bone healing. Exploring the key regulators and pathways associated with the development of HO might provide novel insights into fibrocartilage restoration in the TBI.
The single-cell technique was used in this study to investigate the homogeneity between HO and the TBI. The single-cell transcriptome, which is a recent, novel research method, enables comprehensive exploration of tissue microstates and thus the characterization and changes in the cellular condition [36][37][38][39]. The scRNA-seq transcriptome sample integration analysis has been used to assess cellular homology across different sequencing platforms at different centers [40,41]. We speculate that the homogeneity between the two cells can be explored by integrating the analysis of HO tissue and fibrocartilage cells by single-cell techniques.
Transcription factors are targets for the regulation of cell differentiation [42]. In recent years, significant progress was achieved in the construction of gene regulatory networks (GRNs) based on single-cell transcriptome expression data. This data would help us understand the transcriptional regulatory networks of cells and uncover key TFs [43,44]. On the other hand, transcription factormediated changes in cellular communication and pathway signaling are important manifestations of cellular homogeneity [45]. Cellular communication is defined as the transmission of information from one cell to another through a medium to produce a corresponding response. Ligandreceptor complex-mediated cell-cell communication is essential for coordinating various biological processes, such as development, differentiation, and inflammation [46]. CellPhoneDB serves as an integrated algorithm for identifying biologically relevant ligand-receptor interactions from the scRNA-seq data, which could be used to investigate new targets for disease treatment, while ligand proteins in biological signaling are regulated via transcription [47]. Numerous studies have demonstrated that ligand manipulation therapy is a viable method of targeted therapy [31]. Therefore, after constructing GRNs based on single-cell transcriptome expression data, we propose the use of that technique for ligand-receptor analysis, which may contribute to the identification of key transcription factors for cellular developmental targets in the HO and fibrocartilage regions, thus providing a theoretical foundation for future ligand therapy.
In this study, three single-cell sequencing datasets (GSE168153, GSE138515, and GSE102929) were available for the cells in the tendon fibrochondral region and HO tissue [46][47][48]. After the GRNs were constructed, ligandreceptor analysis was performed to identify the core transcription factors. We aimed to identify common intercellular differentiation and transcriptional regulation pathways associated with fibrocartilage development in the TBI and HO.

Materials and Methods
2.1. Data Retrieval. All data were obtained from the Gene Expression Omnibus database, which mainly includes the GSE94683, GSE144306, GSE168153, GSE138515, GSE102929, and GSE110993 datasets (Table 1). GSE94683 includes transcriptional profiles of mesenchymal stromal cells in ectopically ossified tissues from nine healthy donors and seven patients with HO. GSE168153 is the only singlecell sequencing dataset available for the fibrocartilaginous region of mouse tendons. GSE144306 and GSE168153 are datasets of the same study, and GSE144306 is an assay for transposase-accessible chromatin with high-throughput 2 Oxidative Medicine and Cellular Longevity sequencing (ATAC-seq) dataset. GSE138515 contains singlecell transcriptome data from the Achilles tendon of 6-weekold male C57Bl/6 mice. GSE102929 reveals single-cell RNA transcriptome data from 4-week-old mouse tail tendon cells under physiological and simulated HO pathological conditions. Lactating rat tail tendon has less extracellular matrix, more number of cells per tissue volume, and strong cell growth and differentiation ability. Primary and passaged tendon cells can be obtained in large quantities relatively quickly, facilitating subsequent experiments. Therefore, GSE102929 can be used to reveal the differentiation trajectory characteristics and heterogeneity of tendon-derived cells under pathological HO conditions. GSE110993 is a dataset reflecting microRNA changes following acute ischemic stroke. This dataset was used to screen for microRNAs with altered expression in the serum after acute ischemic stroke. The target mRNAs of the miRNAs were predicted using the miRWalk database (http://mirwalk.umm.uniheidelberg.de/).

Differential Expression Analysis of TFs.
The TTRUST dataset is a database that records the regulatory relationships of TFs and contains not only the target genes corresponding to TFs but also the regulatory relationships among TFs [49].
It was used to screen for differentially expressed TFs in HO tissues. Differential analysis was performed by the R software limma package, and p values less than 0.05 were defined as significantly different. Differentially expressed genes were ranked according to |logFC|, and 70 genes that were highly expressed in HO tissues were obtained.
2.3. Single-Cell RNA Sequencing (RNA-seq) Analysis. RNAseq analysis was performed using the Seurat package based on the R software (version 4.0.2) [50]. Canonical correlation analysis (CCA) was also used for dataset integration. Principal component analysis (PCA) and uniform manifold approximation and projection/t-distributed stochastic neighbor embedding (UMAP/TSNE) methods were used for the scRNA-seq dimensionality reduction analysis. Cell annotation of GSE168153 and GSE138515 was completed according to the information provided in the original article, and the CellMarker database and SingleR package were used for further correction [51,52]. The "FindAllMarkers" function in Seruat was further used for the differential expression analysis, as described in the previous analysis. Monocle 2 and Monocle 3 were used to perform pseudo-temporal trajectory analysis [53]. Cellular ligand-receptor analysis was performed according to the instruction manual using the Cell-PhoneDB database and software [54]. Raw data of all analysis processes were saved.

Hub TF Analysis.
To identify specific transcriptional regulators associated with ectopic ossification, a modified version of the SCENIC method was used to construct a regulatory network between HO-associated targets and TFs [55,56]. Thus, each TF has its corresponding downstream target gene. Currently, SCENIC only supports transcriptional positive regulation analysis. Based on the above results, the regulon activity score was calculated for each cell, and the method was repeated three times to confirm the stability of the regulatory relationships [57]. Then, by performing Pearson correlation coefficient analysis for each pair of regulatory relationships, the similarity of the overall regulatory activity of different TFs was assessed, and this was used to quantify the regulatory correlations among cell types. The top 5 genes in each cell were screened for correlation analysis.
2.5. ATAC-seq Analysis. The Fastq sequencing data for the ATAC-seq dataset of the tendon-bone union region were obtained from GSE144306. Filtering and quality control were performed with the trim_galore software, and bowtie2 was used for comparison and calculation of statistical comparison rates. After finding peaks using MACS2, the length of the insert fragments, fraction of reads in peak values, and irreproducibility discovery rate were used to calculate duplicates. The MEME tool was then used to annotate motifs and verify potential transcriptional regulatory relationships [58].  [61]. False discovery rate < 0:25 and p:adjust < 0:05 were considered significantly enriched.

Differential Expression and Enrichment Pathways of TFs in Ectopically Ossified
Tissues. Differential analysis was performed on 17076 genes from GSE94683 (including nine healthy donor samples and seven samples of ectopically ossified tissues from patients with HO). Among them, 1203 genes satisfied the threshold of |log 2ðFCÞ | >1 and p:adj < 0:05, 456 genes were highly expressed in HO tissues, and 747 genes were lowly expressed in HO tissues. Then, 70 genes with the highest expression in HO tissues were screened, and the differential expressions of these genes were visualized in different samples by drawing a heat map (Figure 1(a)). Moreover, 795 TFs identified were obtained from the TTRUST database. By plotting Venn diagrams, five genes, namely, GLI1, EN1, XBP1, SHOX, and TBX5, were found to be the intersection of the 795 TFs and 70 genes highly expressed in HO tissues (Figure 1(b)). The characteristics of these five TFs, namely, GLI1, EN1, XBP1, SHOX, and TBX5, differentially expressed in GSE94683 are shown in Figure 1(c). Thus, to explore the biological pathways differentially expressed in HO tissues, we performed GSEA. The analysis suggested that the WNT, Notch, Hh, and transforming growth factor (TGF) beta signaling pathways in the KEGG pathway were highly expressed in HO tissues (p < 0:05).

Heterogeneity and Proposed Temporal Differentiation
Trajectory of an Ectopic Ossification Cell Model. At present, the differentiation trajectory of ectopically ossified cells remains unknown, and the associated single-cell transcriptome has not been fully characterized. Silencing Mkx is shown as a viable method of constructing ectopic ossification models [15,64,65]. GSE102929 reveals scRNA transcriptomic data from 4-week-old mouse tail tendon cells under physiological and simulated HO pathological conditions, all filtered, normalized, and clustered. UMAP clustering analysis of single-cell RNA-seq revealed that single cells of mouse tail tendon cells could be distinguished into four clusters under physiological and simulated HO pathological conditions (Figure 2(a)). These cell clusters were found to correlate with the cell model groups. Among them, clusters 0 and 3 belong to the scRNA-seq data of mouse tail tendon cells under physiological conditions, while clusters 1 and 2 belong to the scRNA-seq data of mouse tail tendon cells under pathological conditions (Figure 2(b)). Furthermore, GSVA was used to quantify pathway activity within every single cell, all differential pathways (p < 0:05) were extracted, and a pathway heat map (Supplementary Figure 1A) was drawn according to the pathway activity of each sample. Cluster 2 has a more pronounced degree of HO differentiation than cluster 1. Therefore, the pathway activities of clusters 2 and 0 were compared to further identify the pathways that were  Hh and Notch signaling pathways, were highly expressed, and the oxidative phosphorylation pathway was also highly expressed in cluster 2 (Supplementary Figure 1E-1G). Mkx was significantly upregulated in cluster 4 (Supplementary Figure 2A, 2B). Then, heat maps were plotted to further demonstrate the degree of differential expression between the Notch and the Hh signaling pathways (Figure 2(c)).    Figure 2C). The results of the mimetic timing analysis suggested that HO2 was a precursor cell for HO1 differentiation (Figures 2(e) and 2(f)). To further characterize the pathological differentiation of ectopically ossified cells, the proposed chronological expression profile of genes associated with the ectopic ossification phenotype was revealed (Figure 2(g)). The gene expression profiles of Xbp1, E2f7, Fgfr1, Icam1, Jag1, and Ncam1 are shown in Supplementary Figure 2D. Mgp, Bgn, Sox9, Col11a1, Col1a1, Col1a2, Tnmd, Col11a2, Scx, Hivep2, Bcl2, Col2a1, Mkx, Myc, and Dcn were enriched in HO2-like cells. Col5a1, Foxc1, Gli3, Gli2, Pecam1, Klf2, Klf4, Zbtb48, Ptch2, Smad7, Sufu, Igfbp5, Tcf, Fbn1, and Jag1 were enriched mainly in HO1-like cells (Figure 2(g)). Among them, Gli3, Cli2, Pecam1, and Ptch2 are the star genes of the Hh signaling pathway, and Notch2 and Jag1 are the star genes of the Notch signaling pathway, which suggested that the activities of the Notch and Hh signaling pathways may be upregulated in HO1-like cells. According to the expression characteristics of these genes, HO1-like cells are more similar to chondrocyte-like cells, and HO2-like cells are more similar to tenocyte-like cells. Therefore, HO2 cells were defined as "tenocyte-like cells."

Endoplasmic Reticulum (ER) Stress Effector XBP1 Was
Found to Be a Key Transcriptional Regulator in Ectopically Ossified Tissues. In this study, key TFs in HO1-like cells (chondrocyte-like cells) and HO2-like cells (tenocyte-like cells) were identified using comprehensive network analysis (Supplementary Figure 3). The network analysis identified Ebf1, Usf2, Klf3, Zmiz1, and Sin3a as key transcriptional regulators in HO1-like cells, and Ep300, Smarcc, Foxj2, Xbp1, and E2f5 are key transcriptional regulators in HO2like cells. To systematically compare the transcriptional profiles in different cells, we compared the expression correlation of each TF based on the connection specificity index (CSI) (Figure 3(a)). In HO1-like cells (chondrocytelike cells) and HO2-like cells (tenocyte-like cells), Xbp1 and Klf3 were found to be key transcriptional regulators of MH1 module cells, and TF Ebf1 was found to be a key transcriptional regulator of MH2 module cells. In addition, the MH1 module was found to be more active in HO2-like (tenocyte-like cells) cells, while the MH2 module was more active in HO1-like (chondrocyte-like cells) cells (Figures 3(b) and 3(c)). This also demonstrates that Xbp1 and Klf3 are key transcriptional regulators in HO2 (tenocyte-like cells), while Ebf1 is a key transcriptional regulator in HO1 (chondrocytelike cells). The Venn diagram reveals that Xbp1 is a common key TF for HO tissue bulk-RNA-seq and scRNA-seq analyses (Figure 3(d)).

Ectopically Ossified Cells Have Similar Transcriptional
Characteristics to Cells in the Fibrocartilaginous Zone of Tendons. Based on the proposed temporal analysis characteristics of the gene expression in Figure 2(g), combined with the results of previous studies, the gene expression characteristics of ectopic ossification demonstrated strong similarities to those of the fibrocartilaginous region of the tendonbone union [8,[66][67][68]. Therefore, we hypothesized that ectopically ossified cells have similar transcriptional characteristics to cells in the fibrocartilaginous zone of the tendonbone junction region. Then, we integrated scRNA-seq transcriptome data from mouse Achilles tendon, scRNA-seq transcriptome data from mouse tendon junction region, and scRNA-seq transcriptome data from mouse tail tendon cells under simulated HO pathological conditions for scRNA-seq integration analysis.
First, the cell types of single-cell transcriptome data were annotated for GSE138515 and GSE168153. GSE138515 contains single-cell transcriptome data from the Achilles tendon of 6-week-old male C57Bl/6 mice, and PCA grouping was performed after quality control of these cells (Supplementary Figure 4A, 4C). With the UMAP method, GSE138515 was found to have 12 clusters of cell types (Figure 4(a) and Supplementary Figure 4B). The cell differentiation trajectory of GSE138515 is shown in Supplementary Figure 4D. A preliminary annotation of these cells was performed based on SingleR (Supplementary Figure 4E). Among the cell clusters, clusters 0, 7, and 11 were initially classified as endothelial cells; cluster 2 was initially classified as erythrocytes; clusters 1, 3, 4, 8, 9, 10, and 12 were initially classified as fibroblasts; cluster 5 was tentatively classified as macrophages; cluster 6 was tentatively classified as oligodendrocytes. In addition, these cells were further annotated regarding the results of previous studies (Figure 4(a)) [69]. As shown above, the quality of scRNAseq data of GSE168153 was controlled and grouped into nine clusters (Supplementary Figure 5A, 5D). The cell differentiation trajectory of GSE168153 is shown in Supplemental Figure 5B. Moreover, these cells were further annotated based on the results of previous studies ( Figure 4(b)) [66]. Accordingly, the proposed temporal expression profile of genes associated with the phenotypic cartilaginous region of tendon fibers is shown in Figure 4(c). The cell annotation method combines the results of the pseudo-temporal analysis to classify cells into tenocyte, tenocyte-like attachment cell, attachment cell, chondrocytelike attachment cell, and chondrocyte, which are the five major cell clusters (Figure 4(c) and Supplementary Figure 5C). The trajectory of cell differentiation is presented in the following order: tenocytes, tenocyte-like attachment cells, attachment cells, chondrocyte-like attachment cells, and chondrocytes (Supplementary Figure 5C). To the best of our knowledge, this study provided the most detailed subcellular classification of the tendon junctional area to date. Moreover, this classification demonstrates the spatial distribution characteristics of cells in the fibrocartilage layer. Previous analyses identified Xbp1 as a common key TF for HO tissue bulk-RNA-seq and scRNA-seq analyses. In the present study, the results of the pseudo-temporal analysis suggested that Xbp1 was enriched in tenocytes, tenocyte-like attachment cells, and attachment cells. "Chondrocyte-like attachment cells" and "Chondrocytes" were less expressed (Figures 4(d) and 4(e)). 8 Oxidative Medicine and Cellular Longevity CCA was also used for dataset integration. Three datasets (including GSE168153, GSE138515, and GSE102929 under simulated HO pathological conditions) were integrated using the CCA approach. These datasets were then grouped into nine cell clusters (Supplementary Figure 6A, 6B). The top 5 marker genes in each cell cluster were flagged (Figure 4(f) and Supplementary Figure 6B). In addition, the UMAP graph was used to show all cell clusters, and the original annotations of all cells were also identified (Supplementary Figure 6C). Figure 4(g) shows the distribution of datasets in the UMAP plot. Before reannotating the cells, HO2 was found to be enriched in attachment cell cluster 1 and tendon fibroblasts, while HO1 was enriched near chondrocyte-like attachment cell 3 (Figure 4(g) and Supplementary Figure 6C). After reannotation of the cells, HO2 was found to be enriched in tenocyte-like attachment cell 1, while HO1 was enriched in chondrocyte-like attachment cell 1 (Figure 4(h)). This suggests that ectopically ossified cells have similar transcriptional characteristics to cells in the fibrocartilaginous zone.
Cluster      For the review of cell fractionation, we found the relative proportions of the cell subpopulations, and the average cell numbers in each tissue are shown in Figure 6(f). Tenocytelike attachment cells and chondrocyte-like attachment cells 1 and 2 were found to be the major cell types in tendon tissues and mouse tail tendon cells under simulated HO pathological conditions. The core marker genes and sample preferences for each cluster are shown in Figure 6(g). Anxa2 Pseudo time      Figure 6(i). Then, the cell division cycles of all cells were explored using the Seruat package, and more significant heterogeneity in the mitotic cell divisions of each cell subpopulation was found (Figures 6(j) and 6(k)). These series of    Intersection analysis of these key TFs with previously studied TFs revealed that TF Xbp1 was again found to be a common key TF for HO bulk-RNA-seq, HO scRNA-seq, and the integrated scRNA data.
To explore the potential regulatory function of Xbp1 in the onset of ectopic ossification, the Venn diagram was used to screen 160 genes that collectively serve as downstream regulatory targets of XBP1 in the integrated scRNA data and HO scRNA data (Figure 8(a)). These 160 genes were subjected to GO and KEGG enrichment analysis (Figure 8(b)). GO enrichment analysis included three modules, namely, biological process (BP), cell component (CC), and molecular function (MF). The Golgi vesicle transport, protein localization to the ER, and ER to Golgi vesiclemediated transport were enriched in BP; the Golgiassociated vesicle membrane, COPI-coated vesicle, and Golgi apparatus part were enriched in CC, and the glucocorticoid receptor binding, primary active transmembrane transporter activity, and P-P-bond-hydrolysis-driven transmembrane transporter activity were enriched in MF. The pathways enriched in the KEGG were the circadian rhythm, protein export, and protein processing in the ER. In summary, genes that were downstream regulatory targets of Xbp1 were found to play important roles mainly in Golgi vesicle transport, ER protein processing, cellular communication, and transport. To visualize the interaction relationships of Xbp1 downstream regulatory target genes, string sites were used to map the protein-protein interaction (PPI) network analysis of XBP1 downstream regulatory genes.

Ligand-Receptor and ATAC-seq Analyses Reveal That
Xbp1 Upregulate the Notch Signaling Pathway by Promoting Jag1 Transcription. Cellular communication is an important pathway for inducing cell differentiation. We hypothesized that XBP1 may affect intercellular communication by regulating ligand transcription. To investigate the commonalities and differences in cellular signaling exchange for each cell subtype, the CellPhoneDB was used to infer intercellular communication. Based on the results of the CellPhoneDB analysis, five ligand genes (i.e., Fgfr1, Grin2b, Icam1, Jag1, and Ncam1) were identified as potential target genes of Xbp1 (Figure 9(a)). The Sankey plot was used to demonstrate the potential receptor genes (including AREG, BDNF, FGF23, FGF3, FGF5, FGF7, FGF8, FGF9, and FGFR2) for the five ligand genes (i.e., Fgfr1, Grin2b, Icam1, Jag1, and Ncam1), GDNF, IL16, KL, NCAM1, NOTCH2, and NOTCH3. Accordingly, the XBP1 functional pattern was mapped (Figure 9(c)). This suggests that XBP1 mediates the regulation of tenocyte-like attachment cell, chondrocytelike attachment cell 2, and chondrocyte-like attachment cell 1 phenotype by promoting the transcriptional regulation of these five ligand genes, namely, Fgfr1, Grin2b, Icam1, Jag1, and Ncam1, by receptor genes. The Pearson correlation analysis based on the GTEX database also revealed that XBP1 was significantly and positively correlated (p < 0:05) with the expression of JAG1 (A), NOTCH2 (B), and NOTCH3 (C) (Supplementary Figure 9A-9C). Because NOTCH2 and NOTCH3 are key effectors of the Notch signaling pathway, we proposed that Xbp1 may upregulate the Notch signaling pathway through Jag1 transcription.
GSE144306 and GSE168153 are datasets of the same study, and GSE144306 is an ATAC-seq dataset belonging to the tendon coupling region. The different transcriptional open lengths of Xbp1 and Jag1 in the cartilaginous region of tendon fibers are shown in Figure 10(a). This suggests that Xbp1 is common chromatin developed in different samples, while chromatin opening of Jag1 is significantly heterogeneous in GSE144306. A global analysis of GSE144306 was also performed, and a heat map was drawn for the analysis of chromatin opening in different samples of the fibrocartilaginous region of the tendons (Figure 10(b)). As a result  17 Oxidative Medicine and Cellular Longevity of the ATAC-seq analysis, the MEME Suite database predicted that the promoter region of Jag1 might be regulated by Xbp1 [58,70]. The promoter of Jag1 was significantly more open in the chondrocyte region (Figure 10(d)). By contrast, Xbp1 did not show this phenomenon (Figure 10(e)). This suggests that the transcriptional activity of Jag1 is significantly upregulated in the chondrocyte region, which may be due to the promoter activation of Jag1. Taken together, we revealed that Xbp1 might upregulate the Notch signaling pathway by promoting Jag1 transcription based on the ligand-receptor analysis and ATAC-seq analysis.
3.9. Crosstalk in the Notch/Hh Signaling Pathway Is Associated with Xbp1 Expression. The GSVA based on GSE94683 revealed that upregulated activities of Hh and Notch signaling pathway activities were upregulated in HO tissues (Supplementary Figure 9D, 9E). Interestingly, the activity of the Hh signaling pathway was more expressed in HO2 than in HO1, and the activity of the Hh signaling pathway was more expressed in HO1 than in HO2 (Figures 11(a)-11(c)). By contrast, the proposed temporal expression profile of the Hh and Notch signaling pathways suggests that the Hh signaling pathway may be activated before the Notch signaling pathway. The GSEA of XBP1 expression concerning the pathway activity in different datasets was also performed (Figures 11(e)-11(h)). The upregulations of Xbp1 in GSE94683, GSE102929, GSE168153, and GSE138515 were associated with the upregulation of the WNT, Notch, Hh, and TGF beta signaling pathways. This suggests that Xbp1 upregulation may induce the occurrence of ectopic 19 Oxidative Medicine and Cellular Longevity ossification by upregulating the activities of WNT, Notch, Hh, and TGF beta signaling pathways.

Discussion
In this study, we used a single-cell approach to integrate and analyze the developmental lineage of the fibrocartilaginous zone of tendons and the ectopic ossification cell model.
Ectopic ossified cells demonstrated similar transcriptional characteristics to cells in the fibrocartilaginous region of tendons. This study also identified Xbp1 as a common key TF for ectopically ossified cells and cells in the fibrocartilaginous zone of the tendons. To the best of our knowledge, this is the first study to propose that Xbp1 may upregulate the Notch signaling pathway by transcribing Jag1. These results also suggest a potential association of endoplasmic reticulum stress with tendon bone healing and HO.
This study was performed based on the integration analyses of scRNA-seq, bulk RNA-seq, and ATAC-seq data. The CCA in the Seurat package was used for dataset integration, whereby we found that ectopically ossified cells had a similar transcriptional profile to cells in the fibrocartilaginous region of the tendons. Based on the pseudo-temporal trajectory analysis of Monocle 3, this study revealed a potential model of PC differentiation in ectopically ossified tissues. To identify specific TFs associated with ectopic ossification, we used a modified version of the SCENIC method to construct a regulatory network between HO-related targets and TFs. The results showed that Xbp1 was a common key TF for ectopically ossified cells and cells in the fibrocartilaginous region of the tendons. We then performed cellular ligand-receptor analysis using the CellPhoneDB database and software according to the instruction manual, and after further pathway screening, Xbp1 was found to possibly upregulate the Notch signaling pathway by promoting Jag1 transcription. GSE94683 (including nine samples of healthy donors and seven samples of ectopically ossified tissues from patients with HO) was subjected to differential analysis. Five
We first classified pathological cells in the ectopic ossification model into HO1 and HO2 cell subpopulations. The proposed chronological analysis suggests that HO2 is a precursor cell of HO1. Mgp, Bgn, Sox9, Col11a1, Col1a1, Col1a2, Tnmd, Col11a2, Scx, Hivep2, Bcl2, Col2a1, Mkx, Myc, and Dcn were enriched in HO2 cells. Sox9 and Scx were key TFs in the migrating connective tissue between the tendon and bone [34,71]. Diabetes was found to contribute to the development of tendinopathy through the downregulation of Scx, Mkx, Col1a1, Col1a2, and Bgn [72], and the upregulation of Col1a1, Scx, and Dcn contributed to the control of tendinopathy [73]. Tnmd was expressed at high levels in tendon cells and ligament cells and was found to have great potential in the treatment of tendon repair, which might be related to its coordinated expression of Col2a1 synergistically promoting angiogenesis [74][75][76]. Mkx controlled the differentiation of tendon and ligament cells, while Bcl2 was associated with the apoptosis and regeneration of tenocytes [77,78]. Col5a1, Foxc1, Gli3, Gli2, Pecam1, Klf2, Klf4, Zbtb48, Ptch2, Smad7, Sufu, Igfbp5, Tcf, Fbn1, and Jag1 were mainly enriched in HO1 cells. Gli2/3, Tcf, Fbn1, and Col5a1 were molecular markers in cartilage tissues [79][80][81]. The upregulation of Klf2 and Klf4 was found to promote cartilage stability and resist osteoarticular oxidative stress conditions [82,83]. Foxc1, Ptch2, and Sufu were found to be essential for the activation of the Indian Hh signaling pathway [84,85]. Smad7 and Jag1 were found to affect cartilage development through the MAPK or Notch signaling pathway [86,87]. Igfb5 expression was found to be highly concentrated in developing limb muscle tissue and was highly correlated with osteogenesis [88]. The expression profiles of HO1 and HO2 cells were more similar to that of chondrocytes and tenocyte-like cells, respectively. Some information also implied the developmental potential of HO2 cells to become HO1 cells. For example, (1) Sox9, Col1a1, and Col2a1 have been shown to play important functions in cartilage development [89,90]; (2) Mgp was found to be related to the onset of osteoarthritis during cartilage development [91]; (3) Bcl2 was a key protein of the PI3K-AKT signaling pathway, which was closely associated with chondrocyte proliferation [92]; and (4) COL5A1 polymorphism was associated with susceptibility to tendon disease [93]. Notably, Gli3, Cli2, Pecam1, and Ptch2 were found to be the star genes of the Hh signaling pathway, and Notch2 and Jag1 are the star genes of the Notch signaling pathway, suggesting that the activity of the Notch and Hh signaling pathways may be upregulated in HO1 cells [94][95][96][97].
In this study, we performed an integrative analysis of single-cell transcriptome data of mouse Achilles tendon, single-cell transcriptome data of mouse tendon junction region, and single-cell transcriptome data of mouse tail tendon cells under simulated HO pathological conditions (including GSE168153, GSE138515, and GSE102929 datasets under simulated HO pathological conditions) was also performed in this study. In Figure 6(b), tenocyte-like attachment cells (Matn4 and Col2a1 are cell markers of this cell type) were found to have more HO2 cells, while more HO1 cells were found in chondrocyte-like attachment cell 1 (Anxa2 and Col14a1 are cell markers of this cell type). Matn4 and Col2a1 are cell markers of tenocyte-like attachment cells. Col2a1 is a crossover gene located between tenocyte and chondrocyte differentiation [76,90,98,99]. While Matn4 expression was tightly regulated and served as a

22
Oxidative Medicine and Cellular Longevity marker of cell differentiation, it played a key role in maintaining the stability of articular cartilage [100]. Anxa2 and Col14a1 are cell markers of chondrocyte-like attachment cell 1. The upregulation of Anxa2, a biomarker of cartilage development, may promote tendon healing [101,102]. The results showed that Col14a1 was tendon-associated collagen and was associated with the development of osteoarthritis [103]. These data suggest that ectopically ossified cells have similar transcriptional and developmental characteristics to cells in the fibrocartilaginous zone.
In this study, we first identified the endoplasmic reticulum (ER) stress effector Xbp1 as a common TF for HO bulk-RNA-seq, HO scRNA-seq, and the integrated scRNA data (including GSE168153, GSE138515, and GSE102929 datasets under simulated HO pathological conditions) using a comprehensive network analysis combined with a modified SCENIC approach. XBP1 was found to be an important regulatory pathway for the unfolding protein response that occurred in response to ER stress and played an important role in both physiological and pathological conditions, and its activity had a significant effect on the progression and prognosis of numerous diseases [104]. In this study, 160 genes were identified as core downstream regulatory targets of XBP1, with significant interactions among them, and they mainly played important roles in the Golgi vesicle transport, ER protein processing, cellular communication, and translocation. Lei et al. found that mechanical stress can activate the ER stress and MAPK signaling pathways in the posterior longitudinal ligament, which played an important role in the onset of ossification of the posterior longitudinal ligament [105,106]. The effect of ER stress on cartilage differentiation is regulated by the microenvironment [107]. Recent studies have shown that ER stress regulates the Notch and Hh signaling pathways [108,109]. As a key effector protein of ER stress, Xbp1 is positively correlated with the activity of the Notch signaling pathway. The upregulation of Xbp1 was found to inhibit the Notch pathway effector protein pecanex, resulting in ER enlargement [110][111][112]. In addition, the deletion of XBP1 signaling was found to cause chondrodysplasia and delayed ossification [113]. In the present study, Xbp1, a key effector protein of ER stress, was found to be a key TF for the onset of ectopic ossification. These results also suggest a potential association between endoplasmic reticulum stress and tendon bone healing and HO.
The prevalence of HO in stroke is 0.5%-1.2%, which has led to a lack of understanding of the epidemiology and pathogenesis of HO after ischemic stroke [117][118][119]. Furthermore, the occurrence of HO is dependent on the upregulation of endoplasmic reticulum stress pathways [105,106,[120][121][122]. Meanwhile, the activation of endoplasmic reticulum stress-related pathways and the upregulation of XBP1 are key features of the widespread presence of various brain-related diseases, including stroke [123]. Therefore, the coupregulation of endoplasmic reticulum stress-related pathways in HO and ischemic stroke has attracted our attention. The current study found that the expression of microRNAs carried by exosomes in the blood is altered after ischemic stroke, which may affect the function of organs and tissues in the distant compartment [124][125][126]. We propose the hypothesis that changes in microRNA levels in the peripheral blood exosomes following ischemic stroke may lead to upregulation of stem cell XBP1 expression in muscle tissue and may thus induce neurogenic ectopic ossification. In this study, 24 microRNAs were screened and were found to be potentially associated with the upregulation of XBP1 expression after acute ischemic stroke. However, the effect of circulating blood exosomes on the occurrence of HO remains to be further investigated.
In this study, cell clusters of ectopic ossification phenotypes were isolated in a single-cell dataset of ectopic ossification. Three single-cell transcriptomes (namely, fibrocartilaginous zone single cell, Achilles tendon single-cell transcriptome, and ectopic ossification model transcriptome) were integrated into one dataset using the CCA method. Ectopically ossified cells were found to colocalize with the fibrocartilaginous zone cells, suggesting that ectopically ossified cells share similar transcriptome characteristics to fibrocartilaginous zone cells. The pseudo-temporal trajectory analysis further supported the above findings. Further TF coexpression network analysis revealed that Xbp1 is a hub TF for HO formation. The ligand-receptor analysis reveals potential pathways through which Xbp1 exerts its transcriptional regulatory role. The single-cell analysis revealed that the Xbp1 and Notch pathway expressions were upregulated in ectopic ossified tissues, and they were positively correlated with each other. Ultimately, we proposed a potential mechanism by which Xbp1 may upregulate the Notch signaling pathway by promoting Jag1 transcription. In summary, this study integrated six datasets, namely, GSE94683, GSE144306, GSE168153, GSE138515, GSE102929, and GSE110993. The differentiation trajectory and key TFs of ectopic ossification occurrence were systematically analyzed by integrating scRNA-seq, bulk RNA-seq, and ATAC-seq data of the fibrocartilaginous region and HO.
This study revealed similarities between ectopically ossified cells and cells of the fibrocartilaginous zone of tendons and suggested a potential Xbp1/Jag1/Notch signaling pathway. However, a theoretical work with data from databases with no experimental verification may limit the application value of this research. Further functional experiments are still needed to be performed to explore the specific functional mechanisms of Xbp1. In this study, we obtained data from an integrated analysis of multiple databases and planned to conduct further analysis using ectopic ossified tissues and tendon fibrocartilaginous zone cells from the same model. In addition, the regulatory relationships between the Hh and Notch signaling pathways and other key pathways need to be further validated to confirm the activation patterns of the pathways in the single-cell transcriptome. We believe that further collection of clinical samples is also necessary, which will help in clinical diagnosis and translation.

Conclusion
A systematic analysis of the differentiation landscape and cellular homogeneity facilitates a molecular understanding of the potential homology between tendon fibrocartilage and ectopic ossified tissue. Furthermore, by identifying Xbp1 as a hub regulator and through ligand-receptor analysis, we propose a potential Xbp1/Jag1/Notch signaling pathway. The findings provide a novel concept for the formation of ectopically ossified tissues, i.e., HO cells may be differentiated from misplaced fibrocartilaginous zone stem cells.

Data Availability
The corresponding author agreed to provide the reader with the original data for a suitable reason.

Conflicts of Interest
The authors declare no conflict of interest.