Effect Of XBP1 Deficiency In Cartilage On The Regulatory Network Of LncRNA/circRNA-miRNA-mRNA

X-box binding protein 1(XBP1) is a critical component for unfolded protein response (UPR) in ER stress. According to previous studies performed with different XBP1-deficient mice, the XBP1 gene affects mouse cartilage development and causes other related diseases. However, how the complete transcriptome, including mRNA and ncRNAs, affects the function of cartilage and other tissues when XBP1 is deficient in chondrocytes is unclear. In this study, we aimed to screen the differentially expressed (DE) mRNAs, circRNAs, lncRNAs and miRNAs in XBP1 cartilage-specific knockout (CKO) mice using high throughput sequencing and construct the circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA regulatory networks. DE LncRNAs (DE-LncRNAs), circRNAs (DE-circRNAs), miRNAs (DE-miRNAs), and mRNAs [differentially expressed genes (DEGs)] between the cartilage tissue of XBP1 CKO mice and controls were identified, including 441 DE-LncRNAs, 15 DE-circRNAs, 6 DE-miRNAs, and 477 DEGs. Further, 253,235 lncRNA-miRNA-mRNA networks and 1,822 circRNA-miRNA-mRNA networks were constructed based on the correlation between lncRNAs/circRNAs, miRNAs, mRNAs. The whole transcriptome analysis revealed that XBP1 deficiency in cartilage affects the function of cartilage and other different tissues, as well as associated diseases. Overall, our findings may provide potential biomarkers and mechanisms for the diagnosis and treatment of cartilage and other related diseases.


Introduction
Abnormal chondrocyte differentiation is an important cause of common orthopedic diseases, such as rickets, dwarfism, and osteoarthritis [1,2]. During endoplasmic reticulum stress, IRE1α undergoes autophosphorylation and activation, and the XBP1u mRNA is spliced in the cytoplasm to produce the transcription factor, XBP1s, with multiple regulatory functions. This transcription factor enters the nucleus to regulate the transcription and expression of its associated gene, and participates in the differentiation and development of chondrocytes [3][4][5][6].
The unfolded protein response (UPR) signaling pathway, inositol-requiring enzyme 1 (IRE1)-XBP1, which is mediated by XBP1s, is an important pathway in the determination of cell fate and a key pathway for plasma cell differentiation [3][4][5]. According to our previous studies, XBP1s is expressed in the proliferation and hypertrophy stages of chondrocyte differentiation and is involved in the regulation of chondrocyte differentiation and cartilage development in chondrocytes. As a BMP2-inducible transcription factor, XBP1s positively regulates endochondral bone formation by activating GEP chondrogenic growth factor [6]. Cameron et al. reported that the cartilage-specific deficiency of XBP1 in mouse leads to a chondrodysplasia characteristically defined by reduced chondrocyte proliferation and delayed cartilage maturation and mineralization. XBP1 was inferred to control chondrocyte proliferation and the timing of cartilage maturation and matrix mineralization during endochondral ossification [7,8]. As reported in many literatures, in eukaryotic genomes that transcribe up to 90% of genomic DNA, only 1-2% of the transcripts encode proteins, with most being transcribed as non-coding RNAs (ncRNAs). The number of ncRNAs increases significantly with the complexity of the organism. Although the expression level of most predicted ncRNAs is markedly lower than that of mRNAs, they do not affect their regulatory functions. Based on increasing evidence, ncRNA has a regulatory role in the development process and in the response to stress and environmental stimuli [9][10][11][12][13][14].
Non-coding RNAs play an important regulatory role in cells. In fact, these RNAs, can regulate several biological processes, including epigenetic modification, variable RNA splicing, and protein stability [15][16][17]. For instance, rRNAs, tRNAs, small nuclear RNA (snRNAs), small nucleolar RNAs (snoRNAs) are involved in mRNA translation, splicing, modification of rRNAs, respectively [10]. High-throughput sequencing (NGS) technology and bioinformatics analysis have led to the discovery of a remarkable number of circRNAs, lncRNAs, miRNAs, mRNAs and their networks [18]. Long-noncoding RNAs (LncRNAs) are ncRNAs longer than 200 nt [19]. Previous studies have shown that lncRNAs play a crucial role in the regulation of disease processes and gene expression [20,21]. Further, lncRNAs are reported to participate in cell differentiation, migration, proliferation and apoptosis [22]. For example, a novel osteogenesis-associated lncRNA (lncRNA-OG) significantly promotes BM-MSC osteogenesis. LncRNA-OG interacts with heterogeneous nuclear ribonucleoprotein K (hnRNPK) protein to regulate the activation of the bone morphogenetic protein signaling pathway. hnRNPK positively regulates the transcriptional activity of lncRNA-OG by promoting H3K27 acetylation of the lncRNA-OG promoter [23]. CircRNAs is a class of ncRNAs widely expressed in eukaryotic cells. CircRNAs have a covalent closed circular configuration without 5' and 3' polarities and poly A tail structure [24]. CircRNA can act as an miRNA sponge to affect the functions of miRNAs [25], and has a regulatory effect on several diseases, such as cancer, diabetes, cardiovascular diseases and certain immune diseases [26].
MiR-124-3p is an essential microRNA (miRNA), and its expression changes are related to the proliferation ability of bone mesenchymal stem cells (BMSCs) [27]. MiR-124-3p is also reported to be associated with osteoporosis (OP) and fragility fractures [28]. RNAs have become a novel hotspot in research in cartilage development. Although some ncRNAs, including lncRNAs, circRNAs and miRNAs, have been identified in bone development and associated diseases, the ncRNAs related to the regulation of XBP1 have not been revealed.
The coordinated regulation between different tissues and organs is performed by the full spectrum of transcriptome RNAs, including mRNAs and ncRNAs. In the current study, we attempted to elucidate whether and how XBP1 deficiency affects the full spectrum of ncRNA functional characterization and the normal physiological functions of ncRNA and their role in human diseases. First, we conducted whole transcriptome sequencing of cartilage tissue obtained from XBP1 conditional knockout mice, and carried out a differential expression profile assay of circRNAs, lncRNAs, miRNAs and mRNAs. Thereafter, based on bioinformatics analysis, two competitive endogenous RNA (ceRNA) regulatory networks (the circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA) were established, and GO and KEGG enrichment analyses were performed. XBP1 deficiency in cartilage leads to differences in the complete transcriptome RNA of cartilage and other different tissues, including various ncRNAs with differential expression and regulatory profiles, as well as differences in protein expression and abnormalities of various different secreted factors. Our results support the identification of both the full spectrum of the function and characterization of ncRNAs induced by XBP1 deficiency in cartilage and the role of these XBP1-related ncRNAs in normal physiological functions and human diseases.

Total RNA Extraction and Quality Control
Total RNA from 6 samples was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The concentration of total RNA was checked by using Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer, its purity was checked by using a NanoPhotometer spectrophotometer (IMPLEN, CA, USA), and its integrity was checked by Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

RNA-seq Library Preparation and Sequencing
For constructing the library of mRNA or circRNA, 5 μg of RNA for per sample was used as input. For mRNA library construction, Epicentre Ribo-Zero™ rRNA Removal Kit (Illumina, USA) was used to remove rRNA from the total RNA. For circRNA library construction, 3 U of RNase R (Epicentre, USA) per microgram was added to the rRNA removal system and incubated to remove linear RNA. Subsequently, the sequencing libraries were prepared in accordance with by using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA). The Agilent Bioanalyzer 2100 system was used to detect the quality of the RNA library. After the library was qualified, the Illumina PE150 sequencing was performed.
Considering that the Ribo-Zero library contains both mRNA and lncRNA, we will select lncRNA by establishing a series of strict screening conditions in the data analysis part. For miRNA library construction, 5 μg of total RNA for per sample was used as input using the NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA) according to the manufacturer's protocol. Similarly, library quality was detected on the Agilent Bioanalyzer 2100 system. the Illumina PE150 sequencing was performed.

Data Analysis
The libraries were sequenced using an Illumina sequencing platform. Then adapter sequences and reads with poor quality were filtered using FASTX toolkit pipeline (version 0.0.13). Clean reads were aligned to the UCSC (mus_musculus_Ensembl_97) using Hisat2 (http://ccb.jhu.edu/software/hisat2) [10]. For transcriptome assembly, the mapped reads of each sample were assembled by StringTie (v1.3.1) [32], which has specific parameters for different libraries, and can accurately splice transcripts then achieve transcript quantification [33]. Then, GC contents, Q20 and Q30 of the clean data were calculated.
We use soft of find_circ [34] and CIRI2 [35] to detect and identify circRNA.The workflow of find_circ is as follows: 1) The reads that are continuously aligned with the genome are filtered out, and the spliced reads are retain; 2) the end of each candidate read are mapped to the genome to find unique anchor position; 3)candidate circRNAs are confirmed when the 3'end of the anchor sequence is aligned to the upstream of the 5'end of the anchor sequence, and the presumed breakpoint is flanked by GU/AG splicing signals. CIRI was first analyzed the CIGAR value in the sam file, scans the PCC signal (paired chiastic clipping signals) from the sam file, then filters the Junction reads based on the PEM and GT-AG signals, and finally detects unbalanced junction reads based on the DM algorithm. Filter and prevent false positives caused by homologous gene similarity and repeated sequences. We will choose two software for joint analysis to improve the accuracy of circRNA identification.
Small RNA tags were mapped onto the assembled genome of mus_musculus_Ensembl_97 using Bowtie [36,37]. The mapped small RNA tags were screened in the miRBase 20.0 database to search for known miRNA. In addition, modified software mirdeep2 and srna-toolscli were used to obtain the potential miRNA and draw the secondary structures [38]. The hairpin structure of miRNA precursors can be used to predict novel miRNA. We integrate miREvo [39] and mirdeep2 [38] to predict the novel miRNAs. Meanwhile, custom scripts were used to obtain the identified miRNA counts and base bias on the first position with certain length and on each position of all identified miRNA, respectively [40]. miRNA expression was estimated using TPM (transcript per million) Density distribution [41]. Differentially expressed miRNAs were identified using the threshold: P-value < 0.05.

Differential expression gene analysis (DEG analysis) and Enrichment Analysis
The gene expression level of RNA-seq is usually measured with FPKM(Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) [42]. We used a cutoff value of FPKM > 1 to define the gene expression. Differential expression gene analysis was performed using edgeR [43]. We use p-value or corrected p value (padj) to determine the significance level. The significant difference standard is padj less than 0.05. Differentially expressed genes are screened and subjected to hierarchical cluster analysis. GO (Gene Ontology) is a comprehensive database describing gene functions, which can be divided into three parts: molecular function, biological process, and cellular component. GO analysis was performed by GOseq [44] (Bioconductor, USA). The p-value of GO term and the false discovery rate (FDR) of p-value (q-value) were calculated to find out the most relevant GO term for the genes with differential expression. GO enrichment takes padj less than 0.05 (default) or pvalue less than 0.05 as significant enrichment. Pathway enrichment analysis was performed using KEGG database [45] (Kanehisa Laboratories, Japan). Hypergeometric test was applied to identify the most significant enriched pathways among the candidate differently expressed genes. Similarly, for KEGG pathway enrichment, pvalue or padj is less than 0.05 as significant enrichment.
The lncRNA-miRNA-mRNA ceRNA regulation network was based on the theory that lncRNAs can directly interact by invoking miRNA sponges to regulate mRNA activity [46]. "GDCRNA Tools" (http://bioconductor.org/packages/devel/bioc/html/ GDCRNATools.html) package in R software were used [47]. We use Cytoscape 3.6.1 software to plot lncRNA-miRNA-mRNA network (Biological network exploration with Cytoscape 3). The candidate miRNA-circRNA regulation network were predicted using miRanda (MicroRNA targets in Drosophila). Similarly, the target mRNAs of DE-miRNAs were predicted by scanning for conserved miRNA target sites with MiRanda(MicroRNA targets in Drosophila), and a miRNA-mRNA regulation network was constructed. We use Cytoscape 3.6.1 software to combine circRNA-miRNA network and miRNA-mRNA network to generate circRNA-miRNA-mRNA network (Biological network exploration with Cytoscape 3).

Validation of the Differential Expression of lncRNA and circRNA
To order to verify the results of transcriptome sequencing, we performed real-time PCR to determine the RNA levels of five lncRNAs and three circRNAs in chondrocyte and cartilage tissue. Quantitative real-time PCR amplification was performed by the ABI 7500 Real-Time PCR Systems (Applied Biosystems, USA). The first-strand cDNA was synthesized with PrimeScript™ RT Master Mix (Takara, Japan) and 2 μg total RNA. The amplification procedure was one cycle of 30s at 95 °C for pre-degeneration followed by 40 cycles of 30s at 95 °C and 34 s at 60 °C. The Ct values of the reference gene and the target genes were obtained according to the amplification curve. The relative quantification of gene expression was obtained by using the 2 −ΔΔCt method. GAPDH was used to normalize mRNA levels. The primer sequences used for real-time PCR are shown in Table1

Identification of DE mRNAs between XBP1 CKO mice and control mice
To identify the expression levels of mRNAs in XBP1 CKO mice and control mice, 6 mRNA libraries were constructed. Compared with the wildtype control, 477 genes, including 307 (74.2%) upregulated genes and 107 (25.8%) downregulated genes, were differentially expressed in XBP1 CKO mice, with a q value < 0.05 (Fig 1A and S1 Table). Based on hierarchical cluster analysis, some of the 477 DEGs may participate in the regulation of the same pathway or possess similar functions ( Fig 1B and  and the CC category is associated with contractile cytoplasm, cytoplasmic part, myofibril, sarcomere, contractile fiber and contractile fiber part. The DEGs were also aligned in the KEGG pathway database. As shown in Fig 1D and Table S4, the pathways of PPAR signaling, hypertrophic cardiomyopathy (HCM), dilated cardiomyopathy (DCM), arrhymogenic right ventricular cardiomyopathy (ARVC) and calcium signaling were mainly activated.

Prediction and function analysis of the lncRNA target genes
With reference to HUGO Gene Nomenclature Committee (HGNC), the newly screened lncRNAs were divided into four types according to their positional relationship with known mRNAs (Fig 2A). To verify whether the new lncRNA meets the general characteristics, we compared the length of the transcript, the number of exons, and the length of the ORF between the lncRNA and the mRNA (Fig 2B and Fig S1A and S1B). Based on the figure, the new lncRNA conforms with the general characteristics of the transcript (transcripts with a transcript length longer than 200 nt and number of exons greater than or equal to 2 (Fig 2B and Fig S1A and S1B)). The DEGs in the lncRNA expression patterns between XBP1 knockout mice and controls were analyzed, including 334(75.7%) up-regulated genes and 107(24.3%) downregulated genes (Fig 2C and Table S5). To explore the expression patterns of these lncRNAs, a hierarchical cluster analysis of 441 differentially expressed lncRNAs was performed between XBP1 knockout mice and controls, which classified the expression patterns of the XBP1 knockout mice and wildtype mice into different clusters (Fig 2D and Table S6). We carried out target gene prediction for all lncRNAs, (i.e., predict the target genes of lncRNA through the co-location and co-expression of lncRNA and protein-coding genes). Thereafter, function enrichment analysis (GO/KEGG) of the target genes of the differential lncRNA was carried out to predict the main function of lncRNA. The prediction results of co-expression and colocation are shown in Fig 2E, Fig S2 and Tables S7 and S8). The 841 core GO terms and 45(colocation) were extracted (Fig 2E, S2 Fig, S7 and S8 Tables). In the BP category, GO terms, such as cellular metabolic process, primary metabolic process, metabolic process, organic substance metabolic process, cellular macromolecule metabolic and nitrogen compound metabolic process, were functionally enriched. The CC category was found to be associated with cell, cell part, intracellular part, intracellular, organelle, intracellular organelle and cytoplasm which the MF category was found to be associated with binding, protein binding, ion binding, catalytic activity, organic cyclic compound binding and heterocyclic compound binding. To further focus on the function of these lncRNAs, the KEGG pathways were analyzed. Based on the results, the differentially expressed lncRNA is involved in 4(coexpression) (Fig 2F and Table S9) and 4(co-location) (Fig S3 and Table S10). KEGG metabolic pathways, such as Parkinson's disease (PD), Huntington's disease (HD), carbon metabolism and dilated cardiomyopathy (DCM), were identified ( Fig 2F). The prediction result of colocation is shown in Fig S4 and Table S10.

Overview of the circRNA sequencing data
To understand the roles of XBP1-associated circRNAs in cartilage development, we performed circRNA sequencing using the rRNA-depleted chondrocyte samples of wildtype mice (CON) and XBP1 CKO mice. Clean reads from the 6 libraries were used for circRNA identification.
After a series of selection, 1651 novel circRNAs, which were widely distributed on 21 chromosomes were obtained and termed as mmu_circ_0000004 to novel_circ_0003429 ( Fig   3A, Fig S4 and Table S11). A size distribution analysis revealed that the length of most circRNAs ranged from 200 to 400bp (Fig 3B). By counting the sources of circRNA in all samples, we found that 95.08% of circRNAs were composed of exons, while 3.34% and 1.58% were located in the intronic and intergenic regions, respectively (Fig 3C). To identify circRNAs that potentially participated in cartilage development, the differences in circRNA expression patterns between XBP1 CKO mice and wildtype mice were analyzed. As shown in Table S12, a total of 15 differentially expressed circRNAs were observed relative to the control, including 10 upregulated genes and 5 downregulated genes (Fig 3D). To explore the expression patterns of these circRNAs, a hierarchical cluster analysis of 15 differentially expressed circRNAs was performed between the XBP1 CKO mice and wildtype mice, which classified the expression patterns of samples from XBP1 CKO mice and wildtype mice into different clusters (Fig 3E).  Table S13). KEGG pathways analysis, which revealed the function of these circRNAs, showed that the pathways of cell cycle, and metabolic pathway were activated ( Fig 3G and Table S14).

Overview of the miRNA sequencing data
We identified miRNAs in XBP1 CKO mice and explored their expression patterns. We found that 90.01% to 91.27% of small RNA reads were mapped onto the mouse (mus_musculus_Ensembl_97) genome (Table S15). The length distribution of sRNA was counted, and was found to range from 18 to 35 nt, with a peak of 22 nt, followed by 21 nt and 23 nt, respectively (Fig 4A and Fig S5). The statistics and annotations of all small RNAs and various types of RNAs are summarized in Fig 4B and S6 Fig. We proceeded to predict the new miRNA. M_K_1 as an example, a total of 12,993,994 miRNAs were identified, including 6,126,475 known miRNAs and 830 novel miRNAs (Fig 4B). The predicted new miRNA and the comparison of each sample sRNA are presented in Table S16. The expression levels of new and known miRNAs in each sample were determined, and the TPM (number of transcripts per million clean tags) density distribution was used to verify the gene expression pattern of the sample as a whole (Fig 4C). The expression of miRNAs in different groups was compared, and the overall expression patterns of miRNAs in these groups were found to be highly consistent ( Fig 4C). Compared with the miRNA expression levels of the control group, 6 miRNAs showed significantly differential expression (p < 0.05), including 1 significantly up-regulated miRNA and 5 significantly down-regulated miRNA (Fig 4D). To explore the expression patterns of these miRNAs, a hierarchical clustering analysis of 6 differentially expressed miRNAs was performed. Based on the results of this analysis, the cartilage samples of XBP1 CKO mice and control mice were divided into different clusters (Fig 4E). After further evaluation, we found that these miRNAs target 43,617 transcripts. According to the correspondence between miRNAs and their targets, we performed GO and KEGG enrichment analysis on each set of differentially expressed miRNA target genes. GO analysis showed that these target genes were mainly enriched in 468 GO term processes (p < 0.05) (Fig 4F). Analysis of the KEGG pathways revealed that these miRNA target genes were involved in the cancer pathway, Hepatitis B and caffeine metabolism, thyroid hormone signaling pathway, and mucin type O-glycan biosynthesis (Fig 4G).

Construction of the potential circRNA-miRNA and lncRNA-miRNA network
Competitive endogenous RNA (ceRNA) analysis is directly based on the ceRNA hypothesis, centered on miRNA, and analyzed based on the correlation between the expression levels of different molecules in multiple samples. Further, this assessment does not rely on difference analysis [48]. To clarify the functional relationship between lncRNA and its target miRNA, and circRNA and its target miRNA, we conducted a comprehensive analysis of the interaction between lncRNA and its target miRNA, and circRNA and its target miRNA. According to the expression values of miRNA and circRNA, 5,882 miRNA-circRNA interactions were identified (absolute correlation coefficient > 0.7 and p-value <0.05), with 974 miRNAs and 1,627 circRNAs, respectively (Table S17). Among these interactions, 5,055 were determined to be positively correlated, while the remaining 827 were negative (Table S17). Similarly, 272,329 miRNA-lncRNA interactions were identified (absolute correlation coefficient > 0.7 and p-value <0.05), including 1,076 miRNA and 37,320 lncRNAs, respectively (Table S18).

Construction of the potential miRNA-mRNA network
We used the miRNA software to construct a regulatory network of miRNA and its corresponding target mRNA. Collectively, we identified 14141 miRNA-mRNA interactions, involving 907 and 8758 miRNA and mRNAs, respectively (Table S21). Most of these miRNAs were found to target multiple mRNAs. For example, mmu-miR-709, mmu-miR-149-3p, mmu-miR-1249-5p, mmu-miR-7648-3p and mmu-miR-328-5p were found to have 177, 195, 166, 174, and 68 target mRNAs, respectively. Further, 136 types (14.99%) of miRNAs were found to only target one mRNA. Many mRNAs were related to more than one miRNA. For example, ENSMUST00000108883 was targeted by 2 miRNAs, including mmu-miR-3154 and mmu-miR-370-5p. Such findings indicate a complex miRNA-mRNA regulatory network in the process of cartilage development, the involvement of some miRNAs in cartilage development through the directly targeting of the mRNA of some cartilage growth related marker genes.
Other miRNAs were found to participate in cartilage development through the network of miRNA-lncRNA or miRNA-circRNA.

GO and KEGG analysis of lncRNA and circRNA regulatory networks
We performed GO and KEGG analyses to evaluate the function of DEGs in the network. In the circRNA regulatory networks, GO analysis revealed that there were 119, 24, and 4 enriched GO terms with statistical significance (p < 0.05) in the BP, CC, and MF, respectively( Fig. 5A and Table S24). For BPs, the DEGs were enriched in cellular component organization, negative regulation of biological process, location and macromolecule modification. For CC, the top enriched items were intracellular, intracellular organelle, cytoplasm and membrane-bounded organelle. Finally, for MF, binding and protein binding were enriched. KEGG pathway enrichment analysis was also conducted to characterize the target genes. Accordingly, the top significantly enriched pathways were identified to be related to inositol phosphate metabolism and phosphatidylinositol signaling system (Fig 5B and Table S25). In the lncRNA coexpression networks, GO analysis revealed 298, 67, and 24 enriched GO terms with statistical significance (p < 0.05) in the BP, CC, and MF, respectively( Fig 5D and Table S27).
We selected four circ-RNAs, including mmu_circ_0000172 (log2 fold change=2.   (A)XBP1 deficiency in cartilage leads to differences in various of secreted factors and complete transcriptome RNA of Chondrocytes, including all kinds of ncRNAs with differential expression and regulatory profiles [49]. (B) The coordinated regulation between different tissues and organs is performed by various of secreted factors and full spectrum of transcriptome RNA, including mRNA and ncRNAs. The differential ncRNAs and secreted factors caused by XBP1 deficiency in chondrocyte may cause abnormalities in different signaling pathways of different organs and tissues through blood and body fluid circulation [50][51][52], including Wnt, insulin, AMPK, and adipocytokine signal pathway. Then these differential signal pathways lead to abnormalities in the functions of other different tissues and organs, including the nervous system, heart, bones, thyroid, and intestines. Such abnormalities, in turn, result in many kinds of different diseases. The black marks in the joint cavity represent various different ncRNAs and secreted proteins from chondrocytes.

Discussion
Abnormal differentiation of chondrocytes will lead to several cartilage or bone developmental diseases, which are closely related to the molecular mechanism of their differentiation and development [53]. Previous studies have shown that the IRE1/XBP1s signaling pathway is involved in cartilage formation and bone growth [54]; however, the regulation mechanism of XBP1 and it's associated non-coding RNA involved in cartilage development remain unclear.
In the present study, we established XBP1 -/cartilage-specific knockout mice (XBP1 CKO) and conducted whole transcriptome sequencing using cartilage tissue obtained from XBP1 With the development of molecular biotechnology, non-coding RNA, a type of RNA that cannot encode proteins, was considered to be noise in the genome in the past, and it plays a vital role in various biological processes [55]. Some studies found that the expression disorder of lncRNA or circRNA was related to the occurrence of bone development disease, such as lncRNA -OG [23] and circRNA_0001052/miR-124-3p [27]. However, whether and how XBP1 splicing regulates the formation and function of ncRNAs remain unknown. Further, whether and how XBP1 modulates cartilage development through these ncRNAs-mediated regulatory networks need to be clarified. Herein, some lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks related to XBP1 deficiency were constructed. These networks indicated that circRNA and lncRNA may play a nodal regulatory role. A single circRNA and lncRNA can be associated with multiple mRNAs and multiple identical miRNAs. For example, 14 circRNAs (mmu_circ_0001376, mmu_circ_0001632, mmu_circ_0000713, mmu_circ_0000714, novel_circ_0002622, novel_circ_0002301, novel_circ_0002887, novel_circ_0002720, novel_circ_0000181, novel_circ_0000182, novel_circ_0000199, novel_circ_0000280, novel_circ_0000108 and novel_circ_0000140) were found to have binding sites for mmu-miR-99b-5p, which targets ENSMUSG00000031849 (log2 fold change=4.59) (Table S19). Hsa circ_0000729 was reported as a potential prognostic biomarker in lung adenocarcinoma, and hsa-miR-3154 serves as a potential prognostic biomarker for cervical cancer [56,57]. MicroRNA is a type of small RNA molecule, that is usually combined with a short complementary sequence located in the 3'_UTR region of mRNA to regulate the expression of target mRNA. MiRNA can directly regulate gene expression by interacting with not only lncRNA, but also circRNA and mRNA. Furthermore, both lncRNAs and circRNAs can regulate microRNAs through the "sponge" effect. LncRNAs, miRNAs, and circRNAs can also modulate gene expression by binding to the mRNA or mRNA splicing [12,[58][59][60].
LncRNAs are a highly heterogeneous class of RNAs that can be transcribed from any part of protein-coding genes and intergenic regions in sense or antisense orientation. The functions of lncRNA in eukaryotic cells are very diverse, including roles in high-order chromosome dynamics, telomere biology, and subcellular structure organization [10,11,61]. We found that lncRNAs (ENSMUST00000206990 (log2 fold change=0.98)) possessed binding sites for mmu-miR-23a-5p, which targets ENSMUSG00000000126 (log2 fold change=3.19) (Table   S20). MiRNAs are thought to modulate the translation process of more than 60% of proteincoding genes in eukaryotic cells, which are involved in the regulation of many biological processes, including proliferation, apoptosis, differentiation, and development. Different miRNAs have different regulatory processes; some miRNAs control specific individual targets, whereas others serve as the main regulators of one process. Therefore, pivotal miRNAs can simultaneously regulate the expression levels of multiple genes, and different types of miRNAs coordinate their targets [10,11,61]. In XBP1 KO chondrocytes, mmu-miR-709, mmu-miR-149-3p, mmu-miR-1249-5p, mmu-miR-7648-3p and mmu-miR-328-5p had 177, 195, 166, 174 and 68 target mRNAs, respectively. In summary, these RNAs play an important role in cartilage development by modulating XBP1 expression. In addition, the biological function of target genes was identified based on GO annotation and KEGG pathway enrichment analyses.
All is known that the coordinated regulation between different tissues and organs is performed Prof. Guo FJ designed the manuscript and had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. All authors approved the final manuscript prior to submission.

Data availability
There are no restrictions on sharing the data. All data are available from the corresponding author upon reasonable request.