Evolutionary and developmental characteristics of IBs
Based on the information provided by Patterson and Johnson [14], Yang et al. [25] and Nelson [26], we summarized 5 IB distribution patterns in teleost (Fig. 1A). In Clupeiformes, there are 2 patterns including four-categories-IB model (epineural bones-ENBs, epipleural bones-EPBs, epicentral bones-ECBs and myorhabdois-MBs) and three-categories-IB model (ENBs, EPBs and ECBs); There is no IBs in Siluriformes; Cypriniform fishes have 2 categories of IBs, ENBs and EPBs; Salmoniform fishes have only ENBs and Perciform fishes have only ECBs. IBs are diverse in distribution and found only in teleost fishes [14], which may be attributed to teleost-specific whole-genome duplication (Ts3R).
Cypriniformes constitute the largest monophyletic group of freshwater fishes on the earth [25], hence the IB distribution pattern in Cypriniformes is the most widespread. Furthermore, we focus on the developmental process of IBs in Cypriniformes. As a representative species of Cypriniformes, Wuchang bream was used to observe the developmental process of IBs. We explored the histological structure of IB at different developmental stages (Fig. 1B-D, Additional file 2: Table S1) (stage 1: Only unossified tendon distribution without distribution of IB; stage 2: Both immature IBs and unossified tendons distribution; stage 3: Both mature IBs and unossified tendons distribution) by alizarin red staining. The results are consistent with Chen et al. [6, 18].
Gene expression characteristics of IBs in teleosts
The evolutionary origin of any new organ typically depends on the recruitment of genes that were originally expressed in other tissues [27]. IB developed later relative to other tissues and can be considered as a novel tissue. To identify genes recruited in ENBs and EPBs, we compared 13 transcriptomic data from Wuchang bream obtained in this study (Additional file 2: Table S2). In addition, to assess the expression signature of tissue-specific genes at three different developmental stages of IBs, 6 transcriptomic data were used from Liu et al. [28]. Genes specifically expressed in ENBs or EPBs (hereafter ENB-specific genes or EPB-specific genes) were defined as those that have a Tau index exceeding 0.8 [29] and are expressed most strongly or second most strongly in ENBs or EPBs. 130 ENB-specific genes (Fig. 2A, Additional file 2: Table S3) and 224 EPB-specific genes (Fig. 2B, Additional file 2: Table S4) were identified. Main sources of ENB were EPB (70.00%), neural spines (6.92%), operculum (6.92%), suboperculum (3.85%) and vertebra (3.85%) (Fig. 2A); Main sources of EPB were ENB (40.62%), skin (14.73%), rib (10.27%), operculum (9.38%) and neural spines (6.70%) (Fig. 2B). Importantly, 91 intermuscular bone-core genes were shared by both ENBs and EPBs and enriched in biological processes related to IB development (Fig. 2C, Additional file 2: Table S5, Additional file 2: Table S6). To elucidate the temporal expression characteristics of core genes, time-series analysis was performed. The result of time series analysis showed that high expression of most core genes at early stages of IB formation (stage1 and stage2) which suggests that the IB formation was designed during the early development of fish (Additional file 1: Fig.S1). For the purpose of examining the impact of Ts3R on IB development, we screened for genes undergoing Ts3R in the Wuchang bream genome and found 8 core genes (dbx1b, fgf8b, LOC125252216, LOC125264141, LOC125278619, plg, slc16a1a, tfa) undergoing Ts3R (Fig. 2D, Additional file 2: Table S7). These genes are all highly expressed in the early stages of IB development. Plg was enriched in blood circulation (GO:0008015 pvalue = 8.418093e-04), extracellular matrix organization (GO:0030198 pvalue = 2.142292e-03). Tfa was enriched in positive regulation of bone remodeling (GO:0046852 pvalue = 6.775198e-04), positive regulation of MAPK cascade (GO:0043410 pvalue = 2.380401e-03) and regulation of bone mineralization (GO:0030500 pvalue = 7.032019e-03) (Fig. 2D). In summary, ENB-specific genes and EPB-specific genes were recruited from similar tissues and shared 91 same genes, such results suggested that they share similar gene expression characteristics and tissue origins. Subsequent analysis only focused on the genes shared by ENBs and EPBs. Meanwhile, Ts3R may had an impact on the occurrence of IBs in teleost.
Identification of genes related to tendon ossification across species
Primitive teleosts with IBs contain only ENB, then they develop EPBs [14]. However, ENBs and EPBs occur in Cypriniformes but are absent in Siluriformes. In Cypriniformes, Nie et al.[10] found that stem cells in myosepta could differentiate into osteoblasts and tenocytes and revealed the role of osteoblasts in IB formation. Therefore, we believed that myoseptal stem cell of IB-containing fish had two differentiation directions: from stem cell to osteoblast and from stem cell to tenocyte (Fig. 4B). Myoseptal stem cell of IB-lacking fish can only differentiate into tenocyte (Fig. 4B). With a goal of identifying the genes responsible for the formation of IB in the evolution of teleost, we conducted a comparative transcriptomic analysis of myoseptal tissues in fish with and without IBs. IB-containing group included Wuchang bream (M. amblycephala) and grass carp (Ctenopharyngodon idella); IB-lacking group included yellow catfish (Tachysurus fulvidraco) and channel catfish (Ictalurus punctatus) (Fig. 4A).
To conduct cross-species differential expressed analysis, 19,997 one-to-one orthologs in IB-containing group (Additional file 2: Table S8) and 15,818 one-to-one orthologs in four species (Additional file 2: Table S9) were identified. Furthermore, 3,302 specific genes in IB-containing group were identified (Fig. 3A, Additional file 2: Table S10). We constructed comparative groups (ENB vs. epineural muscles (ENMs include muscles and epineural tendons), EPB vs. epipleural muscles (EPMs include muscles and epipleural tendons) within IB-containing group. After performing differential expressed analysis using GFOLD [30], we obtained a intersection of 623 IB-core genes (differential expressed genes, DEGs) with the same expression trend (Fig. 3B, Additional file 2: Table S11). They were enriched in upregulated biological processes (extracellular matrix organization (GO:0030198), ossification (GO:0001503), regulation of angiogenesis (GO:0045765), connective tissue development (GO:0061448) and others) and downregulated biological processes (protein ubiquitination (GO:0016567) and others) (Fig. 3B, Additional file 2: Table S12). The GSEA (Gene set enrichment analysis) results indicated that upregulated biological processes may be associated with promoting IB development, while downregulated biological processes may be associated with inhibiting IB development. We performed similar differential expressed analysis using GFOLD and SCBN [31] between IB-containing group and IB-lacking group (details in Methods). Finally, 209 IB-core genes (DEGs) with the same expression trend were identified in IB-containing vs. IB-lacking group (Fig. 3C, Additional file 2: Table S14). Divide the two gene sets (623 IB-core genes and 209 IB-core genes) into 4 categories, including 69 specific genes in IB-containing fish, 79 intersection genes shared by two gene sets, 475 other genes in 623 core DEGs and 130 other genes in 209 core DEGs, based on homology and gene expression trends (Fig. 3C). In IB-containing fish group, 11.08% of 623 IB-core genes were IB-containing-fish specific genes; 12.68% of 623 IB-core genes were shared by IB-containing group and IB-containing vs. IB-lacking group; 76.24% of 623 IB-core genes were genes that show no significant differences in IB-containing vs. IB-lacking group or did not belong to one-to-one orthologous gene pairs. In IB-containing vs. IB-lacking group, 37.80% of 209 IB-core genes were shared by IB-containing group and IB-containing vs. IB-lacking group; 62.20% of 209 IB-core genes were genes that show no significant differences in IB-containing fish group (Fig. 3C). These four categories of gene sets exhibited two different expression trends during the formation of IBs in three developmental stages, with one being highly expressed in early IB development (Cluster2 and Cluster3) and the other being highly expressed in late IB development (Cluster1) (Fig. 3D).
28 upregulated and 41 downregulated specific genes in IB-containing fish were shown in Fig. 3E. Clec3bb, cltcl1, six4a, fgfr1b, fgf13a, chrna1 and LOC125277479 were enriched in up-regulated biological processes functioning to IB formation (Fig. 3F, Additional file 2: Table S13), and amfrb was down-regulated specific gene in IB-containing fish which involved in protein ubiquitination (Additional file 2: Table S13). 37 upregulated and 37 downregulated intersection genes in IB-containing vs. IB-lacking group were also shown in Fig. 3E (tnnt3b, rtn4a, fxr1, zgc:91999 and pfkfb3 were down-regulated in IB-containing fish group but up-regulated in IB-containing vs. IB-lacking group). Loxa, timp2a, col5a1, lum, thbs2a, mmp2, ryr1b, aspn, tnfrsf11a and mef2aa involved in up-regulated biological processes functioning to IB formation (Fig. 3F, Additional file 2: Table S13). In short, there were significant changes in gene expression related to IB formation during the evolution of teleost. Only a fraction of core genes exhibited the same expression trends, and these genes may be associated with the process which was from myoseptal stem cell to osteoblast.
Candidate genes show altered expression in IB-lacking zebrafish
Based on comparative analysis of multiple species and tissues, we identified 79 differentially expressed genes that are conserved across species, 69 differentially expressed genes exclusive to IB-containing fish, and 91 IB-specific genes calculated by Tau analysis (Fig. 4A). In order to screen for genes associated with tendon ossification (Fig. 4B), we further investigated the expression of the above candidate genes in runx2b−/− zebrafish (IB-lacking zebrafish by genome editing, Nie et al. [10]). Upon reanalysis of the transcriptome data (runx2b−/− and runx2b+/+ muscle tissue from tail tissue) from Xiao et al. [11], a total of 1,840 differentially expressed genes were identified (Fig. 4A, Additional file 1: Fig.S3, Additional file 2: Table S15). It was found that 4 IB-containing fish specific DEGs (clec3bb, myl2b, actn3b, and ptx3b) (Fig. 4C-4E), 4 conserved intersection DEGs (krt18b, si:dkey − 282h22.5, si:ch73 − 86n18.1 and si:ch211 − 243a20.3) (Fig. 4F-4H) and 2 IB-specific genes (prokr1a and cp) (Fig. 4I-K) exhibit opposite expression patterns in the group of genome editing zebrafish through intersection analysis.
Clec3bb shows specific and elevated expression in tenocytes of zebrafish without IBs
In order to elucidate the single-cell expression atlas of genes associated with IB formation, we reanalyzed the single-cell transcriptome data from runx2b+/+ and runx2b−/− zebrafish tail tissues [10]. The uniform manifold approximation and projection (UMAP) plotting of the scRNA-seq data identified 14 distinct cell clusters (Additional file 1: Fig.S4, Additional file 2: Table S16). To reveal the potential differentiation relationship between tenocyte/fibroblast cell and osteoblast, we investigated the diversity tenocytes, fibroblasts and osteoblasts in more detail (Additional file 1: Fig.S4, Additional file 2: Table S16). Sub-clustering analysis was performed, which revealed an unexpectedly large degree of diversity, with 10 transcriptionally distinct clusters (Fig. 5A, Additional file 2: Table S17). There were significant differences in the proportion of most clusters between runx2b+/+ and runx2b−/− zebrafish. The proportions of clec3bb + tenocyte, fibroblast2 and osteoblast1/2 were higher in runx2b+/+ zebrafish than runx2b−/− zebrafish. Conversely, the proportions of tenocyte2, fibroblast2 and mstnb + tenocyte were higher in runx2b−/− zebrafish than runx2b+/+ zebrafish. There was no significant difference between runx2b+/+ and runx2b−/− zebrafish in tenocyte1, myoseptal stem cell and thbs4b + tenocyte cluster (Fig. 5B). The tenocyte supercluster consisted of four distinct subpopulations, including differentiating tenocyte clusters (tenocyte1 and tenocyte2, which were marked by scxa), clec3bb + tenocyte, thbs4b + tenocyte and mstnb + tenocyte (Fig. 5B). Clec3b gene is associated with bone mineralization, and its deficiency leads to skull defects in mice [32]. Thbs4b controls matrix assembly during development and repair of myotendinous junctions in zebrafish [33]. Spontaneous and targeted mstn loss-of-function mutations lead to double muscle phenotypes in zebrafish [34]. Mstnb + tenocyte may be located in the myotendinous junction, serving as the boundary between muscle and tendon. There were 3 subclusters in fibroblast supercluster, including fibroblast1 (nt5e), myoseptal stem cell (cxcl12a) and fibroblast2 (fn1a) (Fig. 5B). Cxcl12, a marker commonly used for hematopoietic stem cells [35], osteoblasts [36] and mesenchymal stem cells of HO [9]. Cells in myoseptal stem cell cluster may be stem cells in myosepta contributing to tendon ossification. Osteoblast supercluster included osteoblast1 and osteoblast2, which was marked by runx2b. Cd81b, a marker of exosomes [37], was specific expressed in osteoblast2. The result suggested that osteoblasts2 interact with other cells via exosomes.
With the hope of uncovering the effect of runx2b deficiency on cell differentiation, RNA velocity analysis was performed in runx2b+/+ and runx2b−/− single-cell data. Myoseptal stem cell could differentiate into osteoblast, tenocyte, and fibroblast in runx2b+/+ zebrafish (Fig. 5C). Due to runx2b deficiency, cell differentiation trajectory was affected in runx2b−/− zebrafish. The ability of myoseptal stem cell to differentiate into osteoblast1/2 was weakened, and the ability to differentiate into fibroblast1/2 and clec3bb + tenocyte was enhanced (Fig. 5D).
The expressions of candidate genes which were from the previous section of the results were also altered in runx2b−/− zebrafish, and the results were consistent with the bulk RNA-seq (Fig. 4). Clec3bb was involved in ossification (GO:0001503) and was up-regulated in clec3bb + tenocyte of runx2b−/− zebrafish. Fgfr1b, involving in connective tissue development (GO:0061448), was up-regulated in osteoblast1 and osteoblast2 of runx2b−/− zebrafish. Cnn2 was up-regulated in all clusters except osteoblast2 in runx2b−/− zebrafish. Col2a1b was down-regulated in all clusters in runx2b−/− zebrafish. Lum participated in extracellular matrix organization (GO:0030198) and connective tissue development, which was down-regulated in myoseptal stem cell, and osteoblast supercluster of runx2b−/− zebrafish. Thbs4b participated in connective tissue development, which was down-regulated in fibroblast2 and osteoblast1 of runx2b−/− zebrafish. Aspn participated in ossification, which was down-regulated in fibroblast2 and osteoblast2 of runx2b−/− zebrafish. Loxa, involving in extracellular matrix organization, was up-regulated in all clusters except fibroblast1, fibroblast2 and osteoblast2 of runx2b−/− zebrafish.
IB-specific genetic changes and gene expression analysis in tetranectin family
Clec3b belongs to tetranectin gene family [38]. This gene family includes: clec3ba, clec3bb, clec3a and clec11a. To explore IB-specific genetic changes of the tetranectin gene family in teleosts, we compared conserved motifs, gene structures and protein structures from four IB-containing fishes and three IB-lacking fishes. Osteolectin (also known as clec11a) is an osteogenic factor that promotes the maintenance of the adult skeleton by promoting the differentiation of LepR + cells into osteoblasts [39]. Clec11a-deficient mice exhibits accelerated bone loss during aging, reduced bone strength, and delayed fracture healing [40]. Motif10 was specific to IB-lacking fish (Additional file 1: Fig.S5). CLEC3A is identified as cartilage-specific antimicrobial peptides in septic arthritis, and contributes to cell proliferation through HIF1α signaling pathway [41]. Motif7 was specific to IB-lacking fish (Additional file 1: Fig.S5). C-Type Lectin Domain Family 3 Member B (CLEC3B), known as tetranectin (TN), is detected in the cytoplasm, extracellular matrix (ECM) and exosomes [42]. Teleosts had clec3ba and clec3bb under the action of Ts3R, but clec3bb was lost in IB-lacking fish. There was no difference in CLEC3BA conserved motifs of IB-lacking fish and IB-containing fish. Otherwise, clec3ba genes of IB-containing fish had longer UTR sequence than IB-lacking fish’s (Additional file 1: Fig.S5). We further compared the differences in amino acid sequences between CLEC3BA and CLEC3BB in IB-lacking fish and IB-containing fish (Fig. 6A). There were four IB-specific amino acid changes (K26T, N35H, Q60H and H158R) and a missing fragment of IB-containing fish in CLEC3BA. The presence of conserved amino acids (such as A61, K84 and W152) also implied the conserved nature of gene function between CLEC3BA and CLEC3BB (Fig. 6A).
Furthermore, Clec3bb, as an IB-containing fish specific gene, was expressed more highly in ENBs and EPBs than in any other tissues (Fig. 6B) and was expressed in myoseptas of zebrafish embryo (Fig. 6C) and in ossified tendons (IBs) of adult zebrafish (Fig. 6D).
Clec3bb promotes osteoblast2 mineralization in myosepta of IB-containing fish
C-Type Lectin Domain Family 3 Member B (CLEC3B) is a transmembrane Ca2+-binding protein, locating in cell plasma, extracellular matrix and exosomes [43]. In runx2b+/+ zebrafish, clec3bb was specifically expressed in clec3bb + tenocytes, and its expression products functioned to osteoblast2 through exosomes to promote osteoblast2 mineralization (Fig. 7A). Osteoblast2 specifically expressed exosome marker cd81b (Fig. 5B), suggesting that osteoblast2 may interact with exosomes. In runx2b−/− zebrafish, runx2b deficiency inhibited the formation of osteoblast1 and osteoblast, which also led to the loss of exosomal-clec3bb-targeted cells, which further led to redundant expression of clec3bb gene (Fig. 7B). Slingshot analysis revealed that runx2b deficiency resulted in differentiation from thbs4b + tenocytes to clec3bb + tenocytes (Fig. 7A-B).