Expression Pattern of CircRNA in Naringin-induced Fate Dicision of BMSCs


 Objective

To investigate the effect of Naringin (NG) on the differentiation of bone marrow-derived mesenchymal stem cells (BMSCs) in vitro.
Methods

BMSCs with different concentrations of NG were treated. The cell morphology, proliferation, toxicity analysis of each group, the osteogenic differentiation of BMSCs, and the expressions of markers related to osteogenic differentiation of BMSCs, osteogenic marker protein and CircRNAs were observed by Cell Counting Kit-8(CCK-8)method, Reverse Transcription-Polymerase Chain Reaction༈RT-PCR༉ method, etc. We used Miranda to predict the microRNA targets of CircRNAs and the targets of CircRNA-miRNA interaction, and constructed a CircRNA-microRNA co-expression network by analyzing the correlation between differentially expressed CircRNAs and microRNAs. The diagram of CircRNA-microRNA network was created with Cytoscape.
Results

The research results showed that NG of 1, 5, 10, 50, 100, 200µg/mL played an important role in accelerating the proliferation of BMSCs in vitro. It is found that osteocalcin (OCN), Sp7 transcription factor(SP7) and runt-related transcription factor 2 (RUNX2) were significantly enhanced in mRNA and protein expression levels. Furthermore, 633 CircRNAs showed an upward trend of expression, while 762 CircRNAs showed a downward trend of expression. Compared to those in the control group, the expression level of CircTll 1 decreased significantly (P < 0.01), the expression level of CircFkbp5 increased (P < 0.05), and the expression level of Circzbtbl6-2 and CircCped1 increased significantly (P < 0.001). In addition, we analyzed the target mRNAs of miR-342-3p by GO and Panther pathways and found that Circ14046, Circ10410 and Circ7511 were not only closely related to calcium signaling pathway, but also to the regulation of actin cytoskeleton involved in cell growth and differentiation.
Conclusion

This paper confirmed that Circ14046, Circ10410 and Circ7511 were related to osteogenic differentiation. NG may induce osteogenic differentiation of BMSCs through Circ14046/Circ10410/Circ7511 targeting miRNA-mRNA axis.


Introduction
Drynaria rhizome is a drug in traditional Chinese medicine which is widely used in the treatment of orthopedic diseases [1]. Naringin (NG) is the main active component of Drynaria rhizome [2]. Modern pharmacology and our team's previous research have shown that NG can be used to treat a variety of orthopedic diseases, such as fracture, femoral head necrosis and osteoporosis by promoting the proliferation and osteogenic differentiation of bone marrow mesenchymal stem cells (BMSCs), up regulating the expression of osteocalcin and bone morphogenetic pritein2 BMP-2 , promoting osteoclast apoptosis and proliferation, stimulating angiogenesis, inhibiting in ammatory reaction and so on [3,4] .
Therefore, NG is a promising drug.Circular RNAs (CircRNAs) are endogenous, non-coding, covalently closed circular RNAs that are different from linear RNAs [5]. They are highly stable, evolutionarily conservative and speci c in tissue expression [6]. At the same time, as a sponge functional molecule of miRNAs, CircRNAs play an important regulatory role, which is mainly re ected in the post transcriptional regulation. Studies have shown that CircRNAs may participate in the pathological process of a variety of orthopedic diseases, such as osteoarthritis, osteoporosis, and so on [7]. They can promote the proliferation of osteoblasts and the coupling process of osteogenesis and angiogenesis to treat fracture and osteoporosis [8][9]. Therefore, CircRNAs may become a biomarker for diagnosis and prognosis and a potential therapeutic target in orthopedics.
The purpose of this study is to detect the expression Pattern of CircRNAs in osteogenic differentiation of BMSCs under the intervention of NG.

BMSCs culture
The 4-6 week old male sprague-dawley(SD) rats used in the experiment were provided by Laboratory Animal Center of Sun Yat-sen University. SD rats were anesthetized, disinfected, aseptically treated with femur and tibia, immersed in PBS containing 2% double antibody; repeated rinsing of the medullary cavity, collection of medullary lavage uid; centrifuged at 1500 rpm for 5 minutes; The inoculation density was kept at 6-8×10 6 /mL into a 10 cm 2 culture dish, and cultured in an incubator at 37 ° C, 5% CO 2 saturated 90% humidity; After 48h, the cells were replaced by liquid exchange, and the cells were replaced by 2-3d until the cells reached 80%-90% fusion. After digestion with 0.25% trypsin containing 0.02% ethylenediamine tetraacetic acid (EDTA), they were passaged at a ratio of 1:3. Recorded as P1. After the passage, the liquid was changed every 2-3d, and the cells could continue to pass after reaching the basic fusion [10].

Cell proliferation assay
Cell proliferation was estimated by CCK-8 assay. BMSCs were inoculated into 96-well plates (8000/well), and the plates were placed in an incubator (37 ° C, 5% CO 2 ) for pre-cultivation for 24 h. The control group was not treated. The experimental groups were added with NG (0.001, 1, 5, 10, 50, 100, 200 μg /ml), cultured in an incubator for 48 h. Add 20 μL CCK-8 solution to each well; incubated the plates for 2 hours in the incubator; The absorbance was measured at 450nm by spectrophotometer.

ALP vitality test
The activity of alkaline phosphatase (ALP) was detected by ALP staining (ALP assay kit, KGI Bio-KGA353) on the 7th day of NG induction.

RT-PCR maintains the principle of timing and quanti cation
The total RNA of cultured cells (with 10μg/ml, 100μg/ml or without NG for 10d) were extracted by TRIzolRNA extraction reagent (Invitgen, CarlsbadCA). The integrity and concentration of RNA were detected by high purity total RNA rapid extraction kit (rotary column) (BioTeke, S Hangzhou, CH) and Nano-200 super nucleic acid analyzer (Ausheng, Hangzhou, CH)). M-MLV reverse transcriptase (Promega) and random primers were used for reverse transcription. Using the primers described in Table  1. Useing HiScriptIIQ RT SuperMix for qPCR (+ gDNAwiper) (Vazyme, Nanjing, CH)) for all qPCR reactions.
GAPDH was an internal control. Statistical analysis of the data using the ΔΔCt method, the relative expression level of circRNAs was replaced by 2-ΔΔCT. In addition, RNaseR restriction endonuclease digestion reaction was also carried out. Total RNAs (5Mg) with or without RNaseR (3U/mg) (Epicentre Bio Technologies) was incubated at 37 ℃ for 15 min, and then RNA, was puri ed with phenol-chloroform for RT-qPCR. BMSCs were seeded in a 6-well plate (5×10 3 /well) and treated with 100μg/mL NG for 27 days. Adding Ob xative to the cells and x for 15 minutes at room temperature. Washed with 2mL phosphate buffer saline(PBS) 3 times, 3 minutes each time. The Alizarin Red S staining solution was added at room temperature for 30 minutes, and the dyeing solution was washed with distilled water. Under a 40x inverted microscope, 5 elds were randomly selected to take photos and observe the staining.
2.6. RNA-seq Total RNA was extracted by TRIzol (Invitrogen, Carlsbad, CA, USA) and then puri ed by RNeasy Mini Kit (Qiagen, Valencia, CA, USA). In short, the rst strand of Poly (A) mRNA was prepared by fragmentation of cDNA with random hexamers. After the second strand cDNA synthesis, terminal repair, adding a single A base, joint connection, agarose gel separation and ampli cation of about 200bp cDNA, the cDNA was sequenced by IlluminaHiSeq4000 sequencing system. Using the R software package to the quantile normalization and process subsequent data. The dissimilarity display circRNAs were screened according to the critical value of folding change and the statistically signi cant q value (q value ≤0.05, fold change ≥2.0) or by volcano map.

Analysis of gene ontology and its pathways
Using DAVID to elucidate the recessive role of circRNAs parent genes (Database for Annotation, Visualization and Integrated Discovery). The functions of parental genes were predicted by Gene Ontology (GO)functional annotation. Next, the role of genes was summarized into three types by type: biological process (BP), molecular function (MF) and cellular component (CC). The datas calculated by GO were shown in the column classi cation chart. Using Kyoto Encyclopedia of Genes and Genomes (KEGG) to further explore the related pathways of different expressed circRNAs parent genes.

Description of circRNA/miRNA interactions
The interactions between miRNA targets and circRNA/miRNA were based on Miranda, and the circRNAs were predicted by Arraystar's miRNA target prediction software [11]. Constructing a CircRNA-miRNA coexpression network which based on the correlation analysis of differentially expressed circRNAs and miRNAs. In order to construct a circRNA-miRNA network, we used Arraystar's tools for predicting miRNA targets to search for miRNAs response elements to circRNAs, and then selected miRNAs according to seed matching sequences. Drawing the circRNA-miRNA network diagram by Cytoscape3.7.1 [12].

The Osteogenic effects of NG on BMSCs
NG promoted the proliferation of BMSCs in vitro. Compared with the control group, there was no statistically difference in 0.001μg/mL treatment group P>0.05), While 1, 5, 10, 50, 100, 200μg/mL treatment groups optical density(OD) values were higher than the control group, there were very signi cant differences, P <0.001 Fig.1A). NG increased the ALP activity of BMSCs in vitro. Compared with the control group, the precipitation of brown-black cobalt sul de particles in BMSCs supplemented with 10μg/ml and 100μg/ml NG were signi cantly increased, and the amount of cobalt sul de precipitated were proportional to the ALP content Fig.1B).
Under the action of NG, the mRNA and protein expression levels of OCN, Runx2, and SP7 in BMSCs increased signi cantly Fig.1C). NG could induce osteogenic differentiation of BMSCs [13][14]. In order to further explored the effect of NG on the expression of bone-related mRNA and protein, we uesd the qRT-PCR to detect the mRNA and protein expression levels of OCN, Runx2 and SP7 in BMSCs. Our results showed that mRNA and protein expression levels of OCN, Runx2 and SP7 in the10/100μg/ml experimental group increased on the 10th day, while the control group did not increase.Alizarin red staining was performed after 21 days of NG treatment. Compared with the control group, the number and size of orange calci cation nodules in BMSCs increased signi cantly after NG treatment Fig.1D .

Distinguish differentially expressed circRNAs
We used RNA-seq technique to study the CircRNA expression characteristics of NG-treated BMSCs. A total of 15,282 CircRNA targets were found in two pairs of test articles. The reading range of linear and circular RNAs can be clearly seen in the gure ( Fig.2A), and the speci c arrangement of circRNAs in human chromosomes can also be found (Fig. 2B). Most circRNAs are transcribed from protein-coding exons, some from introns, and some are intergenomic. The differentially expressed circRNAs were folded and ltered and expressed in a scatter diagram (Fig. 2C). In addition, volcanic map ltering found signi cant differences in circRNAs between the two groups (Fig.2D).
The results showed that 1395 circRNAs were regulated by Q value ≤ 0.05 and folding change ≥ 2.0. Among them, 633 circRNAs showed an upward trend, while 762 circRNAs showed a downward trend. There were 129 more down-regulated circRNAs than up-regulated circRNAs. Based on log 2 FC, we listed in Table 2 the top 10 circRNAs with increased and decreased expression in this experiment. CircRNA can adsorb microRNA, so we could nd out which microRNA was bound by circRNA through prediction, namely, the correlation analysis between circRNA and microRNA. Fig. 3 shows the top400(A) and top500(B)_diff_circRNA_target_miRNA_network. Note: CircRNA ID: CircBase contains CircRNA ID (http://circbase.mdc-berlin.de). CircRNA type: CircRNAs can be expressed in 5 different forms: "exon", "intron", "antisense", "intra-gene" and "inter-gene". P-value: P-value calculated according to paired t-test. FDR: Analysis of FDR through Benjamini-Hochberg FDR. Fold change: the absolute ratio of normalized intensities between two conditions (without logarithmic scale).

Analysis of differential expression level of CircRNAs by qRT-PCR
Six CircRNAs( CircT111, CircBmper, CircAnkrd44, CircFkbp5, CircCped1 and circZbtb16-2)were veri ed by qRT-PCR in the samples of NG treatment group and control group. In cDNA, CircT111, CircBmper, CircAnkrd44, CircFkbp5 and CircCped1 Both circular primers and linear primers could amplify a single band. circZbtb16-2 circular primers had two bands, indicating that there might be other cyclization sites (other cyclization sites have been veri ed by sequencing), and linear primers could amplify a single band. In gDNAs, ring primers had no ampli cation bands, while linear primers could amplify a single band, indicating that the ring was expressed in the sample. The veri cation results shown that CircT111, CircBmper, CircAnkrd44, CircFkbp5, CircCped1 and circZbtb16-2 were expressed in the samples (Fig. 5A).
In addition, the Sanger sequencing of RT-PCR products ampli ed with different primers further con rmed the post-splice joint (Fig.5B). The relative quantitative results of the six CircRNAs were listed in Fig.8C. Compared to the control group the expression of circZbtb16-2 and circCped1 increased signi cantly p<0.001) , the expression of circFkbp5 increased (p<0.05 , the expression of circTll1 was signi cantly decreased p<0.01 and the expression of circBmper and circAnkrd44 showed no signi cant changes. 3.9. Prediction and annotation of circ14046(circZbtb16-2)/ circ10410(circCped1)/ circ7511(circFkbp5) targeted miRNA-mRNA network The three upregulated CircRNAs veri ed in the previous experiment were selected, and the predicted top20 miRNAs were selected for intersection. One miRNA was found to be the downstream regulated by the three CircRNAs (rno-miR-342-3p) (sig. 6A). To further analyze the mRNA that rno-miR-342-3p may target, we used two bioinformatics websites to make predictions (miRanda / miRWalk). the prediction results were taken as the intersection, which had 39 genes in total miRanda (http://www.microrna.org/microrna/home.do) 515 genes and miRWalk (http://mirwalk.umm.uniheidelberg.de/) 479 genes (sig. 6B). Subsequently, the GO and KEGG pathway enrichment analysis of these 39 genes were found to be signi cantly enriched in the Calcium signaling pathway and Regulation of actin cytoskeleton (sig. 6C/6D). We selected four genes enriched in these two pathways to make interaction network diagrams (sig. 6E).

Discussion
CircRNA is a non-coding RNA with the ability to regulate gene expression, mainly including CircRNA from introns, CircRNA from exons and CircRNA from introns. CircRNA exists widely in the biological world, and regulates various life activities after transcription or transcription [15][16][17].Recent studies have found that certain CircRNAs can act as miRNA sponges, participate in the regulation of alternative splicing, competitively inhibit miRNA transcriptional regulation, and increase the complexity of RNA regulatory network [18], and play an important role in the occurrence and development of orthopedic diseases.Including intervertebral disc degeneration, osteoarthritis, osteosarcoma and so on CircRNAs have been proven to regulate the expression of related genes and the production of important proteins through the interaction of CircRNA-miRNA, and participate in the osteogenic differentiation of BMSCs, osteoblasts and osteoclasts Differentiate, regulation of bone metabolism, etc. [19][20].Similar to miRNA and LncRNA, CircRNA has gradually become a research hotspot in the RNA eld, and is considered a new biomarker and therapeutic target [21][22].
Our experiment results showed that 1395 CircRNAs were regulated by Q value ≤ 0.05 and folding change In addition, we also found that the expression of ALP was signi cantly increased in BMSCs added with NG (10μg/ml and 100μg/ml) (P<0.05). At the same time, compared with the control group, the number and size of orange calci ed nodules in BMSCs increased signi cantly after NG treatment (P <0.05).At present, it has been reported that ALP activity and calcium content play an important role in the reproduction, differentiation and osteogenic ability of osteoblasts, which are usually used to evaluate the function of osteoblasts [23]. At the same time, when osteoblasts are cultured in vitro, they will cause the deposition of calcium salts, thereby forming calci ed nodules. The more calcium salts are deposited, the better the osteogenic differentiation [24].
So far, there were few studies on the role of CircRNAs in the osteogenic differentiation of BMSCs. CircRNAs may regulate linear RNA transcription and protein synthesis by acting as miRNA sponges to regulate gene expression [25][26].To date, the role of circRNAs in BMSCs osteoblastic differentiation has been rarely investigated. CircRNAs can act as miRNA sponges to regulate gene expression, thereby regulating linear RNA transcription and protein synthesis [25][26].Compared with linear isoforms, CircRNAs has a higher expression level and contains abundant miRNA binding sites. Therefore, the role of CircRNA in the osteogenic differentiation of BMSCs may be related to the miRNA-mediated role. QRT-PCR con rmed that the four CircRNAs in the NG-treated cells were all changed, and the expressions of CircRNA14046, CircRNA10410 and CircRNA7511 were all up-regulated (P<0.001/P<0.05).Selecting CircRNA14046, CircRNA10410 and CircRNA7511 veri ed in previous experiments, and selecting the top 20 predicted miRNAs for crossover. It was found that Rno-miR-342-3p is downstream regulated by three CircRNAs. Rno-miR-342-3p could be targeted by CircRNA14046, CircRNA10410 and CircRNA7511 in the CircRNA-miRNA network.CircRNA14046 was predicted to interact with rno-miR-342-3p, rno-miR-3541, rno-miR-874-3p and rno-miR-351-5p, while CircRNA10410 was predicted to interact with rno-miR-343 and rno-miR-351-5p. CircRNA7511 was predicted to interact with rno-miR-134-3p and rno-miR-298-5p.According to the assumption , CircRNA14046, CircRNA10410 and CircRNA7511 may serve as effective miRNA sponges to predict miRNA binding sites. Through the analysis of GO and Panther pathways, we found that the target mRNAs of CircRNA14046, CircRNA10410 and CircRNA7511 were closely related to the regulation of Calcium signaling pathway and Regulation of actin cytoskeleton. These signaling pathways are related to cell growth and differentiation. We selected four genes enriched in these two signaling pathways to form a network diagram of interaction. Meanwhile, compared with the control group, the mRNA expression levels of OCN, SP7 and Runx2 related to bone differentiation were signi cantly increased after NG treatment. Therefore, we believe that NG can participate in osteogenic differentiation induction of BMSCs through circRNA-miRNA-mRNA axis, but its speci c mechanism still needs further experimental veri cation.

Conclusion
In summary, this study con rmed that CircRNA14046, CircRNA10410 and CircRNA7511 were related to osteoblast differentiation of BMSCs. NG may induce BMSCs to differentiate into osteoblasts through CircRNA14046/CircRNA10410/CircRNA7511 targeting miRNA-mRNA axis. It will be a very valuable research hotspot in the future to study CircRNAs, as miRNA sponges to regulate the speci c molecular