- Open Access
Integrative analysis reveals a lineage-specific circular RNA landscape for adipo-osteogenesis of human mesenchymal stem cells
Stem Cell Research & Therapy volume 13, Article number: 106 (2022)
The balance between osteogenesis and adipogenesis of mesenchymal stem cells (MSCs) is critical to skeletal development and diseases. As a research hotspot, circular RNAs (circRNAs) have expanded our understanding of a hidden layer of the transcriptome. Yet, their roles during adipo-osteogenesis remain poorly described.
The identity of human MSCs derived from bone marrow and adipose were first determined by flow cytometry, cellular staining, and quantitative polymerase chain reaction (qPCR). Multi-strategic RNA-sequencing was performed using Poly A, RiboMinus and RiboMinus/RNase R methods. Integrative analysis was performed to identify lineage-specific expressed circRNAs. The structural and expressional characteristics were identified by Sanger sequencing and qPCR, respectively. The regulatory effects of adipogenesis-specific circ-CRLF1 were confirmed using siRNA transcfection and qPCR.
We generated a whole transcriptome map during adipo-osteogenesis based on 10 Poly A, 20 RiboMinus and 20 RiboMinus/ RNase R datasets. A total of 31,326 circRNAs were identified and quantified from ~ 3.4 billion paired-end reads. Furthermore, the integrative analysis revealed that 1166 circRNA genes exhibited strong lineage-specific expression patterns. Their host genes were enriched in distinct biological functions, such as cell adhesion, cytokine signaling, and cell division. We randomly selected and validated the back-spliced junction sites and expression patterns of 12 lineage-specific circRNAs. Functional analysis indicated that circ-CRLF1 negatively regulated adipogenesis.
Our integrative analysis reveals an accurate and generally applicable lineage-specific circRNA landscape for adipo-osteogenesis of MSCs and provides a potential therapeutic target, circ-CRLF1, for the treatment of skeleton-related disease.
Mesenchymal stem cells (MSCs) are multipotent clonogenic cells that existed in diverse tissues, including bone, fat, and muscle. Due to the potentials to differentiate into stromal cells and their immunoregulatory properties, MSCs are increasingly recognized as crucial contributors to tissue homeostasis and ideal candidates for cell-based therapies . The normal bone microenvironment requires the balance of osteogenesis and adipogenesis of MSCs, the imbalance of which has been linked to pathophysiologic processes of bone such as aging, obesity, osteopenia, and osteoporosis [2,3,4]. Therefore, clarifying the adipo-osteogenesis of MSCs will facilitate the application of MSC in regenerative medicine, especially for bone-related diseases.
Circular RNAs (circRNAs) are a novel and large family of non-coding RNAs (ncRNAs) with covalent loop structures formed by head-to-tail splicing [5, 6]. CircRNAs have long been considered to be by-products of transcription due to their inability to be translate into proteins and aberrant splicing events . However, accumulating evidence indicates that circRNAs are involved in the development and pathogenesis inmultiple tissues, including bone [7, 8], by functioning as an inhibitor of microRNA or proteins [6, 9], as structural components of protein complexes , or as signaling molecules to recruit proteins . For example, circRNA CDR1as could compete for miR-7-5p, which controls the expression of WNT5B by targeting its 3′-untraslated region (3’UTR), to promote adipogenesis and suppress osteogenesis of bone marrow-derived MSCs (BMSCs) . Similarly, hsa_circ_0074834 sponges miR-942-5p, which could target the mRNA of ZEB1 and VEGF, to promote the osteogenesis and angiogenesis of BMSCs; circRNA-vgll3 could regulate miR-326-5p/Itga5 axis to promote osteogenesis of adipose-derived MSCs (AMSCs) [13, 14]. CircFUT10 and hsa_circH19 could target let-7c and PTBP1 proteins, respectively, to control adipogenesis [15, 16].
Although the roles of circRNAs in the lineage determination of MSCs have been becoming more evident, only a few publications have reported the expressions with potential functions of circRNAs at a genome-wide level during adipo-osteogenesis in humans. For instance, the circRNA expression levles during osteogenesis have been reported in BMSCs using microarray or RNase R-treated RNA-sequencing (RNA-seq) [17, 18], and in AMSCs and stem cells from the periodontal ligament and apical papilla using RiboMinus-treated RNA-seq [14, 19, 20]. Adipogenesis-related circRNAs have been reported in human visceral and subcutaneous fat tissue and mouse pre-adipocytes using RiboMinus-treated RNA-seq . However, due to their low abundance compared with linear RNAs, a limited number of circRNAs were detected by only removing ribosomal RNAs during RiboMinus-treated RNA-seq , and a biased estimation of circRNA levles may be generated during RNase R treatment and multiple purification steps during RNase R-treated RNA-seq . Integrative analysis of multi-strategic RNA-seq data using appropriate algorithms could improve the accuracy of circRNA identification and quantification with a significantly reduced false discovery rate .
In the present study, we produced 50 high-quality RNA-seq datasets, using RiboMinus, Poly A, and RiboMinus/RNase R methods, at five time points during adipo-osteogenesis of human BMSCs and AMSCs. Our integrative analysis identified 31,326 circRNAs shared by MSCs from different tissues, 1,166 host genes of which displayed strong lineage-specific expression patterns. Twelve circRNAs were further validated by both quantitative polymerase chain reaction (qPCR) and Sanger sequencing. Moreover, we found that levels of the circRNA, circ-CRLF1, which is derived from the CRLF1 locus, were markedly increased during adipogenesis, but not osteogenesis. Moreover, the knock-down of circ-CRLF1 significantly promoted adipogenesis. These datasets will be useful for further investigation of the roles of circRNAs during bone development and dysfunction.
Materials and methods
Human bone marrow and adipose-derived MSCs were purchased from American Type Culture Collection (Catalog PCS-500-012, PCS-500-011, ATCC, MD, USA). They were cultured in MesenCult-ACF Plus Medium (catalog number 05448, STEMCELL Technologies, Vancouver, Canada) supplemented with 1% penicillin–streptomycin. They were cultured at 37 ℃ with 5% CO2. The cells were passaged at a seeding density of 4 × 103 viable cells/cm2 when reaching 80% confluence using Animal Component-Free Cell Dissociation Kit (catalog number 05426, STEMCELL Technologies, Vancouver, Canada).
Osteogenic and adipogenic differentiation
MSCs were orientally induced under passage 8 when they reached 80% confluence. For osteogenic induction, the cells were washed with 1 × dPBS and then cultured with OsteoMAX-XF Differentiation Medium (catalog number SCM121, EMD Millipore, MD, USA) supplemented with 1% penicillin–streptomycin. The osteogenic medium was half-changed on day 2, 4, 7, 9, and 11, according to the previous study . MesenCult Adipogenic Differentiation Kit (catalog number 05412, STEMCELL Technologies, Vancouver, Canada) supplemented with 1% penicillin–streptomycin was used for adipogenic induction, and all the culture medium was replaced at the same time point mentioned above. The samples were collected on day 0, 7, and 14 for the experiments.
The flow cytometry was performed with antibodies for positive markers: anti-CD73 FITC with related isotype (Catalog Number 561254, 554679, BD Pharmingen, CA, USA), anti-CD105 APC with related isotype (catalog number FAB10971A, IC002A, R&D systems, MN, USA), anti-CD166 PE with related isotype (Catalog Number 343904, 400111, BioLegend, CA, USA); and negative markers: anti-CD14 PE with related isotype (Catalog Number 555398, 555574, BD Pharmingen, CA, USA), anti-CD34 APC with related isotype (Catalog Number 555824, 555751, BD Pharmingen, CA, USA), anti-CD45 APC with related isotype (Catalog Number 555485, 555751, BD Pharmingen, CA, USA). The MSCs at passage 8 were digested into single-cell suspension and washed with 1 × dPBS, and then incubate with those fluorescent-conjugated antibodies for one hour at room temperature. After washing with 1 × dPBS, the cells were detected by flow cytometer (Beckman Coulter DxFLEX, CA, USA) with at least 10,000 events. Data were analyzed by FlowJo (TreeStar, OR, USA).
Alizarin Red S and oil Red O staining
The cells were fixed in 4% paraformaldehyde for 20 min at room temperature and then washed twice with 1 × dPBS. To detect the matrix mineralization during osteogenesis, the cells were stained with 1% Alizarin Red S (pH 4.2; catalog number G1452, Solarbio, Beijing, China) for 15 min. To visualize the intracellular lipid deposits, the cells were incubated with 60% isopropanol for 5 min, then stained with Oil Red O solution (catalog number G1262, Solarbio, Beijing, China) for 15 min. The stained cells were examined by light microscope (EVOS XL Core, 10 × magnification, Life Technology, PA, USA).
Cells grown in 6-well plates were washed 3 times with PBS, and fixed for 15 min in 4% paraformaldehyde at room temperature and treated with 0.1% Triton X-100 for 15 min and then blocked for 60 min with 5% goat serum in PBS. After that, cells were incubated with rabbit anti-RUXN2 (catalog number 12556, Cell Signaling Technology, MA, USA; dilution 1:800) or rabbit anti-CEBPA (catalog number 8178, Cell Signaling Technology, MA, USA; dilution 1:200) primary antibody at 4 °C overnight and secondary antibodies (Alexa Fluor 488, catalog number 4412, Cell Signaling Technology, MA, USA; dilution 1:500) for one h at room temperature, followed by washing and staining with DAPI (2 μg/ml, catalog number D8417, Sigma-Aldrich, MO, USA) for 5 min. Then stained cells were washed two times with PBS, and images were taken using a Leica microscope (Leica, DMi8, Wetzlar, Germany) with its software LAS X. The average optical density (AOD) of immunoreactivity were normalized by AOD of DAPI in the same field using Image J.
RNA isolation and quantitative PCR
Total RNA was isolated using TRIzol (Catalog Number 15596026, Invitrogen, CA, USA) and RNeasy Mini Kit (Catalog Number 74104, QIAGEN, CA, USA) according to manufacturers’ instructions. For RNase R treatment, 2 μg of total RNA were incubated with or without 3 U/μg RNase R (Catalog Number RNR07250, Epicentre Technologies, WI, USA) for 20 min at 37 ℃, for 5 min at 70 ℃. Subsequently, complementary DNA was produced using the PrimeScript RT Reagent Kit with gDNA Eraser (catalog number RR047A, Takara, Shiga, Japan), and quantitative real-time reverse transcription polymerase chain reaction (qPCR) was performed using the Real-Time PCR detection system (Thermo Scientific SteponePlus, MA, USA) with TB Green Premix Ex Taq™ II kit (catalog number RR820A, Takara, Shiga, Japan). Relative mRNA expression was calculated using the delta threshold cycle (ΔΔCT) method and normalized to GADPH expression . The PCR primers are listed in Additional file 7: Table S6. In addition, the PCR products for circRNAs were performed with sanger sequencing to validate the back-spliced junction sequences.
For RiboMinus-treated RNA-seq, the ribosomal RNAs were removed by Epicentre Ribo-zero rRNA Removal Kit (Epicentre, WI, USA). For Poly A-treated RNA-seq, the poly-T oligo-attached beads were used to purify poly A-containing RNAs. For RiboMinus/RNase R RNA-seq, the rRNA-depleted RNAs were digested with RNase R (Epicentre, WI, USA) to remove linear RNAs. Whereafter, the strand-specific libraries were performed with the treated RNAs using the rRNA-depleted RNA by NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, MA, USA) following the manufacturer’s recommendations. The sequencing was carried out using the Illumina HiSeq X Ten instrument by the commercial service of Novogene Co., Ltd. (Tianjin, China).
Low-quality sequencing reads and the reads with adapters or with > 10% of ploy-N bases were filtered out. After filtering, the clean reads were used to calculate gene expression values and identify circular RNAs. The read count and TPM (transcripts per million) values of each type of gene were determined using the fast pseudoaligner Kallisto (v0.46.0) , which pseudoaligned reads to reference gene models (GENCODE v35) . Then, the read counts and TPMs at the transcript level were summarized to the gene level. Differentially expressed genes (DEGs) were identified using the DESeq2 method . DEGs with fold-change more than two and adjusted P value less than 0.01 were retained and used in the following analysis. Gene Ontology (GO) enrichment analysis was performed using DAVID bioinformatics resources 6.8 . According to RNA-seq data and reference gene models, circRNAs were identified using CIRI2 . The expression levels of circRNAs were estimated by CIRIquant with RNase R effect correction model based on RNase R treated and untreated samples . CircAltas 2.0 database was used as the source of the reference circular RNAs .
A predesigned siRNA targeting human circ-CRLF1 (si-circRNA) and negative control (scramble siRNA, si-scramble) were purchased from GenePharma Co., Ltd. (Shanghai, China). The sequences used are listed in Additional file 7: Table S6. MSCs were transfected with siRNA at a final concentration of 20 nM on day 0, 2 and 4 during adipogenesis, using Lipofectamine RNAiMAX (Invitrogen, Cat.No. 13778075), according to the manufacturer’s protocol. Cells were harvested on day 7 post-transfection.
Prism software 5.01 (GraphPad Software, Inc., San Diego, USA) were used to for statistical analysis. A two-tailed Student’s t-test was used to compare mean values between two groups. The statistical significance in the comparison of multiple sample sets versus control was performed with Bonferroni’s multiple comparisons test after one-way ANOVA test.
Adipo-osteogenic characteristics of inducted human BMSCs
Human BMSCs were expanded and induced using osteogenic or adipogenic medium for 14 days to establish an adipo-osteogenesis system. The expanded BMSCs maintained the typical expression pattern of their surface markers: high CD73, CD105 and CD166 expression levels, and low CD14, CD34 and CD45 expression levels (Fig. 1b; Additional file 1: Fig. S1a). Histological staining indicated the accumulation of extracellular matrix mineralization during osteogenesis and the accumulation of intracellular lipid droplets during adipogenesis (Fig. 1a). Expression levels of osteogenic marker genes, including RUNX2, ALPL, and MMP1, and adipogenic marker genes, including PPARG, CEBPA, and FABP4, were significantly increased during osteogenesis and adipogenesis, respectively (Fig. 1c, d). Therefore, human BMSCs displayed typical adipogenic and osteogenic lineage characteristics after oriented induction.
Advantages of multi-strategic RNA-seq analysis of induced BMSCs
To generate the comprehensive RNA landscape of induced BMSCs, we sequenced 30 samples from five time points during adipo-osteogenesis of human BMSCs, using three RNA-seq library types: (1) RiboMinus-treated library for total RNA detection, (2) Poly(A) enrichment library as standard RNA-seq library, (3) RiboMinus/RNase R-treated library for circRNA detection (10 libraries for each type) (Fig. 2a; Additional file 2: Table S1). For the three types of RNA-seq libraries, an average of 62 million clean reads were obtained, 94% of which were mapped to the human genome (version hg38) (Additional file 1: Fig S2a; Additional file 2: Table S1). The expression profiles of both protein-coding and lncRNA genes of two replicates from each library type and time point exhibited high correlations, indicating excellent repeatability of our experiments (Additional file 1: Fig S2b, c). The expressional changes of adipo-osteogenic marker genes showed similar trends when detected by RNA-seq or qPCR (Fig. 1c; Additional file 1: Fig S2d). Furthermore, we found that an average of 87% of protein-coding genes and 68% of lncRNA genes annotated by the GENCODE database were detected in RiboMinus data, whereas these proportions were 80% and 44% in Poly A data, 74% and 34% in RiboMinus/RNase R data (Fig. 2b). An average of 97% of the detected genes in Poly A data belonged to protein-coding genes, whereas the proportions were 54% in RiboMinus data, and 47% in RiboMinus/RNase R data (Fig. 2c). In addition, the highest portions of back-spliced junction (BSJ) reads were detected in RiboMinus/RNase R data (Additional file 1: Fig S2e; Additional file 2: Table S1). Therefore, both RiboMinus- and RiboMinus/RNase R-treated RNA-seq methods preserved the balanced RNA proportions of protein-coding and lncRNA genes, while the RiboMinus method detected slightly more annotated genes and the RiboMinus/RNase R method detected extraordinarily more BSJ reads, which is valuable for circRNA identification. The Poly A data displayed a strong bias for the enrichment of RNAs from protein-coding genes. These results demonstrated that the three library types exhibited distinct advantages in providing an overall RNA landscape during adipo-osteogenesis of human BMSCs.
Integrative identification of circRNAs in inducted BMSCs
To improve the accuracy of circRNA identification, an integrative pipeline was performed to utilize the advantages of multiple dataset types. The predicted circRNAs were detected using the CIRI2 tool for each library type, while predicted circRNAs from Poly A-treated libraries were considered as false-positive signals and filtered out. The expression values of the remaining circRNAs in the RiboMinus and RiboMinus/RNase R data were corrected for RNase R efficiency using the CIRIquant tool (Fig. 3a). According to the circRNA identification results, at least eight times more circRNAs were identified in each RiboMinus/RNase R data set (16,767 circRNAs) than those in RiboMinus data (2,063 circRNAs), while an average of 14 circRNAs were identified and filtered from the Poly A data (Fig. 3b). After correction and merging, a total of 48,152 circRNAs were obtained from 7366 genes (Additional file 3: Table S2), of which 30% were identified in more than five samples and thus may be considered conserved circRNAs (Fig. 3c). Among these circRNAs, 96.2% were matched with circAtlas 2.0 annotation, 95.9% were generated from protein-coding genes, and 88.5% were from exonic regions (Fig. 3d–f).
Next, we evaluated the relationship between host gene characteristics and their circRNA numbers. The results showed that there was a higher correlation between the number of host gene exons and the number of their circRNAs than between the host gene expression levels and the numbers of their circRNAs (Fig. 3g; Additional file 1: Fig. S2f). For example, BIRC6 has 74 exons and generated 96 circRNAs. Furthermore, the circRNA correlation analysis, based on the expression levels of 7,366 host genes, showed a high correlation between two replicates at each time point from each type of library and among different lineages. This suggested that circRNA-related genes exhibited lineage-specific expression characteristics (Fig. 3h).
Lineage-specific expression of circRNA genes during adipo-osteogenesis of BMSCs
The efficient differentiation of MSCs requires genome-wide gene-specific expression at different developmental stages. To characterize the distinct patterns of circRNA expression more comprehensively, we performed differential gene expression analysis, using the circRNAs we identified as a reference, during adipo-osteogenesis of BMSCs. In total, 2,077 circRNA genes were found to be significantly differentially expressed (adjusted P value < 0.05 & |log2 FoldChange|> 1) and were further grouped into three distinct expression patterns (Additional file 4: Table S3). Especially, 969 circRNA genes were significantly up- or down-regulated during adipogenic process and were classified as pattern I (AD-specific), represented by PTK2B , TGFBR3 , LEPR  and AEBP1 , whose host genes were key regulators of adipogenesis (Additional file 4: Table S3; Fig. 4). On the other hand, 635 circRNA genes with significantly different expression levels during osteogenesis were classified as pattern II (OS-specific), represented by COL4A , SULF1 , SLC39A8 , IGFBP5 , HOMER2  and ALCAM , whose host genes were involved in osteogenesis (Additional file 4: Table S3; Fig. 4). Besides, 1173 significantly changed circRNA genes in BMSCs were classified as pattern III (BMSC-specific), represented by LPAR1  and EPAS1 , whose host genes were related to the maintenance of MSC homeostasis (Additional file 4: Table S3; Fig. 4). For those representative circRNA genes detected by RiboMinus/RNase R method, their corresponding RNAs detected by RiboMinus and Poly A methods showed similar expressional trends (Fig. 4). Taken together, the results suggested that the circRNAs identified for each pattern may be orchestrated and served as the functional program to regulate MSC differentiation.
Comparison analysis of circRNAs between induced BMSCs and AMSCs
AMSCs are important source for stem cell therapies in addition to BMSCs [45, 46]. The universal clinical application of MSCs from diverse sources requires an understanding of their common mechanisms. Therefore, we established an adipo-osteogenesis system of AMSCs using the same induction medium as the medium used for BMSCs (Additional file 1: Fig. S3). Subsequently, we sequenced 20 samples from differentiated AMSCs at five time points using RiboMinus- and RiboMinus/RNase R-treated methods. We generated an average of 77 million clean reads, of which 92% were precisely mapped to the human genome. BSJ reads were highly abundant in RiboMinus/RNase R data (Fig. 5a; Additional file 1: Fig. S4; Additional file 2: Table S1). After a comparison analysis, 31,326 circRNAs were collectively identified in AMSC- and BMSC-derived datasets, among which 1,166 circRNA genes were significantly changed among different lineages and grouped into AD-specific (613 genes), OS-specific (326 genes) and MSC-specific patterns (168 genes) (Fig. 5b; Additional file 4: Table S3; Additional file 5: Table S4). Gene Ontology functional analysis showed that the host genes in AD-specific pattern were related to cell adhesion , glucose [48, 49] and lipid metabolism [50, 51], which are crucial for adipogenesis. The host genes in OS-specific pattern were involved in cytokine signaling , cell migration  and growth factor signaling , which are important in the regulation of bone formation. The host genes in the MSC-specific pattern were mainly related to cell division related processes , which are important for controlling the pluripotency and differentiation of stem cells (Fig. 5d; Additional file 6: Table S5). Thus, a considerable number of lineage-specific circRNA genes were shared by MSCs from both bone marrow and adipose, and their host genes were involved in distinct functional sets related to differentiation and stemness, which provided the potential common regulators of adipo-osteogenesis processes.
Validation of circRNAs during adipo-osteogenesis of MSCs
To validate the integrative analysis results of multi-strategic RNA-seq data for human BMSCs and AMSCs, we performed qPCR for twelve lineage-specific circRNAs, including osteogenesis-specific (circ-CYP24A1 and circ-SLC39A8, etc.) and adipogenesis-specific circRNAs (circ-SETBP1 and circ-TTC39C, etc.). Divergent primer pairs were designed to amplify regions that included BSJ sites (Additional file 7: Table S6). The GAPDH was used as a reference gene for data normalization. The relative expression results of those twelve circRNAs in induced BMSCs and AMSCs showed similar trends when detected by qPCR or RNA-seq, which indicated that the integrative expression analysis of RNA-seq data was of high quality (Fig. 6; Additional file 1: Fig. S5). Besides, the PCR products were further purified and performed with Sanger sequencing to conform the circRNA sequences. As showed in Fig. 6, the sequencing results precisely matched the BSJ sites of those circRNAs. Therefore, the expression data of circRNA during adipo-osteogenesis of human MSCs were reliable and may be used as a resource for studying the balance between the osteogenesis and adipogenesis of human MSCs.
Circ-CRLF1 was lineage-specially expressed and regulated adipogenesis
To elucidate the roles of lineage-specific circRNAs in adipo-osteogenesis of MSCs, we performed functional validation of the circRNA gene, CRLF1, one of the top ten significantly high-expressed adipogenesis-specific circRNA genes (Fig. 4). Circ-CRLF1 (position, chr19:18596622–18597049; spliced length, 327 nt) was originated from the fifth and sixth exons of the CRLF1 gene (Fig. 7a). The expression levels of CRLF1 mRNA decreased sharply after RNase R treatment, while RNase R failed to degrade circ-CRLF1 (Additional file 1: Fig. S6), which confirmed the stability of circ-CRLF1. To detect circ-CRLF1 and its related mRNA at the same time without interference, divergent primers for circ-CRLF1 were designed to cover the BSJ of the fifth and sixth exons, while convergent primers were designed to cover the exon-spliced junction (ESJ) between the third and fourth exon. The specificity of those primers was successfully validated by Sanger sequencing (Fig. 7a). The qPCR results indicated that the levels of both circ-CRLF1 and its corresponding mRNA were significantly increased during adipogenesis, but not osteogenesis in human BMSCs and AMSCs, which were similar to the RNA-seq results (Fig. 7b). Furthermore, we designed a siRNA targeting the BSJ site to knock down circ-CRLF1. After 7 days of adipogenesis induction, circ-CRLF1 expression levels were significantly reduced by the siRNA in both BMSCs and AMSCs, compared with its expression levels after treatment with scramble siRNAs. Meanwhile, the targeted siRNA did not reduce the expression levels of CRLF1 mRNA (Fig. 7c). The knock-down of circ-CRLF1 significantly increased the expression levels of adipogenesis marker genes, including PPARG, CEBPA and FABP4, in human BMSCs and AMSCs after 7 days of adipogenesis induction (Fig. 7d). Taken together, these results indicated that circ-CRLF1 was lineage-specially expressed during adipo-osteogenesis and regulated the adipogenesis process.
Recent studies have indicated that epigenomic regulators, including ncRNAs, are essential for MSC homeostasis . Compared with other ncRNAs, circRNAs are resistant to degradation due to their circular structure [57, 58]. They have become stable biomarkers and long-acting therapeutic targets of various diseases [59, 60]. The genome-wide expressional studies of circRNAs in osteogenesis or adipogenesis, using single-strategic high-throughput analysis, have been accumulating and helpful for elucidation of the functional roles of circRNAs in bone tissues[17,18,19,20,21]. In the present study, we systematically analyzed the adipo-osteogenesis of both BMSCs and AMSCs using multi-strategic methods. We provided an accurate and generally applicable transcriptomic profile of circRNAs to allow further studies of bone development and dysfunction.
In adipo-osteogenic differentiated MSCs, RiboMinus-treated RNA-seq method generated the highest proportion of coding (87%) and non-coding (68%) transcripts, which were consistent with previous results . Nevertheless, this method could only detect limited numbers of circRNAs, due to the low abundance of circRNAs [24, 62]. RNase R treatment could improve the detection of circRNAs by degrading linear RNAs . We detected 16,767 circRNAs per library in differentiated human BMSCs using the RiboMinus/RNase R-treated strategy. This was at least eight times higher than the number of detected circRNAs using the RiboMinus-treated strategy. However, RNase R treatment may generate bias by preferentially preserving highly-expressed circRNAs [23, 24]. The integrative analysis of RiboMinus data and RiboMinus/ RNase R data could correct RNase R efficiency to improve the quantification of circRNAs . The Poly A method removes most circRNAs during the Poly A selection step, and these data are often used as a blank background . In the present research, we generated the multi-strategic RNA data from RiboMinus, RiboMinus/RNase R, and Poly A libraries for differentiated human BMSCs, and performed integrative analysis using CIRI2 and CIRIquant tools. At last, 48,152 circRNAs from 7366 host genes were identified after systemic correction, and 96.2% of them were overlapped with reference circRNAs, which indicated the high reliability of our data. Similar to previous studies , most of the circRNAs that we identified in differentiated human BMSCs were from protein-coding genes and their exons. Thus, our integrative data provided high-quality circRNA annotations for differentiated human BMSCs.
Public RNA-seq datasets have suggested that circRNAs are dynamically expressed in spatial- and temporal-specific manners, and may have essential roles in developmental processes [62, 65]. Therefore, it is crucial to characterize the expression patterns of circRNAs for studies of bone formation and homeostasis. In the present study, we identified 969 adipogenic, 635 osteogenic and 1,173 BMSC-specific host genes of circRNAs. The top ten up- or down-regulated lineage-specific circRNAs shared similar expressional trends with their corresponding total and linear RNAs. Although the correlation between circRNAs and their related linear RNAs have differed from various studies, this may be due to the different selected clusters of circRNAs [66, 67]. Moreover, several host genes of those circRNAs with differential expression levels were involved in lineage-specific activities [33,34,35,36,37,38,39,40,41,42,43,44]. Thus, these lineage-specific circRNAs may provide regulatory candidates of adipo-osteogenesis of human BMSCs.
MSCs are heterogenous populations and differ depending on their tissue origin or the conditions of the donors (e.g., age, diseases, or unknown factors) . Therefore, it is urgent to identify the common regulatory factors of MSCs for their universal application. In the present study, we annotated 31,326 shared circRNAs and identified 613 adipogenic and 326 osteogenic host genes with significantly altered expression levels in both BMSCs and AMSCs. The shared up-regulated host genes in each lineage of MSCs were related to distinct biological processes, including cell adhesion, cytokine signaling, and cell division, which indicated their potential involvement in the regulation of adipo-osteogenesis of MSCs. In particular, the circRNA of CRLF1 was significantly up-regulated during the adipogenesis of both BMSCs and AMSCs. The knock-down of circ-CRLF1 significantly promoted adipogenesis, which suggested that circ-CRLF1 could be a potential target for the prevention or treatment of bone-related diseases induced by an imbalance of adipo-osteogenesis.
In summary, we have annotated over ten thousand of circRNAs, and identified various circRNAs that were lineage-specifically regulated during osteogenesis and adipogenesis of human BMSCs and AMSCs. Their host genes were involved in distinct biological processes, including cell adhesion, cytokine signaling, and cell division. Furthermore, we identified a lineage-specific circRNA, circ-CRLF, that negatively regulates adipogenesis. This study paves the way for further investigation to understand the potential roles of circRNAs in bone development and identify the therapeutic targets of bone diseases.
Availability of data and materials
All the datasets supporting the results have been listed in the article and its additional files. The RNA-seq data have been deposited to the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) under accession code PRJNA688535.
Mesenchymal stem cells
Bone marrow-derived mesenchymal stem cells
Adipose-derived mesenchymal stem cells
- RiboMinus/RNase R:
Ribosome RNAs-depleted and RNase R-treated
Uccelli A, Moretta L, Pistoia V. Mesenchymal stem cells in health and disease. Nat Rev Immunol. 2008;8(9):726–36.
Infante A, Rodríguez CI. Osteogenesis and aging: lessons from mesenchymal stem cells. Stem Cell Res Ther. 2018;9(1):244.
Chen Q, Shou P, Zheng C, Jiang M, Cao G, Yang Q, et al. Fate decision of mesenchymal stem cells: adipocytes or osteoblasts? Cell Death Differ. 2016;23(7):1128–39.
Muruganandan S, Govindarajan R, Sinal CJ. Bone marrow adipose tissue and skeletal health. Curr Osteoporos Rep. 2018;16(4):434–42.
Chen LL. The expanding regulatory mechanisms and cellular functions of circular RNAs. Nat Rev Mol Cell Biol. 2020;21(8):475–90.
Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.
Kristensen LS, Andersen MS, Stagsted LVW, Ebbesen KK, Hansen TB, Kjems J. The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet. 2019;20(11):675–91.
Haddad G, Lorenzen JM. Biogenesis and function of circular RNAs in health and in disease. Front Pharmacol. 2019;10:428.
Abdelmohsen K, Panda AC, Munk R, Grammatikakis I, Dudekula DB, De S, et al. Identification of HuR target circular RNAs uncovers suppression of PABPN1 translation by CircPABPN1. RNA Biol. 2017;14(3):361–9.
Du WW, Yang W, Liu E, Yang Z, Dhaliwal P, Yang BB. Foxo3 circular RNA retards cell cycle progression via forming ternary complexes with p21 and CDK2. Nucleic Acids Res. 2016;44(6):2846–58.
Chen N, Zhao G, Yan X, Lv Z, Yin H, Zhang S, et al. A novel FLI1 exonic circular RNA promotes metastasis in breast cancer by coordinately regulating TET1 and DNMT1. Genome Biol. 2018;19(1):218.
Chen G, Wang Q, Li Z, Yang Q, Liu Y, Du Z et al. Circular RNA CDR1as promotes adipogenic and suppresses osteogenic differentiation of BMSCs in steroid-induced osteonecrosis of the femoral head. Bone. 2020; 133:115258.
Ouyang Z, Tan T, Zhang X, Wan J, Zhou Y, Jiang G, et al. CircRNA hsa_circ_0074834 promotes the osteogenesis-angiogenesis coupling process in bone mesenchymal stem cells (BMSCs) by acting as a ceRNA for miR-942-5p. Cell Death Dis. 2019;10(12):932.
Zhang D, Ni N, Wang Y, Tang Z, Gao H, Ju Y, et al. CircRNA-vgll3 promotes osteogenic differentiation of adipose-derived mesenchymal stem cells via modulating miRNA-dependent integrin α5 expression. Cell Death Differ. 2021;28(1):283–302.
Jiang R, Li H, Yang J, Shen X, Song C, Yang Z, et al. circRNA Profiling reveals an abundant circFUT10 that Promotes adipocyte proliferation and inhibits adipocyte differentiation via sponging let-7. Mol Ther Nucleic Acids. 2020;20:491–501.
Zhu Y, Gui W, Lin X, Li H. Knock-down of circular RNA H19 induces human adipose-derived stem cells adipogenic differentiation via a mechanism involving the polypyrimidine tract-binding protein 1. Exp Cell Res. 2020; 387(2):111753.
Della Bella E, Menzel U, Basoli V, Tourbier C, Alini M, Stoddart MJ. Differential regulation of circRNA, miRNA, and piRNA during early osteogenic and chondrogenic differentiation of human mesenchymal stromal cells. Cells. 2020; 9(2).
Zhang M, Jia L, Zheng Y. circRNA expression profiles in human bone marrow stem cells undergoing osteoblast differentiation. Stem Cell Rev Rep. 2019;15(1):126–38.
Li Z, Li N, Ge X, Pan Y, Lu J, Gobin R, et al. Differential circular RNA expression profiling during osteogenic differentiation of stem cells from apical papilla. Epigenomics. 2019;11(9):1057–73.
Zheng Y, Li X, Huang Y, Jia L, Li W. The circular RNA landscape of periodontal ligament stem cells during osteogenesis. J Periodontol. 2017;88(9):906–14.
Arcinas C, Tan W, Fang W, Desai TP, Teh DCS, Degirmenci U, et al. Adipose circular RNAs exhibit dynamic regulation in obesity and functional role in adipogenesis. Nat Metab. 2019;1(7):688–703.
Xiao MS, Wilusz JE. An improved method for circular RNA purification using RNase R that efficiently removes linear RNAs containing G-quadruplexes or structured 3’ ends. Nucleic Acids Res. 2019;47(16):8755–69.
Li X, Zhang B, Li F, Yu K, Bai Y. The mechanism and detection of alternative splicing events in circular RNAs. PeerJ. 2020; 8:e10032.
Zhang J, Chen S, Yang J, Zhao F. Accurate quantification of circular RNAs identifies extensive circular isoform switching events. Nat Commun. 2020;11(1):90.
Rauch A, Haakonsson AK, Madsen JGS, Larsen M, Forss I, Madsen MR, et al. Osteogenesis depends on commissioning of a network of stem cell transcription factors that act as repressors of adipogenesis. Nat Genet. 2019;51(4):716–27.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.
Frankish A, Diekhans M, Ferreira AM, Johnson R, Jungreis I, Loveland J, et al. GENCODE reference annotation for the human and mouse genomes. Nucleic Acids Res. 2019;47(D1):D766-d773.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, et al. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8(9):R183.
Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Brief Bioinform. 2018;19(5):803–10.
Wu W, Ji P, Zhao F. CircAtlas: an integrated resource of one million highly accurate circular RNAs from 1070 vertebrate transcriptomes. Genome Biol. 2020;21(1):101.
Ni B, Farrar JS, Chen S, Lownik JC, Celi FS. A novel role for PTK2B in cultured beige adipocyte differentiation. Biochem Biophys Res Commun. 2018;501(4):851–7.
Lee MJ, Pickering RT, Shibad V, Wu Y, Karastergiou K, Jager M, et al. Impaired glucocorticoid suppression of TGFβ signaling in human omental adipose tissues limits adipogenesis and may promote fibrosis. Diabetes. 2019;68(3):587–97.
Palhinha L, Liechocki S, Hottz ED, Pereira J, de Almeida CJ, Moraes-Vieira PMM, et al. Leptin induces proadipogenic and proinflammatory signaling in adipocytes. Front Endocrinol (Lausanne). 2019;10:841.
Ro HS, Zhang L, Majdalawieh A, Kim SW, Wu X, Lyons PJ, et al. Adipocyte enhancer-binding protein 1 modulates adiposity and energy homeostasis. Obesity (Silver Spring). 2007;15(2):288–302.
Dussold C, Gerber C, White S, Wang X, Qi L, Francis C, et al. DMP1 prevents osteocyte alterations, FGF23 elevation and left ventricular hypertrophy in mice with chronic kidney disease. Bone Res. 2019;7:12.
Zaman G, Staines KA, Farquharson C, Newton PT, Dudhia J, Chenu C, et al. Expression of Sulf1 and Sulf2 in cartilage, bone and endochondral fracture healing. Histochem Cell Biol. 2016;145(1):67–79.
Song J, Kim D, Lee CH, Lee MS, Chun CH, Jin EJ. MicroRNA-488 regulates zinc transporter SLC39A8/ZIP8 during pathogenesis of osteoarthritis. J Biomed Sci. 2013;20(1):31.
Mukherjee A, Rotwein P. Insulin-like growth factor-binding protein-5 inhibits osteoblast differentiation and skeletal growth by blocking insulin-like growth factor actions. Mol Endocrinol. 2008;22(5):1238–50.
Son A, Kang N, Oh SY, Kim KW, Muallem S, Yang YM, et al. Homer2 and Homer3 modulate RANKL-induced NFATc1 signaling in osteoclastogenesis and bone metabolism. J Endocrinol. 2019;242(3):241–9.
Bruder SP, Ricalton NS, Boynton RE, Connolly TJ, Jaiswal N, Zaia J, et al. Mesenchymal stem cell surface antigen SB-10 corresponds to activated leukocyte cell adhesion molecule and is involved in osteogenic differentiation. J Bone Miner Res. 1998;13(4):655–63.
Li M, Lv Y, Chen F, Wang X, Zhu J, Li H, et al. Co-stimulation of LPAR(1) and S1PR(1/3) increases the transplantation efficacy of human mesenchymal stem cells in drug-induced and alcoholic liver diseases. Stem Cell Res Ther. 2018;9(1):161.
Zhu C, Yu J, Pan Q, Yang J, Hao G, Wang Y, et al. Hypoxia-inducible factor-2 alpha promotes the proliferation of human placenta-derived mesenchymal stem cells through the MAPK/ERK signaling pathway. Sci Rep. 2016;6:35489.
El-Badawy A, Amer M, Abdelbaset R, Sherif SN, Abo-Elela M, Ghallab YH, et al. Adipose stem cells display higher regenerative capacities and more adaptable electro-kinetic properties compared to bone marrow-derived mesenchymal stromal cells. Sci Rep. 2016;6:37801.
Mazini L, Rochette L, Amine M, Malka G. Regenerative capacity of adipose derived stem cells (ADSCs), comparison with mesenchymal stem cells (MSCs). Int J Mol Sci. 2019; 20(10).
Frith JE, Mills RJ, Hudson JE, Cooper-White JJ. Tailored integrin-extracellular matrix interactions to direct human mesenchymal stem cell differentiation. Stem Cells Dev. 2012;21(13):2442–56.
Geoghegan G, Simcox J, Seldin MM, Parnell TJ, Stubben C, Just S, et al. Targeted deletion of Tcf7l2 in adipocytes promotes adipocyte hypertrophy and impaired glucose metabolism. Mol Metab. 2019;24:44–63.
Kolodziej M, Strauss S, Lazaridis A, Bucan V, Kuhbier JW, Vogt PM, et al. Influence of glucose and insulin in human adipogenic differentiation models with adipose-derived stem cells. Adipocyte. 2019;8(1):254–64.
Sun S, Adyshev D, Dudek S, Paul A, McColloch A, Cho M. Cholesterol-dependent modulation of stem cell biomechanics: application to adipogenesis. J Biomech Eng. 2019;141(8):0810051–08100510.
Baugh JE Jr, Floyd ZE, Stephens JM. The modulation of STAT5A/GR complexes during fat cell differentiation and in mature adipocytes. Obesity (Silver Spring). 2007;15(3):583–90.
Sims NA. The JAK1/STAT3/SOCS3 axis in bone development, physiology, and pathology. Exp Mol Med. 2020;52(8):1185–97.
Su P, Tian Y, Yang C, Ma X, Wang X, Pei J et al. Mesenchymal stem cell migration during bone formation and bone diseases therapy. Int J Mol Sci. 2018; 19(8).
Chen J, Yuan K, Mao X, Miano JM, Wu H, Chen Y. Serum response factor regulates bone formation via IGF-1 and Runx2 signals. J Bone Miner Res. 2012;27(8):1659–68.
Soufi A, Dalton S. Cycling through developmental decisions: how cell cycle dynamics control pluripotency, differentiation and reprogramming. Development. 2016;143(23):4301–11.
Sui BD, Zheng CX, Li M, Jin Y, Hu CH. Epigenetic regulation of mesenchymal stem cell homeostasis. Trends Cell Biol. 2020;30(2):97–116.
Enuka Y, Lauriola M, Feldman ME, Sas-Chen A, Ulitsky I, Yarden Y. Circular RNAs are long-lived and display only minimal early alterations in response to a growth factor. Nucleic Acids Res. 2016;44(3):1370–83.
Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, et al. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19(2):141–57.
Lei B, Tian Z, Fan W, Ni B. Circular RNA: a novel biomarker and therapeutic target for human cancers. Int J Med Sci. 2019;16(2):292–301.
Li Y, Zheng Q, Bao C, Li S, Guo W, Zhao J, et al. Circular RNA is enriched and stable in exosomes: a promising biomarker for cancer diagnosis. Cell Res. 2015;25(8):981–4.
Cui P, Lin Q, Ding F, Xin C, Gong W, Zhang L, et al. A comparison between ribo-minus RNA-sequencing and polyA-selected RNA-sequencing. Genomics. 2010;96(5):259–65.
Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, et al. Circular RNAs in the Mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell. 2015;58(5):870–85.
Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):4.
Starke S, Jost I, Rossbach O, Schneider T, Schreiner S, Hung LH, et al. Exon circularization requires canonical splice signals. Cell Rep. 2015;10(1):103–11.
Mahmoudi E, Cairns MJ. Circular RNAs are temporospatially regulated throughout development and ageing in the rat. Sci Rep. 2019;9(1):2564.
Ruan H, Xiang Y, Ko J, Li S, Jing Y, Zhu X, et al. Comprehensive characterization of circular RNAs in ~1000 human cancer cell lines. Genome Med. 2019;11(1):55.
Zhang M, Xu K, Fu L, Wang Q, Chang Z, Zou H et al. Revealing epigenetic factors of circRNA expression by machine learning in various cellular contexts. iScience. 2020; 23(12):101842.
Costa LA, Eiro N, Fraile M, Gonzalez LO, Saá J, Garcia-Portabella P et al. Functional heterogeneity of mesenchymal stem cells from natural niches to culture conditions: implications for further clinical uses. Cell Mol Life Sci. 2020:1–21.
We thank for the technical support by the core facilities of Translational Medicine Collaborative Innovation Center, Shenzhen People's Hospital.
This work was supported by the China Postdoctoral Science Foundation (2018M640885), the Science and Technology Project of Shenzhen (ZDSYS20190902093203727) and the Science and Technology Project of Shenzhen Municipality (JCYJ20190807153413130).
Ethics approval and consent to participate
All aspects of the cell experimental procedures were approved by the ethics committee of Jinan University.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1
. Figure S1: Statistical analysis of human BMSCs. Figure S2. General features of multi-strategic libraries for adipo-osteogenesis of BMSCs. Figure S3. Adipogenesis and osteogenesis of human AMSCs. Figure S4. General features of multi-strategic libraries for adipo-osteogenesis AMSCs. Figure S5. Confirmation of circRNAs by qPCR and Sanger sequencing during adipo-osteogenesis of human AMSCs. Figure S6. The abundant of circ-CRLF1 and CRLF1 mRNA at day 14 during adipogenesis treated with RNase R.
Additional file 2
. Table S1: Detailed information for multi-strategic RNA-seq of adipo-osteogenesis of MSCs.
Additional file 3
. Table S2: Identified circRNAs through integrative analysis of adipo-osteogenesis of BMSCs.
Additional file 4
. Table S3: Significantly differentially expressed circRNA genes in induced BMSCs from RiboMinus/RNase R data.
Additional file 5
. Table S4: Significantly differentially expressed circRNA genes in induced AMSCs from RiboMinus/RNase R data.
Additional file 6
. Table S5: Top ten significantly enriched GO terms for lineage-specific host genes of circRNAs.
Additional file 7
. Table S6: Oligonucleotides used in this study.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Huang, HB., Luo, HT., Wei, NN. et al. Integrative analysis reveals a lineage-specific circular RNA landscape for adipo-osteogenesis of human mesenchymal stem cells. Stem Cell Res Ther 13, 106 (2022). https://doi.org/10.1186/s13287-022-02792-5
- Mesenchymal stem cells
- Circular RNAs