Growth curves and histological analyses of Qingyu pigs
To better understand the growth and development of Qingyu pigs, the body weight of 126 Qingyu pigs was fitted with three nonlinear growth models, i.e., Logistic, Von Bertalanffy and Gompertz curve models (Fig. 1a, see Additional file 1). All three models showed a good fit with a typical sigmoidal curve, although the Von Bertalanffy curve showed the highest R2 value with the best goodness of fit (R2 = 0.9971) (see Additional file 2). The inflection point analysis of the growth curve indicated that Qingyu pigs reached the MGI stage at 156.40 days of age, and the average body weight of these pigs at this time point was 51.73 kg. Similarly, the GRI and RSI stages were reached at 23.82 days of age with 3.14 kg body weight and 288.97 days of age with 107.03 kg body weight, respectively (see Additional file 1). Additionally, the maximum growth rate of Qingyu pigs was 465.61 g per day (Fig. 1b, see Additional file 3). During muscle development, the mean cross-section area (CSA) of the longissimus dorsi muscle increased from 270 μm2 at GRI to 880 μm2 and 1500 μm2 at MGI and RSI, respectively (Fig. 1c).
Expression pattern of mRNAs and lncRNAs
To comprehensively identify transcripts related to the physiological differences in Qingyu pigs among the MGI, GRI and RSI stages, a total of 9 libraries were constructed (three libraries at each stage). A total of 133.79 Gb data were generated (Table 1), with an average of 99.11 million raw reads per sample sequenced at approximately 4× coverage. To explore the differences between lncRNAs and mRNAs, the average lncRNA and mRNA levels were transformed to log2 (CPM + 1). The results showed that the average level of lncRNAs was lower than that of mRNAs (see Additional file 4), consistent with the expression pattern obtained in other tissues [21].
Table 1 Summary of RNA_seq in Qingyu pigs
Samples
|
Raw_Yield(G)
|
Raw Reads (M)
|
Clean_Yield(G)
|
Clean Reads (M)
|
Clean_Q20(%)
|
Clean_GC(%)
|
GRI-1
|
13.708
|
91.38
|
13.347
|
89.96
|
98.04
|
50.19
|
GRI-2
|
14.745
|
98.3
|
14.271
|
95.79
|
97.79
|
49.62
|
GRI-3
|
14.26
|
95.06
|
13.814
|
92.58
|
97.85
|
50.46
|
MGI-1
|
17.549
|
117
|
17.103
|
114.87
|
98.01
|
50.67
|
MGI-2
|
16.739
|
111.59
|
16.341
|
109.99
|
98.06
|
50.46
|
MGI-3
|
14.605
|
97.37
|
14.237
|
95.48
|
97.84
|
50.16
|
RSI-1
|
14.72
|
98.13
|
14.28
|
96.15
|
97.97
|
50.13
|
RSI-2
|
12.683
|
84.56
|
12.332
|
82.72
|
97.9
|
50.63
|
RSI-3
|
14.789
|
98.6
|
14.4
|
96.74
|
97.9
|
51.29
|
GRI, the inflection of gradual increase stage to rapid increasing stage; MGI, the inflection point with the maximum growth rate; RSI, the inflection point of rapid increasing stage to slowly increasing stage.
The mRNAs and lncRNAs are differentially expressed during skeletal muscle development; however, little research has been conducted on skeletal muscle transcriptome based on timepoints according to growth curve for porcine. Therefore, we sought to explore the expression profiles of mRNAs and lncRNAs during muscle development at the GRI, MGI and RSI stages. A total of 14,530 mRNAs and 11,970 lncRNAs were expressed at the three stages (see Additional file 4). Among these, 14,475 mRNAs and 11,955 lncRNAs were detected at the GRI stage, 14,446 mRNAs and 11,949 lncRNAs at the MGI stage and 14,439 mRNAs and 11,942 lncRNAs at the RSI stage.
Because all coding and noncoding transcripts were quantified in parallel, our expression profile also allowed the assessment and comparison of temporal changes in lncRNAs and mRNAs during muscle development. Firstly, we performed hierarchical clustering analyses on transcripts showing maximal expression in three different developmental stages. The mRNA expression profiles readily separated all samples into two distinct groups, as expected, and samples clustered tightly within each stage repetition (Fig. 2). MGI and RSI were clustered together in one branch distinct from the GRI. Interestingly, a nearly identical pattern of sample clustering was observed for regulated lncRNAs (Fig. 2a), indicating that expression profiles of lncRNAs could serve as a developmental signature, similar to protein-coding mRNAs. Consistently, principal component analysis (PCA) of all regulated transcripts, including mRNAs or lncRNAs, readily separated all samples into three distinct groups (Fig. 2b). These patterns suggest that regulated lncRNA and mRNA transcriptomes function coordinately in related physiological processes, and our samples were highly reliable for subsequent analysis.
Functional enrichment analysis of differentially expressed mRNAs (DEGs)
The results of growth curve analysis indicated that Qingyu pigs reached the maximum growth rate at the MGI stage. We then investigated the DEGs and conducted functional enrichment analysis of these DEGs between the GRI vs. MGI (GRI–MGI) group and RSI vs. MGI (RSI–MGI) group to identify the physiological changes before and after reaching the MGI stage. A total of 645 and 323 DEGs were identified in the GRI–MGI and RSI–MGI groups, respectively. Among these DEGs, 318 were up-regulated and 327 were down-regulated in the GRI–MGI group (Fig. 3), whereas 177 were up-regulated and 146 were down-regulated in the RSI–MGI group (Fig. 3). Consistently, a distinct expression pattern was found between GRI vs. MGI because more DEGs were detected in the GRI–MGI group than in the RSI–MGI group. These results also confirmed the results of hierarchical clustering analysis and PCA, indicating that a massive physiological change occurred at the early muscle development stage.
Next, we separately performed Gene Ontology (GO) and KEGG pathway enrichment of DEGs in the GRI–MGI and RSI–MGI groups. As expected, DEGs up-regulated in the GRI–MGI group were enriched in skeletal system development (GO:0001501), myosin light chain binding (GO:0032027) and hallmark myogenesis (M5909). Additionally, DEGs up-regulated in the GRI–MGI group (i.e., genes showing higher expression at the GRI stage) were mainly enriched in immune related terms, such as cell activation involved in immune response (GO:0002263), activation of immune response (GO:0002253), immune response-activating signal transduction (GO:0002757), immune response-regulating signaling pathway (GO:0002764) and hallmark complement (M5921) (Fig. 3a, Additional file 5). Similarly, DEGs down-regulated in the GRI–MGI group (i.e., genes showing higher expression at the MGI stage) were enriched in muscle system process (GO:0003012), muscle contraction (GO:0006936), striated muscle contraction (GO:0006941) and myofibril (GO:0030016). Additionally, these down-regulated DEGs were also enriched in amino acid metabolism and glycogen metabolism, e.g., cellular amino acid catabolic process (GO:0009063), cellular amino acid metabolic process (GO:0006520), metabolism of amino acids and as well as derivatives (R-HSA-71291), glycogen metabolic process (GO:0005977) and glycogen metabolism (R-HSA-8982491) (Fig. 3a, Additional file 5). On the other hand, DEGs up-regulated in the RSI–MGI group (i.e., genes showing higher expression at the RSI stage) were mostly enriched in the regulation of lipid metabolic process (GO:0019216), regulation of lipid biosynthetic process (GO:0046890), fat cell differentiation (GO:0045444), metabolism of lipids (R-HSA-556833) and hallmark adipogenesis (M5905) (Fig. 3b, Additional file 6), whereas DEGs down-regulated in the RSI–MGI group (i.e., genes showing higher expression at the MGI stage) were enriched in the glucose metabolic process (GO:0006006), glycolipid biosynthetic process (GO:0009247), hexose metabolic process (GO:0019318), PPAR signaling pathway (hsa03320) and fatty acid metabolism (hsa01212) (Fig. 3b, Additional file 6).
Functional enrichment analysis of differentially expressed lncRNAs (DELs)
A total of 696 and 706 DELs were identified in the GRI–MGI and RSI–MGI groups, respectively (Fig. 4a). Among these DEGs, 292 were up-regulated and 404 were down-regulated in the GRI–MGI group, whereas 379 were up-regulated and 327 were down-regulated in the GRI–MGI group.
Generally, lncRNAs act in cis, as diffusion or transport to other cellular compartments renders these transcripts too dilute to perform any function [22]. Recent studies focused on potential protein-coding genes affected by lncRNAs located within 100-kb upstream and downstream regions [23]. We thus performed functional enrichment analysis of potential protein-coding genes located near the DELs to explore their functions. A total of 587 and 583 GO terms and pathway categories were significantly enriched, including biological process (BP), cellular component (CC) and molecular function (MF) (see Additional file 7). Notably, target genes of up-regulated DELs in the GRI–MGI group (i.e., lncRNAs showing higher expression at the GRI stage) were found to be primarily involved in muscle system process (GO:0003012), muscle contraction (GO:0006936), myotube differentiation (GO:0014902), regulation of muscle system process (GO:0090257), AMPK signaling pathway (hsa04152) and positive regulation of immune effector process (GO:0002699) (Fig. 4b). By contrast, target genes of down-regulated DELs in the GRI–MGI group (i.e., lncRNAs showing higher expression at the MGI stage) were enriched in amino acid activation (GO:0043038), ATPase activity (GO:0016887), mitochondrion organization (GO:0007005), mitochondrial respiratory chain complex assembly (GO:0033108) and hallmark glycolysis (M5937) (Fig. 4b). In the RSI–MGI group, target genes of up-regulated DELs (i.e., lncRNAs showing higher expression at the RSI stage) were mainly enriched in mitochondrial matrix (GO:0005759), mTOR signaling pathway (hsa04150), muscle system process (GO:0003012), PI3K-Akt signaling pathway (hsa04151) and hallmark glycolysis (M5937) (Fig. 4b), whereas target genes of down-regulated DELs (i.e., lncRNAs showing higher expression at the MGI stage) were enriched in skeletal muscle tissue development (GO:0007519), skeletal muscle organ development (GO:0060538), hexose catabolic process (GO:0019320) and cGMP-PKG signaling pathway (hsa04022) (Fig. 4b).
Dynamic expression of myogenesis genes and lncRNAs
To investigated the changes in gene expression during muscle development, we analyzed the dynamic expression counts of myogenesis related genes at the GRI, MGI and RSI stages. As shown in Fig. 5a, there was more counts per million (CPM) at the MGI stage than at the other two stages. Because of the lack of lncRNA annotation libraries, we could not directly predict the function of lncRNAs. Gene expression correlation across samples can be used as an indicator of functional coregulation [24]. We therefore performed correlation analysis of lncRNAs and myogenesis related genes downloaded from the Molecular Signatures Database (MSigDB) [25]. Intriguingly, the expression pattern of lncRNA G1430 was similar to that of myogenesis related genes (Fig. 5b).
The lncRNA G1430 was up-regulated at the MGI stage, and its expression pattern was confirmed by qRT-PCR (Fig. 5c). Additionally, the expression of 10 myogenesis related genes showed a significant correlation with that of lncRNA G1430, of which six genes (APOD, TNNT2, MYBPH, MYL3, DAPK2, RIT1) showed a significant positive correlation, while four genes (TEAD4, OCEL1, AKT2, APLNR) showed a significant negative correlation (see Additional file 8). To verify the correlation between lncRNA G1430 and myogenesis related genes, the expression of lncRNA G1430 and two myogenesis marker genes (myoD1 and myoG) was analyzed by qRT-PCR in 47 pigs (Fig. 5d), followed by Pearson correlation analysis. The results showed that lncRNA G1430 was significantly positively correlated with myoD1 (r = 0.55; P = 5.9E-05) and myoG (r = 0.43; P = 0.0029) (Fig. 5d). Based on these results, we further analyzed the sequence and function of lncRNA G1430 by bioinformatics analysis and in vitro experiments, respectively. Analysis of lncRNA G1430 using CNIT (http://cnit.noncode.org/CNIT/) suggested a low coding potential of the whole sequence (Fig. 5e), which was consistent with a classic non-coding RNA feature [26]. Subsequently, we performed the RACE assay to identify the full-length sequence of lncRNA G1430 in skeletal muscle, according to the sequence archived in the RNA-seq data. The results of RACE showed that the full-length sequence of lncRNA G1430 is 316 nt (Fig. 5f, Additional file 9). Both prediction and qRT-PCR analysis suggested that lncRNA G1430 is mainly located in the cytoplasm of skeletal muscle cells (Fig. 5g). Given that lncRNA act as a miRNA sponge via its ceRNA activity, thereby regulating the target gene expression of miRNAs [27–29], we next explored the binding of miRNAs of lncRNA G1430. The putative binding sites were identified RNAhybird-based prediction of the lncRNA sequence and miRNA seed region (Fig. 5h) and verified by the dual luciferase assay. The results showed that miR-133a significantly decreased the luciferase activity when co-transfected with miR-133a mimic and pCK-G1430-3'UTR-WT, and recovered the luciferase activity when co-transfected with miR-133a mimic and pCK-G1430-3'UTR-Mut (Fig. 5h). Thus, these results showed that lncRNA G1430 acted as a sponge for ssc_mir-133a-3p, thereby reversing the luciferase activity.
Validation of lncRNAs
Four lncRNAs (G5755, G11155, G8431 and G19619) were randomly selected for validation by quantitative real-time PCR (qRT-PCR) in three replicates, and the relative expression of all four lncRNAs determined by qRT-PCR was compared with their transformed log2(CPM+1) values determined by RNA-seq (Fig. 6a). The qRT-PCR and RNA-seq data of all four lncRNAs were consistent during muscle development. We also investigated the relative expression of lncRNAs G5755 and G8431 in eight other tissues (Fig. 6b). The results showed that both these lncRNAs, especially the lncRNA G8431, were highly expressed in skeletal muscle tissues. Together, these results demonstrate the reliability of our RNA-seq data, thus confirming the accuracy of lncRNAs identified in the present study.