CDR1as related miRNA-mRNA networks in differentiating goat skeletal muscle satellite cells

Background: Myogenesis is a complex process controlled by several coding and non-coding RNAs (ncRNAs) such as circular RNAs (circRNAs) that well-known function as endogenous microRNAs (miRNAs) sponges. Over the past few years, numerous circRNAs have been known and their roles in biological processes have begun to be understood. Cerebellar Degeneration-Related protein 1 antisense (CDR1as), the most spotlighted circRNA as miR-7 sponge that has been blooming circRNAs’ research for a decade, and can potentially sponge several miRNAs in disease and muscle physiology. Nevertheless, the linear-RNAs-differed character that the acute interventions for circRNAs do not affect miRNAs levels, and has retarded the transcriptome-wide discovery of miRNAs sponged by. Therefore, the purpose of this study was to provide the transcriptomic effect of CDR1as during muscle differentiation. Methods: siCDR1as and siDICER1 were transfected into goat skeletal muscle satellite cells (SMSCs). RNA-seq technology and bioinformatics tools were used to analyze genes that are deregulated by siCDR1as and siDICER1. quantitative PCR was used to verify the expression levels of the differentially expressed mRNAs and miRNAs. Results: Here, to systematically identify miRNAs targeting CDR1as, we employed the critical enzyme DICER1 that governs the biogenesis of miRNAs. The deciency of either DICER1 or CDR1as inhibited myogenic differentiation of SMSCs, and knockdown of DICER1 decreased the expression of CDR1as. Moreover, we screened for the targeted messenger RNAs (mRNAs) and miRNAs in SMSCs transfected with siDICER1 or siCDR1as respectively and found out that some well-known muscle-related pathways such as phosphoinositide 3-kinase (PI3K)-AKT signaling pathway, Rap1 signaling pathway, and MAPK signaling pathway were enriched in all groups. Further, regarding the miRNAs identied in siDICER1 and siCDR1as together with the sequence complementary information, we identied 11 miRNAs including miR-1, GAPDH was used as an internal control. KEGG of the downregulated mRNAs with the top 20 enrichment. Bubble color and size correspond to the Q value and gene number enriched in the pathway. The rich factor indicates the ratio of the number of DEGs mapped to a certain pathway to the total number of genes mapped to this pathway.


Cell Transfection
For gain and loss function study, SMSCs with 80%-90% con uency in 6-well plates were transfected with Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) with siCDR1as, siDICER1 or siNC, at a concentration of 50 nM, according to manufacturer's instructions. Moreover, pCDNA-3 was also transfected as a negative control for CDR1as. After 5 hr of transfection, the medium containing transfection reagent was replaced by fresh GM (Growth medium). After culturing for an additional 48 hr, cells were used for further experiments. The siCDR1as and siDICER1 sequence are 5'-GTCTACGATATCCAGGGTT-3' and 5'-GCAGTTACGATTTAGCTAA-3' respectively.

Rna Extraction
Total RNAs were extracted from cells with liquid nitrogen by using RNAiso Plus reagent (TaKaRa Bio, Inc., Japan) and puri ed using a QIAGEN RNeasy Mini Kit (QIAGEN, Chatsworth, CA, USA) according to the manufacturer's instructions. RNA degradation and contamination were monitored on 1.5% agarose gels.
The purity and concentration of the RNAs were measured using NanoPhotometer® spectrophotometer (IMPLEN, Los Angeles, CA, USA) and a Qubit RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), respectively. All RNA samples were stored at -80℃ until further use.

Mirna Sequence
After verifying the concentration and purity, the integrity was assessed using an RNA Nano 6000 Assay Kit in a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Only samples that had RNA Integrity Number (RIN) scores > 7.5 were used for RNA sequencing. Brie y, small RNAs were reversed transcribed and ampli ed by PCR. The PCR products were then puri ed by denaturing polyacrylamide gel electrophoresis (PAGE). A total of 3 µg RNA per sample was used as input material for miRNA sampling preparation. miRNA libraries were constructed and sequencing was performed on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA), and 125-bp long paired-end reads were generated. miRDeep2 software was used to predict the novel miRNAs with trimmed reads. Then the reads were aligned to merged pre-miRNA databases (known pre-miRNA from miRBase v21 plus the newly predicted pre-miRNAs) using Novoalign software (v2.07.11) with, at most, one mismatch. We used the most abundant isomiR, the mature miRNA annotated in miRBase and all isoforms of miRNA (5p or 3p) to calculate miRNA expression. Fold change and p-value were used to calculate the differentially expressed miRNA pro les between two groups. Hierarchical clustering was performed to generate an overview of the characteristics of expression pro les based on values of a signi cant differentially expressed transcripts.

Mrna Sequencing Data Processing
Clean reads were obtained by removing reads containing adapters, reads containing over 10% of poly (N), and low quality reads (> 50% of the bases had Phred quality scores ≤ 10) from the raw data. All downstream and upstream analyses were based on high-quality clean data. Goat reference genome and gene model annotation les were downloaded from NCBI database (CHIR_1.0, NCBI) [33]. Index of the reference genome was built using Bowtie v2.0.6 [34,35] and paired-end clean reads were aligned to the reference genome using TopHat v2.0.14 [36]. The mapped reads from each library were assembled with cu inks v2.2.1 [37]. We used the reference annotation based transcript (RABT) assembly method in cu inks v2.2.1 to construct and identify mRNA transcripts from the TopHat2 alignment results. Hierarchical clustering was performed to generate an overview of the characteristics of expression pro les based on values of a signi cant differentially expressed transcripts.

Kegg Pathway Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs was performed with KOBAS software [38] using a hypergeometric test. KEGG pathways with Q value < 0.05 were considered signi cantly enriched.

Validation Of Rna-seq Data By Qpcr
Total RNA was extracted from SMSCs transfected with siCDR1as siNC using RNAiso Plus reagent (TaKaRa Bio, Inc., Japan). For qPCR of mRNA, all PCR primers were designed at or just outside exon/exon junctions to avoid the ampli cation of residual genomic DNA using the Primer-BLAST on the NCBI website, and speci city was determined using BLASTN. 1 µg of total RNAs (1 mg) was reversetranscribed into cDNA by using PrimeScriptTM RT reagent Kit with gDNA Eraser (TaKaRa, Otsu, Japan). And using these cDNA as templates, expression levels of genes were quanti ed by quantitative real time PCR (RT_qPCR) in a Bio-Rad CFX96 system (Bio-Rad, Hercules, USA) with SYBR Premix Ex TaqTM II (TaKaRa, Otsu, Japan), according to the manufacturer's protocols. Three samples were collected for each treatment and each sample at least triplicates. The PCR protocol was as follows: denaturation at 95 ℃ for 30 s, followed by 40 cycles of 95 ℃ for 20 s and 60 ℃ for 20 s, then 72 ℃ for 30 s. The 2 −△△Ct procedure was used to calculate the relative expression levels of mRNAs, with GAPDH as an internal control [39]. For qPCR of miRNA, reversed transcribed using The First-strand cDNA Synthesize (TaKaRa, Mount view, CA, USA). For real-time PCR, all reactions were performed in triplicate with SYBR Premix Ex TaqTM II (TaKaRa, Otsu, Japan) under the following conditions: 10 s at 95℃ for initial denaturation, followed by 39 cycles of 95 ℃ for 5 s and 60 ℃ for 20 s, then melting curve of 65 ℃ to 95℃ for 5 s. The expression levels of U6 were used to normalize the expression levels of the gene of interest. Primers for the mRNAs and miRNAs were designed using Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primers-blast) (Table S1 a, b, c, and d).

Statistical analysis
Data are expressed as mean ± SEM. All statistical analyses will be performed using GraphPad Prism 6.01.
Unpaired student's tests were used to performed statistical analysis used for two group comparisons. Statistical signi cant difference was evaluated p < 0.05.

Results
Effect of siCDR1as and siDICER1 on SMSCs differentiation and the relationship between CDR1as and DICER1 DICER1 is central to miRNA-mediated silencing [40] and modulates miRNAs during myogenic differentiation [18]. Meanwhile, CDR1as plays a critical role in myogenesis by functioning as a molecular sponge for miR-7 [29]. To identify whether CDR1as and DICER1 are associated in myoblast differentiation, we in vitro cultured SMSCs isolated from Longissimus dorsi (LD) muscles of newborn Jiangzhou Big-Eared goats and knocked down the expression of CDR1as and DICER1 (Fig. 1a, d) as well as mouse myoblast C2C12 (Suppl. Figure 1, 2) by using their respective interfering RNAs. Similar to previous results, de ciency of CDR1as or DICER1 decreased the expression of MyoD mRNA, the master myogenic gene (Fig. 1b, e), and the formation of myotubes (Fig. 1c, h). Moreover, CDR1as was shown to be downregulated in SMSCs transfected with siDICER1 ( Fig. 1f) and co-transfection of pCDNA3 (negative control for CDR1as) and siDICER1 (Fig. 1g), suggesting that the close relationship between DICER1 and CDR1as could help identify miRNAs and their targeted mRNAs regulated by CDR1as. Overview of mRNA and miRNA sequencing data associated with DICER1 and CDR1as In order to predict mRNAs and miRNAs associated with DICER1 and CDR1as in SMSCs, we systematically cultured SMSCs and knocked down the expression of CDR1as and DICER1, with three biological replicates for each treatment. Cells were harvested at 48 h after transfection and the total RNA were extracted to construct the cDNA libraries individually for mRNA-seq and miRNA-seq using an Illumina HiSeq 2500 platform and 125 bp paired-end reads.
After removing low-quality sequences and adapters, considering the siCDR1as results, an average of 59,821,102 and 58,552,463 mRNAs were produced from raw and clean reads respectively (Table 1). In addition, an average of 15,350,087 raw reads and 14,614,375 clean reads miRNAs for siCDR1as were generated respectively ( Table 2). On the other hand, siDICER1 has an average of 128,604,571 and 120,968,505 mRNAs were obtained from raw and clean reads respectively (Table 3). Also, an average of 14,260,428 raw reads and 14,028,682 clean reads miRNAs for siDICER1 were acquired accordingly (Table 4).
To explore the expression relationship of genes between samples, the Pearson's correlation coe cient (PCC) of mRNAs and miRNAs expression levels of siCDR1as-1, 2 and 3, as well as, siDICER1-1, 2, and 3 in SMSCs were calculated and used to generate a correlation chart. As shown, the correlation coe cient of the siCDR1as-1, 2, and 3 as well as, siDICER1; 1, 2, and 3 in SMSCs ranged from 0.96 to 0.99 (average of 0.98), indicating that the samples replicate very well biologically (Fig. 2a, b, c, and d). The abscissa and the ordinate were the respective samples, and the abscissa and the ordinate of each patch represented the correlation of siCDR1as and siDICER1 samples. Importantly, two completely related genomes had a value of 1. The closer to 1 the relative value is, the larger the Pearson correlation coe cient (PCC) for the siCDR1as and siDICER1 samples; conversely, the closer to 0 the relative value was, the smaller the PCC between the siCDR1as and siDICER1 samples.
To verify the RNA-sequencing data, in extension samples of SMSCs with siDICER1 or siCDR1as, ve downregulated as well as ve upregulated mRNA and miRNAs identi ed were randomly selected and quanti ed for their expression levels using qRT-PCR. As shown in (Fig. 3b; Suppl.  Table 2d), all the differential expression tendency were con rmed.
In addition, KEGG pathway analysis was conducted based on the differentially expressed mRNAs (Q value < 0.05). The results showed that among the top 20 most enriched KEGG pathways of the downregulatory genes include PI3K-AKT signaling pathway, Focal adhesion, and ECM receptor interaction (Fig. 3c), which are all associated with muscle development. Moreover, considering the upregulatory genes (Suppl. Figure 3), none of the 20 enriched KEGG pathways was associated with muscle development. These data show that the DICER1 gene is involved in the formation of SMSCs.
Furthermore, using TargetScan, RNAhybrid, and miRanda, a total of 7,194 targets mRNAs were predicted as targets of those upregulated miRNA genes whilst 1,686 target mRNAs of the downregulatory miRNAs were screened. All the mRNAs for the upregulated-miRNAs were enriched in the KEGG pathway such as, Focal adhesion, FoxO signaling pathway, MAPK signaling pathway, PI3K-AKT signaling pathway, Rap1 signaling pathway, and mTOR signaling pathway which are related to muscle development ( Fig. 4c; Suppl. Figure 4). Moreover, this result is consistent with the functional enrichment for DE mRNAs mentioned before. Bubble color and size correspond to the Q value and gene number enriched in the pathway. The rich factor indicates the ratio of the number of DEGs mapped to a certain pathway to the total number of genes mapped to this pathway.
Further, we performed KEGG pathway analysis with KOBAS software based on the differentially expressed mRNAs genes (Q value < 0.05). The results showed that among the top 20 most enriched pathways of the downregulated genes in siCDR1as samples, some well-known muscle-related pathways including PI3K-AKT signaling pathway signaling pathway, Focal adhesion, Rap1 signaling pathway, and MAPK signaling pathway (Fig. 5c) were identi ed. While those genes up-regulated by siCDR1as were mainly enriched in fatty acid biogenesis and metabolism (Suppl. Figure 5). These results indicate that just like DICER1, the downregulated mRNAs caused by the de cit of CDR1as are closely related to muscle development.
Further, we performed KEGG pathway analysis with KOBAS software based on the differentially expressed mRNAs genes (Q value < 0.05). The results showed that among the top 20 most enriched pathways of the downregulated genes in siCDR1as samples, some well-known muscle-related pathways including PI3K-AKT signaling pathway signaling pathway, Focal adhesion, Rap1 signaling pathway, and MAPK signaling pathway (Fig. 5c) were identi ed. While those genes up-regulated by siCDR1as were mainly enriched in fatty acid biogenesis and metabolism (Suppl. Figure 5). These results indicate that just like DICER1, the downregulated mRNAs caused by the de cit of CDR1as are closely related to muscle development. Anchoring The Novel Core Mirnas Regulated By Cdr1as Considering the major function of DICER1 on modulating miRNAs during myogenesis, and its closely positive effect on CDR1as as shown by expression of CDR1as as well as the enrichment results. First of all, we investigated the critical miRNAs that mediate the function of DICER1, using online tool MSigDB (http://www.broadinstitute.org/gsea/msigdb/index.jsp) and the DE mRNAs dataset caused by siDICER1. We found out 100 miRNAs were potentially targeted by these DE mRNA genes altered by siDICER1 (FDR q-value < 0.004), among which 11 miRNAs including miR-1, miR-199b-5p, miR-206, miR-27a-5p, miR-19b-3p, miR-30b-5p, miR-129-5p, miR-128-3p, miR-30e-5p, miR-27b-5p, and miR-424-5p were overlapped with those differentially expressed miRNAs detected by using miRNA-seq in SMSCs transfected with siDICER1.
This indicates that these miRNAs may play important roles in mediating DICER's function in myogenic differentiation of SMSCs.

Discussion
Recent studies indicate that muscle development is associated with several coding and non-coding genes, including mRNAs, circular RNAs, and miRNAs [22]. For instance, circFGFR4 induces myoblast differentiation through the upregulation of Wnt3a via the downregulation of miR-107 [56]. Another example includes the inhibition of miR-203 by circSVIL via increasing expression of MEF2C during muscle differentiation [49]. DICER1 is known to play a signi cant role in vascular smooth muscle development and function by controlling proliferation and contractile differentiation [30]. Therefore studying the transcriptome view of CDR1as and DICER1 in muscle development as well as their relationship during muscle development in SMSCs may contribute to the understanding of genes involved in myogenesis. In this study, we found that CDR1as and DICER1 contribute signi cantly to myogenesis through differentiation and are positively correlated.
Based on the properties of coding and non-coding RNAs, various RNAs are known to be related to cancerous diseases [57,58]. For instance, circRIP2 promotes bladder cancer progression through miR-1305/Tfg-β2/smad3 pathway [59]. Also, miR-135a regulates the proliferation and chemosensitivity of endometrial cancer cells by targeting AKT signaling pathway [60]. However, the number of miRNAs and mRNAs associated with CDR1as and DICER1 in SMSCs remains unknown. In our studies, we knockdown the expression of CDR1as and DICER1 in SMSCs, our interest in the results were based on the downregulated mRNAs and the upregulated miRNAs. A total of 789 mRNAs were differentially expressed, comprising of 401 downregulated mRNAs and 388 upregulated mRNAs. In addition, 27 miRNAs were upregulated, with 16 miRNAs being downregulated after the knockdown of CDR1as in SMSCs. Some of the signi cantly downregulated mRNAs such as MEF2C [6], ANGPT1 [45], E2F2 [46], CCN1 [47], and FGFR1 [48], are shown to regulate muscle development. Besides, after the knockdown of DICER1, there was a total of 1,113 differentially expressed mRNAs consisting of 686 downregulated mRNAs and 427 upregulated mRNA. Also, there were 15 upregulated miRNAs and 7 downregulated miRNAs. In addition, among the signi cantly downregulated genes such as ANGPT1 [45], MEF2D [41], IGFBP5 [61], E2F2 [46], and CCN1 [47] are also known to be associated with muscle development.
Skeletal muscle development is a gradual process that involves proliferation, differentiation, and fusion.
Myoblast differentiation is characterized as the key stage responsible for myogenesis. Moreover, myogenesis is controlled by various signaling regulatory networks, including PI3K-AKT signaling pathway, Focal adhesion, Rap1 signaling pathway, MAPK signaling pathway, and FoxO signaling. For instance, in our study, the KEGG pathway was performed to study the putative functions of DEGs. In siCDR1as pro les, the most enriched pathways of the downregulated mRNAs were (PI3K)-AKT signaling pathway, Focal adhesion, Rap1 signaling pathway, and MAPK signaling pathway. The downregulation of PI3K/AKT signaling pathway by miR-106a-5p was shown to inhibit myogenic differentiation of C2C12 myoblast [62]. Also, according to Jiang et al. the knockdown of PI3-Kinase or its downstream target AKT prevents muscle differentiation in cell culture, while activation of PI3-Kinase and Akt induces myogenic differentiation [63,64]. Previous studies have indicated that focal adhesion together with ECM are signaling centers of several intracellular pathways related to cell proliferation and differentiation [65].
Also, focal adhesion plays a signi cant regulatory role in the biological process including muscle differentiation and striated muscle tissue development [66]. Interestingly, Rap 1 signaling is known to be related beta-adrenergic signaling pathway which has shown to play a crucial role in skeletal muscle growth and development [67]. Moreover, Rap 1 is a small GTPase that regulates different processes, including tightening of cell-cell junction [68], cell polarity, and cell adhesion [69]. MAPK signaling pathway is indicated to regulate the regeneration and proliferation of muscle stem cells [70]. In addition, the MAPK signaling pathway plays a signi cant role during myoblast differentiation [71,72]. FoxO signaling pathway is an essential pathway related to muscle development at different stages [66]. Aside from the above-mentioned signaling pathways, the KEGG pathways of the upregulated miRNAs also include WNT and Hippo signaling pathways which are also known to be among the signi cantly enriched pathways (Fig. 4A). Besides, WNT signaling pathway is involved in the control of satellite cell differentiation and regeneration [73,74], while Hippo signaling pathway is also an essential pathway known to be involved in myogenesis [75,76].
Considering the knockdown of DICER1, we realized that the KEGG pathways of the downregulated mRNAs include PI3K-AKT signaling pathway [62], Focal adhesion [66], and ECM receptor interaction [66], which are also recognized to be associated with muscle development. Moreover, FoxO signaling pathway [66], MAPK signaling pathway [70], Rap1 signaling pathway [68,69], and mTOR signaling pathway [77] were also known to be associated with the upregulated miRNAs. Moreover, all these pathways are related to muscle development.
According to Xie et al. six elevated miRNAs including miR-1a-3p, miR-133a-3p, miR-133b-3p, miR-206-3p, miR-128-3p, and miR-351-5p were being classi ed as myogenesis-associated miRNAs (mamiRs), negatively regulated the JNK/MAPK pathway by inhibiting several factors related to the phosphorylation of the JNK/MAPK pathway in skeletal muscle differentiation [78]. miR-146a-5p has been demonstrated to signi cantly induced VSMC proliferation [55], however, lack of miR-146a did not increase muscular dystrophy in mdx mice but was highly expressed in dystrophic muscle [79]. Decreased in the expression of miR-199b-5p induces differentiation of Bone-Marrow Mesenchymal Stem Cells (BMSCs) toward Cardiomyocyte-Like Cells [51]. The Knockdown of miR-199a-3p inhibits slow-to-fast muscle ber type change in mice and C2C12 Cells [80]. Overexpression of miR-19b-3p prevents VSMC proliferation by targeting CTGF expression [81]. miR-30 family is shown to increase myoblast differentiation in vitro and provide negative feedback on the miRNA pathway [82]. However, according to Zhang et al. miR-30-5p has a negative effect on the differentiation of C2C12 cells by targeting MBNL [50]. Overexpression of miR3405p is known to reduce the expression of Nrf2 protein in the postexercise skeletal muscle of mice [83]. miR-424 and miR-542-3p are among the most abundant miRNAs in goat muscles [20]. Meanwhile, Garros et al. stated that miR-542 overexpression caused muscle wasting in mice, and reduced mitochondrial function [84]. The knockdown of miR-7 is shown to induce myoblast differentiation [29]. Also, miR-7 is downregulated in skeletal muscle of master athletes [85]. An increase in the expression of miR-7 levels has been recognized as a therapeutic target for opposing muscle dysfunction in DM1 [86]. miR-20a-5p and miR-20b-5p are involved in myoblast differentiation by downregulating the expression of E2F1 [87]. miR-17-5p, miR-378b, miR-199a-5p, and miR-7 may have key roles to play in muscle aging [88], and upregulation of miR-17-5p has been demonstrated to contributes to hypoxia-induced proliferation in human pulmonary artery smooth muscle cells through modulation of p21 and PTEN [52]. Huang et al. stated that miR-374a-5p was shown to negatively regulate MAPK6 expression, and miR-374a-5p may have protective effects against cardiac I/R injury in vivo, and H9C2 H/R injury in vitro [89]. Also, miR-129-5p blocks proliferation of vascular smooth muscle cells (VSMC) by decreasing the expression of Wnt5a [90] and cyclin-dependent kinase 6 (CDK6) cardiac myocytes [91]. Decreased in the expression of miR-27a/b resulted in a reduced myoblast proliferation through an increase in myostatin [92]. Overexpression of miR-296 acts as ceRNA for LncRNA CAREL and signi cantly further inhibits Trp53inp1 and ltm2a to induce proliferation of cardiomyocytes [93]. Overexpression of miR-423-3p in C2C12 myogenic differentiation lead to decreased expression of its target gene Cox6a2 as well as decreased in cellular ATP level [94]. miR-99b-3p is closely related to cancer cell proliferation [95,96]. LncRNA UCA1 promotes the progression of cardiac hypertrophy by competitively binding with miR-184 to enhance the expression of HOXA9 [97]. miR-296-3p enhances cardiac differentiation of embryonic stem cells [98]. Moreover, miR-128-3p that promotes differentiation of chicken SMSCs [99], was signi cantly downregulated after the knockdown of DICER1. Based on the reports from different works with regards to the functions of these miRNAs, we can deduce that downregulation of DICER1 increases proliferation but rather decreases differentiation.
In this study, we have made a signi cant contribution regarding the function of CDR1as and DICER1 in myogenesis. According to the sequencing results, there is a positive correlation between CDR1as and DICER1 differentially expressed miRNAs and mRNAs. Moreover, all the results stated indicates that the knockdown of CDR1as and DICER1 are related to the downregulation of myoblast differentiation. Also, we identi ed 11 putative miRNAs that can be sponged by CDR1as in SMSCs.

Conclusion
In summary, we report that CDR1as and DICER1 gene can induce myogenesis in SMSCs through differentiation. Moreover, the knockdown of CDR1as can also decrease the expression levels of DICER1 in SMSCs via vice versa. In future studies, we plan to investigate the functional regulatory pathways of some miRNAs and mRNAs associated with CDR1as and DICER1.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.