Binglin Yue1
Mei Liu
Mingxun Li
Hong Chen- 1College of Animal Science and Technology, Northwest A&F University, Yangling, China
- 2State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources, Guangxi University, Nanning, China
- 3College of Animal Science and Technology, Yangzhou University, Yangzhou, China
- 4Yunnan Academy of Grassland and Animal Science, Kunming, China
There is growing evidence that non-coding RNAs are emerging as critical regulators of skeletal muscle development. In order to reveal their functional roles and regulatory mechanisms, we constructed a lncRNA–miRNA–mRNA network according to the ceRNA (competitive endogenous RNA) theory, using our high-throughput sequencing data. Subsequently, the network analysis, GO (Gene Ontology) analysis, and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis were performed for functional annotation and exploration of lncRNA ceRNAs. The results uncovered a scale-free characteristics network which exhibited high functional specificity for bovine skeletal muscle development: co-expression lncRNAs were significantly enriched in muscle development related biological processes and the Wnt signaling pathway. Furthermore, GSEA (Gene Set Enrichment Analysis) indicated that the risk score has a tendency to associate with myogenesis, and differentially expressed RNAs were validated by qPCR, further confirming the credibility of our network. In summary, this study provides insights into lncRNA-mediated ceRNA function and mechanisms in bovine skeletal muscle development and will expand our understanding of lncRNA biology in mammals.
Introduction
Vertebrates’ skeletal muscles are mostly derived from paraxial mesodermal somites and undergo hyperplasia and hypertrophy process successively (Buckingham et al., 2003; Guo et al., 2015). Vertebrates’ skeletal muscle development is a complex process, which plays a crucial role in overall body metabolism, and a transcriptional hierarchy including the MRFs (myogenic regulatory factors) and MEF2 (myocyte enhancer factor 2) family precisely orchestrates it through coordinating the activities of a series of muscle genes (Pownall et al., 2002; Buckingham, 2006). Furthermore, several studies showed that post-transcriptional regulation also has a significant effect on skeletal muscle development (Cesana et al., 2011). However, the precise molecular mechanisms are still poorly understood (Te et al., 2010; He and Liu, 2013).
NcRNAs (non-coding RNAs) are RNA molecules that have little protein-coding capacity, but rather than being transcriptional noise, serves as master regulators in a variety of biological processes (Huarte, 2009; McHugh et al., 2015). Recently, accumulating evidence has suggested that ncRNAs like miRNA (microRNA) and lncRNA (long non-coding RNA) are emerging as critical regulators of skeletal muscle development (Williams et al., 2009; Matsumoto et al., 2017). Like mRNA-like transcripts longer than 200 nucleotides, lncRNAs have extremely abundant binding sites for microRNAs (Phelps et al., 2016), and LncRNAs such as linc-MD1 (Cesana et al., 2011), H19 (Kallen et al., 2013), lncMyoD (Gong et al., 2015) and lnc-mg (Zhu et al., 2017) showed relevant roles in myogenesis by acting as ceRNA (competing endogenous RNA). The ceRNA function, how noncoding RNAs and mRNAs regulate each other by competing for binding to shared miRNAs using partially complementary sequences known as MREs (miRNA response elements), which induces degradation of mRNA targets or translation repression at the post-transcriptional level (Salmena et al., 2011). Up to now, only few lncRNAs have been functionally well annotated (Butchart et al., 2016; Yu et al., 2017). Considering this, the aim of this study was to further identify a novel regulatory mechanism based on the ceRNA theory in bovine skeletal muscle development.
In recent years, RNA-Seq (RNA sequencing) based on high-throughput platforms has emerged as an efficient tool to screen genomic variation, associated with divergent phenotypic traits and identifies the molecular mechanisms underlying the domestication process (Dong et al., 2015). In the present study, our high-throughput sequencing data were processed to construct the lncRNA–miRNA–mRNA network, which has the potential to highlight specific molecular functions and mechanisms (Sun et al., 2013, 2016; Li et al., 2017; Figure 1). Afterward, functional enrichment analyses and annotation were performed to further explore the importance of lncRNAs in myogenesis. Considering this, our study expands the potential lncRNA ceRNA functions during bovine skeletal muscle development and enhances the understanding of non-coding RNA regulatory networks.
Figure 1. The main frame of data processing. (1) High-throughput sequencing data were pre-selected according to the | log2FC| ≥ 4(| log2FC| ≥ 2 for DEMs) and P-values < 0.001. (2) The DEMis-DEMs pairs were obtained using LncACTdb. (3) According to MFE theory, the DEMis-DELs pairs were screened using the RNAhybrid program. (4) the DEMis-DEMs pairs were merged with the DEMis-DELs pairs to obtain the co-expression RNAs. (5) 46 co-expression lncRNAs, 17 co-expression miRNAs, and 79 co-expression mRNAs were finally selected to construct the lncRNA–miRNA–mRNA network. DELs, differentially expressed lncRNAs; DEMis, differentially expressed miRNAs; DEMs, differentially expressed mRNAs; MFE, minimum free energy.
Results
Processed-Data Analysis
After data was processed, a total of 61 DELs (differentially expressed lncRNAs), 102 DEMis (differentially expressed miRNAs), and 2978 DEMs (differentially expressed mRNAs) were identified between two developmental stages in Chinese Qinchuan bovine longissimus muscles. Among these differentially expressed RNAs, DEMis were mostly down-regulated in samples from the embryonic stage to the adult stage, while DELs and DEMs showed the opposite results (Figures 2A–C). As described in the Materials and Methods section, we identified 46 co-expression lncRNAs, 17 co-expression miRNAs, and 79 co-expression mRNAs, which were selected to construct the lncRNA–miRNA–mRNA network (Figures 2D–F).
Figure 2. The differentially expressional pattern of RNAs from bovine embryonic and adult stages. (A-C) Volcano plot of 61 DELs (A), 102 DEMis (B), and 2978 DEMs (C) for embryo versus adult. Each point indicates an RNA. Green dots denote down-regulated RNAs, red points represent up-regulated RNAs under the same thresholds, and black points indicate RNAs that did not change significantly. The criteria is based on| log2FC| ≥ 4(| log2FC| ≥ 2 for DEMs) and P-values < 0.001. (D-F) Heatmap of 46 co-expression lncRNAs (D), 17 co-expression miRNAs (E), and 79 co-expression mRNAs (F) displaying when comparing embryonic and adult stages. Each column represents one sample, and each row refers to an RNA. The color legend is on the top-right of the figure, and compound relative abundances were standardized (Z score, shown in legend) prior to unsupervised hierarchical clustering of samples (rows). Blue indicates RNAs with a lower expression relative to the geometrical means; red indicates genes with a greater expression relative to the geometrical means.
LncRNA–miRNA–mRNA Network Construction and Visualization
The LncRNA–miRNA–mRNA network was constructed on the basis of these co-expression RNAs and visualized using Cytoscape. As shown in Figure 3A, 35 lncRNA nodes, 14 miRNA nodes, 69 mRNA nodes, and 141 edges composed the lncRNA–miRNA–mRNA network. In the study of networks, the degree of a node indicates the number of edges linked to the given node, and the degree distribution is the probability distribution of these degrees over the whole network. Following the network analysis, indicators of centrality such as DC (degree centrality), BC (betweenness centrality), and CC (closeness centrality) exist to determine the importance of a single node which possesses essential functions in a complex network. The degree of distribution of nodes in the network was investigated, and the power-law distribution with a slope of –1.141 and R-squared = 0.774 was observed (Figure 3B), suggesting that the network displayed scale-free characteristics typical of biological networks. Further comparison analysis then showed that there were significant differences in the DC, BC, and CC among lncRNAs, miRNAs, and mRNAs (P-value < 0.01 for DC, BC, and CC using Kruskal–Wallis test) (Figures 3C–E), indicating that lncRNAs and miRNAs tended to be critical in the context. Overall, the scale-free distribution and module characteristics of the entire network implied the presence of functionally important nodes in bovine skeletal muscle development.
Figure 3. The LncRNA–miRNA–mRNA network and its structural characteristics. (A) Global view of the network in bovine skeletal muscle development. This network consists of 141 edges among 35 lncRNA nods (rhombus), 14 miRNAs (circle), and 69 mRNAs (square). Red nodes represent upregulation, whereas green nodes represent downregulation. The node size typifies the level of edges enrichment: the higher the degree, the bigger the size is. (B) Degree distribution of the network. (C) The difference of degree centrality among lncRNAs, miRNAs, and mRNAs. The lncRNA nodes had a significantly higher degree centrality than mRNA nodes in the network. (D) The difference of betweenness centrality among lncRNAs, miRNAs, and mRNAs. The lncRNA nodes had a higher betweenness centrality than mRNA nodes in the network. (E) The difference of closeness centrality among lncRNAs, miRNAs, and mRNAs. The lncRNA nodes had a higher closeness centrality than mRNA nodes in the network. P-values were calculated based on Kruskal–Wallis test.
Molecular Function and Pathway Prediction
Results of the GO analysis revealed 54 BP (Biological Process) terms (Supplementary Table S1) which were mainly involved in three functional clusters, including nitric-oxide synthase biosynthetic process, gastrulation with mouth forming second, and axis specification. BP terms like muscle cell cellular homeostasis, mesoderm development, axis specification, and dorsal/ventral pattern formation were found in the enrichments, indicating the regulatory roles of these co-expression lncRNAs in bovine skeletal muscle development (Figure 4A). Following this, a pathway analysis demonstrated that these co-expression mRNAs were enriched in five pathways, such as longevity regulating pathway, adherens junction pathway, pancreatic cancer pathway, colorectal cancer pathway, and Wnt signaling pathway (Supplementary Table S2, Figure 4B), in which the Wnt signaling pathway is essential for embryonic muscle development and maintenance of adult skeletal muscle homeostasis (Figure 4C; Kanehisa et al., 2017). For instance, SFRP2 appears to function as a Wnt antagonist to prevent myoblasts from entering the terminal differentiation process in embryonic myogenesis (Ladher et al., 2000). As a downstream effector of the canonical Wnt signaling pathway, β-catenin (CTNNB1) affects muscle mass and slows fiber numbers in mice with conditional deletion of CTNNB1 in the muscle progenitor cells (Hutcheson et al., 2009). In addition, Wnt signaling activation was also suggested to induce satellite cell proliferation during skeletal muscle regeneration (Otto et al., 2008). These results suggested that the LncRNA–miRNA–mRNA network is related to skeletal muscle development.
Figure 4. Enrichment analysis of differentially expressed mRNAs. (A) Advanced bubble chart of GO-BP. Y-axis label represents terms name, and X-axis label typifies rich factor which is defined as the percentage of target genes per term. Size and color of the bubble are measured as count and –log10qvalue to represent the amount of differentially expressed mRNAs enriched in these terms and the enrichment significance, respectively, as well as three shapes of SampleGroup represent three functional clusters. (B) KEGG interaction network. Terms are represented as nodes based on the kappa score level (= 0.4), and node size indicates the term enrichment significant. Edges represent shared genes, the wider the edge, the larger the overlap is. The colors indicating different functional groups are shown at the right top. (C) Bos taurus Wnt signaling pathway map (Kanehisa et al., 2017). Upregulated genes in our network are marked in red.
GSEA Increased the Credibility of the LncRNA–miRNA–mRNA Network
Certain gene sets of the HALLMARK_MYOGENESIS from MSigDB, previously described to be involved in myogenesis, were used in Gene Set Enrichment Analysis (GSEA), and the adult high expression showed a tendency for a high enrichment score of the myogenesis-related gene set (NES: –1.83, FDR: 0.0003) (Supplementary Table S3 and Figure 5A). Moreover, the miR-377 targets were shown to be significantly located in the bottom of the pre-ranked gene-list in adult high expression (NES: –1.76, FDR: 0.006) by subjecting C3 MIR gene sets (Figure 5B), in which the targeted relationship of miR-377-SOD1 was predicted in the LncRNA–miRNA–mRNA network (Supplementary Table S4). These results increased the credibility of the LncRNA–miRNA–mRNA network to a certain degree.
Figure 5. Gene Set Enrichment Analysis (GSEA) of mRNA expression levels in bovine skeletal muscle. (A and B) GSEA result using HALLMARK_MYOGENESIS (A) and C3 MIR gene sets (B) Left panel: X-axis label represents “E”(Embryonic stage)/“A”(Adult stage) level, and Y-axis label represents the enrichment score of these mRNAs. Right panel: heat map of related mRNAs expression in two stages.
Validation of Differentially Expressed RNAs by qPCR
In the present study, a total of 10 lncRNAs, 7 miRNAs, and 7 mRNAs were randomly selected from the constructed network to perform validation by qRT-PCR. We confirmed stage-specific differences in the abundance of certain RNAs when comparing embryonic and adult longissimus muscle tissue samples using qPCR. qPCR expression patterns of the 10 lncRNAs, 7 miRNAs, and 7 mRNAs were basically similar to high-throughput sequencing results, which verified the reliability of sequencing results (Figures 6A–C). In order to explore the function of lncRNAs in myoblasts, 4 lncRNAs were selected to perform the qPCR for different tissues from different stages. We found that these lncRNAs had high expression levels in muscle tissues, and their expression was much higher at the adult stage compared to the embryonic stage (Figures 6D–G). These results implied their regulatory roles in muscle development.
Figure 6. Validation of putative RNAs by quantitative real-time PCR. (A) Absolute quantification for 10 lncRNAs in longissimus muscle tissue from the fetal and the adult Qinchuan cattle. (B) The expression of 7 miRNAs in longissimus muscle tissue of Qinchuan cattle at the embryonic stage compared with the adult stage. (C) The expression levels of 7 mRNAs in longissimus muscle tissue from the fetal and the adult Qinchuan cattle. (D-G) The expression of 4 randomly selected lncRNAs in different tissues from the fetal and the adult Qinchuan cattle. Values are means ± SEM for three individuals. P-values < 0.05.
Discussion
The development of bovine skeletal muscle is a complex process involving prenatal and postnatal patterns, which are mainly due to all muscle fiber formation, the muscle fiber size increases and new muscle fibers regenerate, respectively (Buckingham et al., 2003). In our research, two developmental stages were chosen: the embryonic 90-day stage, a prenatal generation when a large amount of muscle progenitor cell proliferation occurs; and the adult 24-month-old stage, during which myofibers are well established (Hocquette, 2010; Yue et al., 2017). Considering the differential expression of RNAs from these two stages, these RNAs might be associated with myogenesis, particularly lncRNAs and miRNAs that are known to have more developmental stage-specific expression patterns than protein-coding genes (Guttman and Rinn, 2012). According to the study, a series of exquisitely regulated and orchestrated changes happened at multiple levels including mRNA transcription and translation, protein stability, and degradation from the embryonic to the adult stage (Bismuth and Relaix, 2010). CeRNA, a well-known regulatory mechanism sets up an extensive regulatory network among the transcriptome, in which lncRNAs have been proposed to sponge miRNAs and thereby regulate other transcripts containing MREs (Denzler et al., 2016). In spite of their low abundance and/or nuclear localization, thousands of lncRNAs possess developmental stage-specific expression patterns and localizations (Mercer and Mattick, 2013; Koufariotis et al., 2015). Based on this reasoning, a global triple network was visualized in our study to predict potential lncRNA-mediated networks in myogenesis (Shannon et al., 2003).
In this study, natural ceRNA crosstalk interactions, with prominent expression shifts under the two stages, were screened bioinformatically and topological properties of the network verified more accepted ceRNA activity of lncRNAs (Gursoy et al., 2008). These observations suggested that the constructed ceRNA network could serve as a powerful prediction tool to explore the lncRNAs’ ceRNA function and role in bovine skeletal muscle development. Accumulated research showed that lncRNA ceRNAs have been implicated in myogenesis: lincMD1 controls muscle differentiation in human and mouse myoblasts directly by competing for miR-133 and miR-135, regulating the expression of MAML1 and MEF2C mRNAs, respectively (Cesana et al., 2011); acting as a molecular sponge, H19 lncRNA, which is highly expressed in the developing embryo and adult muscle, was found to bind to let-7 and inhibits its myogenesis function (Kallen et al., 2013); Zhu et al. (2017) reported that as a ceRNA for microRNA-125b, lnc-mg facilitates myogenesis with IGF2; while current reports showed that lnc133b/miR-133b/IGF1R axis is a potential pathway for promoting satellite cell proliferation and repressing their differentiation through the ceRNA mechanism (Jin et al., 2017); and Linc-smad7 was reported to promote myoblast differentiation and muscle regeneration via sponging miR-125b (Song et al., 2018). These findings indicate that individual lncRNAs may be potent natural miRNA sponges in certain settings, and lncRNA-associated ceRNA crosstalk will likely shift under different specific conditions.
To further research, we used an efficient way to infer the potential function of lncRNAs by studying miRNAs and/or mRNAs annotated functions (Wang et al., 2016; Zhou et al., 2016). In such a scenario, the functional enrichment analysis of 69 co-expression mRNAs in the network were performed including GO and KEGG pathway analysis. It was found that significant GO categories and pathways of DEMs were mainly enriched in myogenesis, suggesting the possible roles of these lncRNAs in bovine skeletal muscle development. Intriguingly, associated genes in myogenesis-related GO terms are highly consistent with those in a Wnt signaling pathway, which plays an essential role during embryonic muscle development and in the maintenance of adult skeletal muscle homeostasis (von Maltzahn et al., 2012). For instance, APC dampens canonical Wnt signaling to regulate mouse muscle stem cell proliferation and quiescence (Parisi et al., 2015); CTNNB1, the coding gene of β-catenin which is a key component of the canonical Wnt signaling pathway, was reported to alter the expression of some myogenic markers and downregulates myogenesis by regulating the expression of PPARG (Gustafsson et al., 2002; Okamura et al., 2009); as a Wnt antagonist, SFRP2 is upregulated during injury regeneration and exerts a differentiation inhibitory action on fibroblasts (Zhao and Hoffman, 2004; Descamps et al., 2008). However, previous in vitro studies have revealed that miR-199a-5p was affects myoblast differentiation by targeting several myogenic cell proliferation and differentiation regulatory factors of the Wnt signaling pathway, such as FZD4, JAG1, and WNT2 (Alexander et al., 2013). Functional analysis has suggested the possible correlation between the lncRNA-associated ceRNA crosstalk and the effects of Wnt on skeletal muscles.
The limitations of the network should be acknowledged. Early screening of differential expressed mRNAs was so strict that many myogenesis-related mRNAs were missed. Thus, Gene Set Enrichment Analysis was carried out to identify the associated gene effect from the embryonic to the adult stage, and the results determined a significant correlation between the adult high expression and the myogenesis-related gene set. Interestingly, our network also showed the mRNA-upregulated expression. Previous experiments suggested that SMAD4, whose protein levels increased dramatically in myoblasts of aged mice, is directly down-regulated by miR-431 and miR-26a to promote muscle differentiation and regeneration (Dey et al., 2012; Lee et al., 2015); SRF expression is up-regulated after birth and subsequently stabilizes in adulthood, and SRF inhibits muscle cell proliferation and differentiation in vitro and in vivo (Wang et al., 2002; Li et al., 2005). Accordingly, we infer that some mRNAs upregulated at the adult stage may exert a differentiation inhibitory action on fibroblasts, which are regulated by signaling pathways such as the Wnt signaling (Descamps et al., 2008).
To our knowledge, this is the first study that constructs a lncRNA-associated ceRNA network in Chinese Qinchuan bovine skeletal muscle development. Many miRNA-target gene pairs (miR-335-RB1, miR-122-SRF, miR-377-SOD1) were literately evaluated, and the validation of differentially expressed RNAs by qPCR further increases the credibility of our network (Scarola et al., 2010; Wang W. et al., 2015; Zeng et al., 2015). However, it has been suggested that, the ceRNA effectiveness was influenced by various factors including miRNA- and ceRNA- abundance, RBPs (RNA binding proteins) and subcellular localization, as well as the number and binding affinities of the MREs (Mukherji et al., 2011; Wee et al., 2012; Liu and Miao, 2016), so further experimental studies should be conducted to uncover the lncRNAs ceRNA function in bovine skeletal muscle development as a next step.
Materials and Methods
Ethics Statement
The animal care and experiments were conducted according to the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, 2004) China and approved by the Institutional Animal Care and Use Committee (College of Animal Science and Technology, Northwest A&F University, China). All samples (heart, liver, spleen, lung, kidney, and longissimus muscle) from Qinchuan cattle at the embryonic stage (90 days) and the adult stage (24-month-old) were collected at Shannxi Kingbull Animal Husbandry Co., Ltd. (Baoji, China). Adult animals were humanely killed where necessary, to ameliorate suffering and were not fed the night before they were slaughtered.
High-Throughput Sequencing Data
The sequence-based data were obtained from our high-throughput sequencing results, including lncRNA (Supplementary Table S5), miRNA (Supplementary Table S6), and mRNA (Supplementary Table S7) libraries from Chinese Qinchuan bovine longissimus muscles in two developmental stages (embryonic 90-day and adult 24-month-old). In the process, six RNA samples, that passed the quality control, were pooled to obtain an “averaged” transcriptome at each developmental stage (mRNAs were not “averaged”), and the databases were based on the platform of Solexa technology (Beijing Genomics Institute, China) and Illumina HiSeq 2500 Technology (LC Sciences, Houston, TX, United States).
Data Processing
The DEMis, DELs, and DEMs between two developmental stages were selected according to the | log2FC|≥ 4(| log2FC|≥ 2 for DEMs) and P-values < 0.001, which were used to control the number of differentially expressed RNAs within a manageable but useful range, not too many or too few. Subsequently, the DEMis were handled using the LncACTdb (LncRNA-associated Competing Triplets DataBase) to provide competition based lncRNA–miRNA–mRNA interactions, from which we obtained the DEMis-DEMs pairs (Wang P. et al., 2015). Then the DEMis-DELs pairs were predicted using the RNAhybrid program based on the MFE (minimum free energy) of miRNA–lncRNA duplexes (Rehmsmeier et al., 2004). Data for these interactions were downloaded from miRBase (Kozomara and Griffiths-Jones, 2014) and NONCODE (Zhao et al., 2016). In order to acquire high-quality target lncRNAs, DELs that had perfect nucleotide pairing between the 2nd and 8th positions of the 5’end of DEMis sequences were selected. Finally, the DEMis-DEMs pairs were merged with the DEMis-DELs pairs to obtain the co-expression RNAs which were used to construct and visualize the network.
Construction and Analysis of the lncRNA–miRNA–mRNA Network
The lncRNA–miRNA–mRNA network was constructed and visualized using Cytoscape software based on the ceRNA theory (Shannon et al., 2003). Here, nodes and edges were used to represent large biological data in an intuitive way, in which each node represents a biological molecule, the edges stand for the interactions between nodes, and the node degrees indicating the number of edges linked to a given node were calculated to exploit the hub nodes that possess essential biological functions (Kohl et al., 2011). To explore the structure and feature of the lncRNA–miRNA–mRNA competing triplets, a network analysis was performed using a Cytoscape plug-in called “Network Analyzer” (Assenov et al., 2008). Topological parameters like DC, BC, and CC, which are standard measures of centrality in a network, were assessed here (The DC is defined as the number of links incident upon a node, the BC for each node is calculated as the number of these shortest paths that pass through the node, and the CC of a node is the total of the length of the shortest paths between the node and all other nodes in the network).
Functional Enrichment Analysis
In order to assess functional enrichment in the co-expression competing triplet, Cytoscape plug-in ClueGO was used to constitute the interaction network based on the GO and KEGG database (Mlecnik et al., 2017). In this section, functional annotation of the 69 co-expression mRNAs were performed1 to infer the function of each co-expression lncRNA, and the statistical test used for the enrichment was based on the two-sided hypergeometric test with a Bonferroni step-down correction and kappa score of 0.4.
GSEA
Gene Set Enrichment Analysis was performed based on our mRNA-seq data using the GSEA software version 3.02. The mRNA-seq data were pre-ranked to form a gene-list according to their expression levels between two developmental stages. Following this, certain gene sets from MSigDB (H and C3 MIR gene set) were mapped to the pre-ranked gene-list to calculate the ES (enrichment score), in which 1000 permutations were used to calculate significance, and a certain gene set with FDR < 0.01 was considered significant.
Total RNA Extraction and Quantitative Real-Time PCR Analysis of RNAs
Using Trizol reagent (TakaRa; Japan), total RNA was extracted from Chinese Qinchuan bovine tissues. A total of 1000 ng total RNA was reverse transcribed into cDNA with the PrimeScript RT Kit (Takara), and miRNA-specific stem-loop primers (Supplementary Table S8) were used for reverse-transcribed cDNA of miRNAs. For RNA quantification, all measurements were replicated three times in a Bio-Rad master cycler using SYBR Green II Master Mix Reagent Kit (Takara) following the manufacturer’s instructions. qPCR primers are shown in Supplementary Table S8, and U6 (for miRNA) and GAPDH were used as the internal control for normalization of the data. The expression levels of lncRNAs, miRNAs and mRNAs were calculated using the 2−ΔΔCt method.
Data Availability
The datasets generated for this study can be found in NCBI, SRR1818309, SRR1818415, and SRR1818416.
Author Contributions
HC conceived the study idea and designed the work. BY, HL, MeL, JW, and MiL analyzed and interpreted the data. BY wrote the manuscript. BH provided partial samples. CL revised the manuscript critically. All authors read and approved the final manuscript.
Funding
This study was supported by the National Natural Science Foundation of China (No. 31772574), the Program of National Beef Cattle and Yak Industrial Technology System (CARS-37) and the Program of Yunling Scholar.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00091/full#supplementary-material
Footnotes
References
Alexander, M. S., Kawahara, G., Motohashi, N., Casar, J. C., Eisenberg, I., Myers, J. A., et al. (2013). MicroRNA-199a is induced in dystrophic muscle and affects WNT signaling, cell proliferation, and myogenic differentiation. Cell Death Differ. 20, 1194–1208. doi: 10.1038/cdd.2013.62
Assenov, Y., Ramirez, F., Schelhorn, S. E., Lengauer, T., and Albrecht, M. (2008). Computing topological parameters of biological networks. Bioinformatics 24, 282–284. doi: 10.1093/bioinformatics/btm554
Bismuth, K., and Relaix, F. (2010). Genetic regulation of skeletal muscle development. Exp. Cell Res. 316, 3081–3086. doi: 10.1016/j.yexcr.2010.08.018
Buckingham, M. (2006). Myogenic progenitor cells and skeletal myogenesis in vertebrates. Curr. Opin. Genet. Dev. 16, 525–532. doi: 10.1016/j.gde.2006.08.008
Buckingham, M., Bajard, L., Chang, T., Daubas, P., Hadchouel, J., Meilhac, S., et al. (2003). The formation of skeletal muscle: from somite to limb. J. Anat. 202, 59–68. doi: 10.1046/j.1469-7580.2003.00139.x
Butchart, L. C., Fox, A., Shavlakadze, T., and Grounds, M. D. (2016). The long and short of non-coding RNAs during post-natal growth and differentiation of skeletal muscles: Focus on lncRNA and miRNAs. Differentiation 92, 237–248. doi: 10.1016/j.diff.2016.05.003
Cesana, M., Cacchiarelli, D., Legnini, I., Santini, T., Sthandier, O., Chinappi, M., et al. (2011). A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 147, 358–369. doi: 10.1016/j.cell.2011.09.028
Denzler, R., McGeary, S. E., Title, A. C., Agarwal, V., Bartel, D. P., and Stoffel, M. (2016). Impact of MicroRNA levels, target-site complementarity, and cooperativity on competing endogenous RNA-regulated gene expression. Mol. Cell 64, 565–579. doi: 10.1016/j.molcel.2016.09.027
Descamps, S., Arzouk, H., Bacou, F., Bernardi, H., Fedon, Y., Gay, S., et al. (2008). Inhibition of myoblast differentiation by Sfrp1 and Sfrp2. Cell Tissue Res. 332, 299–306. doi: 10.1007/s00441-008-0574-z
Dey, B. K., Gagan, J., Yan, Z., and Dutta, A. (2012). miR-26a is required for skeletal muscle differentiation and regeneration in mice. Genes Dev. 26, 2180–2191. doi: 10.1101/gad.198085.112
Dong, Y., Zhang, X., Xie, M., Arefnezhad, B., Wang, Z., Wang, W., et al. (2015). Reference genome of wild goat (capra aegagrus) and sequencing of goat breeds provide insight into genic basis of goat domestication. BMC Genomics 16:431. doi: 10.1186/s12864-015-1606-1
Gong, C., Li, Z., Ramanujan, K., Clay, I., Zhang, Y., Lemire-Brachat, S., et al. (2015). A long non-coding RNA, LncMyoD, regulates skeletal muscle differentiation by blocking IMP2-mediated mRNA translation. Dev. Cell 34, 181–191. doi: 10.1016/j.devcel.2015.05.009
Guo, B., Greenwood, P. L., Cafe, L. M., Zhou, G., Zhang, W., and Dalrymple, B. P. (2015). Transcriptome analysis of cattle muscle identifies potential markers for skeletal muscle growth rate and major cell types. BMC Genomics 16:177. doi: 10.1186/s12864-015-1403-x
Gursoy, A., Keskin, O., and Nussinov, R. (2008). Topological properties of protein interaction networks from a structural perspective. Biochem. Soc. Trans. 36, 1398–1403. doi: 10.1042/BST0361398
Gustafsson, M. K., Pan, H., Pinney, D. F., Liu, Y., Lewandowski, A., Epstein, D. J., et al. (2002). Myf5 is a direct target of long-range Shh signaling and Gli regulation for muscle specification. Genes Dev. 16, 114–126. doi: 10.1101/gad.940702
Guttman, M., and Rinn, J. L. (2012). Modular regulatory principles of large non-coding RNAs. Nature 482, 339–346. doi: 10.1038/nature10887
He, H., and Liu, X. (2013). Characterization of transcriptional complexity during longissimus muscle development in bovines using high-throughput sequencing. PLoS One 8:e64356. doi: 10.1371/journal.pone.0064356
Hocquette, J. F. (2010). Endocrine and metabolic regulation of muscle growth and body composition in cattle. Animal 4, 1797–1809. doi: 10.1017/S1751731110001448
Huarte, M. (2009). Journal club. A biologist looks at new functions for non-coding RNAs. Nature 459:487. doi: 10.1038/459487e
Hutcheson, D. A., Zhao, J., Merrell, A., Haldar, M., and Kardon, G. (2009). Embryonic and fetal limb myogenic cells are derived from developmentally distinct progenitors and have different requirements for beta-catenin. Genes Dev. 23, 997–1013. doi: 10.1101/gad.1769009
Jin, C. F., Li, Y., Ding, X. B., Li, X., Zhang, L. L., Liu, X. F., et al. (2017). lnc133b, a novel, long non-coding RNA, regulates bovine skeletal muscle satellite cell proliferation and differentiation by mediating miR-133b. Gene 630, 35–43. doi: 10.1016/j.gene.2017.07.066
Kallen, A. N., Zhou, X. B., Xu, J., Qiao, C., Ma, J., Yan, L., et al. (2013). The imprinted H19 lncRNA antagonizes let-7 microRNAs. Mol. Cell 52, 101–112. doi: 10.1016/j.molcel.2013.08.027
Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361. doi: 10.1093/nar/gkw1092
Kohl, M., Wiese, S., and Warscheid, B. (2011). Cytoscape: software for visualization and analysis of biological networks. Methods Mol. Biol. 696, 291–303. doi: 10.1007/978-1-60761-987-1_18
Koufariotis, L. T., Chen, Y. P., Chamberlain, A., Vander, J. C., and Hayes, B. J. (2015). A catalogue of novel bovine long noncoding RNA across 18 tissues. PLoS One 10:e0141225. doi: 10.1371/journal.pone.0141225
Kozomara, A., and Griffiths-Jones, S. (2014). miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi: 10.1093/nar/gkt1181
Ladher, R. K., Church, V. L., Allen, S., Robson, L., Abdelfattah, A., Brown, N. A., et al. (2000). Cloning and expression of the Wnt antagonists Sfrp-2 and Frzb during chick development. Dev. Biol. 218, 183–198. doi: 10.1006/dbio.1999.9586
Lee, K. P., Shin, Y. J., Panda, A. C., Abdelmohsen, K., Kim, J. Y., Lee, S. M., et al. (2015). miR-431 promotes differentiation and regeneration of old skeletal muscle by targeting Smad4. Genes Dev. 29, 1605–1617. doi: 10.1101/gad.263574.115
Li, H., Wei, X., Yang, J., Dong, D., Huang, Y., Lan, X., et al. (2017). Developmental transcriptome profiling of bovine muscle tissue reveals an abundant GosB that regulates myoblast proliferation and apoptosis. Oncotarget 8, 32083–32100. doi: 10.18632/oncotarget.16644
Li, S., Czubryt, M. P., McAnally, J., Bassel-Duby, R., Richardson, J. A., Wiebel, F. F., et al. (2005). Requirement for serum response factor for skeletal muscle growth and maturation revealed by tissue-specific gene deletion in mice. Proc. Natl. Acad. Sci. U.S.A. 102, 1082–1087. doi: 10.1073/pnas.0409103102
Liu, Z., and Miao, H. (2016). Prediction of protein-RNA interactions using sequence and structure descriptors. Neurocomputing 206, 28–34. doi: 10.1016/j.neucom.2015.11.105
Matsumoto, A., Pasut, A., Matsumoto, M., Yamashita, R., Fung, J., Monteleone, E., et al. (2017). mTORC1 and muscle regeneration are regulated by the LINC00961-encoded SPAR polypeptide. Nature 541, 228–232. doi: 10.1038/nature21034
McHugh, C. A., Chen, C. K., Chow, A., Surka, C. F., Tran, C., McDonel, P., et al. (2015). The Xist lncRNA interacts directly with SHARP to silence transcription through HDAC3. Nature 521, 232–236. doi: 10.1038/nature14443
Mercer, T. R., and Mattick, J. S. (2013). Structure and function of long noncoding RNAs in epigenetic regulation. Nat. Struct. Mol. Biol. 20, 300–307. doi: 10.1038/nsmb.2480
Mlecnik, B., Galon, J., and Bindea, G. (2017). Comprehensive functional analysis of large lists of genes and proteins. J. Proteomics 171, 2–10. doi: 10.1016/j.jprot.2017.03.016
Mukherji, S., Ebert, M. S., Zheng, G. X., Tsang, J. S., Sharp, P. A., and van Oudenaarden, A. (2011). MicroRNAs can generate thresholds in target gene expression. Nat. Genet. 43, 854–859. doi: 10.1038/ng.905
Okamura, M., Kudo, H., Wakabayashi, K., Tanaka, T., Nonaka, A., Uchida, A., et al. (2009). COUP-TFII acts downstream of Wnt/beta-catenin signal to silence PPARgamma gene expression and repress adipogenesis. Proc. Natl. Acad. Sci. U.S.A. 106, 5819–5824. doi: 10.1073/pnas.0901676106
Otto, A., Schmidt, C., Luke, G., Allen, S., Valasek, P., Muntoni, F., et al. (2008). Canonical Wnt signalling induces satellite-cell proliferation during adult skeletal muscle regeneration. J. Cell Sci. 121, 2939–2950. doi: 10.1242/jcs.026534
Parisi, A., Lacour, F., Giordani, L., Colnot, S., Maire, P., and Le Grand, F. (2015). APC is required for muscle stem cell proliferation and skeletal muscle tissue repair. J. Cell Biol. 210, 717–726. doi: 10.1083/jcb.201501053
Phelps, M., Coss, C., Wang, H., and Cook, M. (2016). Registered report: Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. eLife 5:e12470. doi: 10.7554/eLife.12470
Pownall, M. E., Gustafsson, M. K., and Emerson, C. J. (2002). Myogenic regulatory factors and the specification of muscle progenitors in vertebrate embryos. Annu. Rev. Cell Dev. Biol. 18, 747–783. doi: 10.1146/annurev.cellbio.18.012502.105758
Rehmsmeier, M., Steffen, P., Hochsmann, M., and Giegerich, R. (2004). Fast and effective prediction of microRNA/target duplexes. RNA 10, 1507–1517. doi: 10.1261/rna.5248604
Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014
Scarola, M., Schoeftner, S., Schneider, C., and Benetti, R. (2010). miR-335 directly targets Rb1 (pRb/p105) in a proximal connection to p53-dependent stress response. Cancer Res. 70, 6925–6933. doi: 10.1158/0008-5472.CAN-10-0141
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Song, C., Wang, J., Ma, Y., Yang, Z., Dong, D., Li, H., et al. (2018). Linc-smad7 promotes myoblast differentiation and muscle regeneration via sponging miR-125b. Epigenetics 13, 591–604. doi: 10.1080/15592294.2018.1481705
Sun, J., Li, M., Li, Z., Xue, J., Lan, X., Zhang, C., et al. (2013). Identification and profiling of conserved and novel microRNAs from Chinese Qinchuan bovine longissimus thoracis. BMC Genomics 14:42. doi: 10.1186/1471-2164-14-42
Sun, X., Li, M., Sun, Y., Cai, H., Lan, X., Huang, Y., et al. (2016). The developmental transcriptome sequencing of bovine skeletal muscle reveals a long noncoding RNA, lncMD, promotes muscle differentiation by sponging miR-125b. Biochim. Biophys. Acta 1863, 2835–2845. doi: 10.1016/j.bbamcr.2016.08.014
Te, P. M., Keuning, E., Hulsegge, B., Hoving-Bolink, A. H., Evans, G., and Mulder, H. A. (2010). Longissimus muscle transcriptome profiles related to carcass and meat quality traits in fresh meat Pietrain carcasses. J. Anim. Sci. 88, 4044–4055. doi: 10.2527/jas.2010-2952
von Maltzahn, J., Chang, N. C., Bentzinger, C. F., and Rudnicki, M. A. (2012). Wnt signaling in myogenesis. Trends Cell Biol. 22, 602–609. doi: 10.1016/j.tcb.2012.07.008
Wang, D., Passier, R., Liu, Z. P., Shin, C. H., Wang, Z., Li, S., et al. (2002). Regulation of cardiac growth and development by SRF and its cofactors. Cold Spring Harb. Symp. Quant. Biol. 67, 97–105. doi: 10.1101/sqb.2002.67.97
Wang, H., Niu, L., Jiang, S., Zhai, J., Wang, P., Kong, F., et al. (2016). Comprehensive analysis of aberrantly expressed profiles of lncRNAs and miRNAs with associated ceRNA network in muscle-invasive bladder cancer. Oncotarget 7, 86174–86185. doi: 10.18632/oncotarget.13363
Wang, P., Ning, S., Zhang, Y., Li, R., Ye, J., Zhao, Z., et al. (2015). Identification of lncRNA-associated competing triplets reveals global patterns and prognostic markers for cancer. Nucleic Acids Res. 43, 3478–3489. doi: 10.1093/nar/gkv233
Wang, W., Ding, X. Q., Gu, T. T., Song, L., Li, J. M., Xue, Q. C., et al. (2015). Pterostilbene and allopurinol reduce fructose-induced podocyte oxidative stress and inflammation via microRNA-377. Free Radic. Biol. Med. 83, 214–226. doi: 10.1016/j.freeradbiomed.2015.02.029
Wee, L. M., Flores-Jasso, C. F., Salomon, W. E., and Zamore, P. D. (2012). Argonaute divides its RNA guide into domains with distinct functions and RNA-binding properties. Cell 151, 1055–1067. doi: 10.1016/j.cell.2012.10.036
Williams, A. H., Valdez, G., Moresi, V., Qi, X., McAnally, J., Elliott, J. L., et al. (2009). MicroRNA-206 delays ALS progression and promotes regeneration of neuromuscular synapses in mice. Science 326, 1549–1554. doi: 10.1126/science.1181046
Yu, X., Zhang, Y., Li, T., Ma, Z., Jia, H., Chen, Q., et al. (2017). Long non-coding RNA Linc-RAM enhances myogenic differentiation by interacting with MyoD. Nat. Commun. 8:14016. doi: 10.1038/ncomms14016
Yue, B., Wu, J., Wang, Y., Zhang, C., Fang, X., and Chen, H. (2017). Expression profiles analysis and functional characterization of microRNA-660 in skeletal muscle differentiation. J. Cell. Biochem. 118, 2387–2394. doi: 10.1002/jcb.25901
Zeng, C., Wang, Y. L., Xie, C., Sang, Y., Li, T. J., Zhang, M., et al. (2015). Identification of a novel TGF-beta-miR-122-fibronectin 1/serum response factor signaling cascade and its implication in hepatic fibrogenesis. Oncotarget 6, 12224–12233. doi: 10.18632/oncotarget.3652
Zhao, P., and Hoffman, E. P. (2004). Embryonic myogenesis pathways in muscle regeneration. Dev. Dyn. 229, 380–392. doi: 10.1002/dvdy.10457
Zhao, Y., Li, H., Fang, S., Kang, Y., Wu, W., Hao, Y., et al. (2016). NONCODE 2016: an informative and valuable data source of long non-coding RNAs. Nucleic Acids Res. 44, D203–D208. doi: 10.1093/nar/gkv1252
Zhou, M., Wang, X., Shi, H., Cheng, L., Wang, Z., Zhao, H., et al. (2016). Characterization of long non-coding RNA-associated ceRNA network to reveal potential prognostic lncRNA biomarkers in human ovarian cancer. Oncotarget 7, 12598–12611. doi: 10.18632/oncotarget.7181
Keywords: bovine, ceRNA, lncRNA, skeletal muscle, enrichment analysis
Citation: Yue B, Li H, Liu M, Wu J, Li M, Lei C, Huang B and Chen H (2019) Characterization of lncRNA–miRNA–mRNA Network to Reveal Potential Functional ceRNAs in Bovine Skeletal Muscle. Front. Genet. 10:91. doi: 10.3389/fgene.2019.00091
Received: 15 October 2018; Accepted: 29 January 2019;
Published: 20 February 2019.
Edited by:
Zhifu Sun, Mayo Clinic, United StatesCopyright © 2019 Yue, Li, Liu, Wu, Li, Lei, Huang and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hong Chen, chenhong1212@263.net; easyman123@sina.com