Muscle Fiber Type Distribution
Hematoxylin and eosin staining was performed to characterize the structural traits of longissimus thoracis between Lantang and Landrace pigs at birth (LT1D and LW1D, respectively) and 90 postnatal days (LT90D and LW90D, respectively). The number of fibers per unit area and average cross-sectional area of myofibers were determined (Fig. S1). We found that the structural traits showed large variations between different growth stages (P < 0.01). To explain, the number of myofibers significantly decreased between birth and 90 postnatal days (P < 0.01), and the average cross-sectional area of myofibers showed an obvious increase during postnatal development (P < 0.01). Further, Lantang pigs showed higher number of myofibers than Landrace pigs at birth (P < 0.01), but there were no significant differences in terms of the cross-sectional area of myofibers. In comparison with Landrace pigs, Lantang pigs showed lower number and cross-sectional area of myofibers at 90 postnatal days (P < 0.01). We then calculated the proportion of different muscle fiber types based on the expression of myosin heavy chain isoforms (MyHCs; Fig. 1). The proportion of MyHC I, IIa, and IIx myofibers at birth was higher than that at 90 postnatal days in both Lantang and Landrace pigs (P < 0.01), while the proportion of MyHC IIb myofibers was higher at 90 postnatal days (P < 0.01). At birth, the expression of MyHC I and IIa in Lantang pigs was significantly higher than that in Landrace pigs (P < 0.01), while the expression of MyHC IIb was higher in Landrace pigs (P > 0.05). Besides, at 90 postnatal days, higher amount of MyHC I was distributed in Lantang pigs (P < 0.01), and MyHC IIb showed the opposite trend between Lantang and Landrace pigs (P < 0.01).
miRNA Expression Analysis
miRNA-seq generated 21.06 ± 1.32 million raw reads with a length of 49 nucleotides from each library. After filtering, approximately 18.68 ± 1.33 million clean reads were obtained, accounting for 89.12% ± 3.47% of total raw sequences (Table S1A). The clean sequences were then annotated and assigned to 470.78 ± 86.52 thousand unique tags in each library by alignments to Rfam and porcine-specific sequences within miRBase, Repeat and Reference mRNA databases (Table S1B). We observed that only 3.16% ± 0.41% (11,464 ± 415 per library) unique reads belonged to known porcine miRNAs, and these unique reads represented 15.27% ± 0.69% (2,871,531 ± 269,401 per library) of total clean sequences. Using the miRDeep2 algorithm, 321, 311, 315, and 316 known miRNAs were identified in LT1D, LT90D, LW1D, and LW90D libraries, respectively, and 302 miRNAs were common across all samples. The 10 most highly expressed miRNAs in each library accounted for 71.26% ± 1.15% of the total count of all identified miRNAs, and six miRNAs (miR-1, miR-206, let-7a, let-7c, miR-10b, and let-7f) were found across all libraries. Of them, ssc-miR-206 showed the highest expression level in LT1D and LW1D libraries, as well as ssc-miR-1 in LT90D and LW90D libraries. We further investigated differentially expressed miRNAs between different breeds and growth stages (Table S1C). In Lantang pigs, 89 miRNAs were differentially expressed between LT1D and LT90D libraries; 34 miRNAs were upregulated and 55 were downregulated in LT90D libraries. In Landrace pigs, 100 differentially expressed miRNAs were identified; 19 miRNAs were upregulated and 81 were downregulated in LW90D libraries. At birth, only seven miRNAs were differentially expressed between LT1D and LW1D libraries, including six up- and one downregulated miRNAs in LW1D libraries, and at 90 postnatal days, 11 miRNAs were differentially expressed between LT90D and LW90D libraries, including four up- and seven downregulated miRNAs in LW90D libraries.
Transcriptome Expression Analysis
In total, 12 muscle tissue samples obtained from Lantang and Landrace pigs at birth and 90 postnatal days (in triplicate) were subjected to Illumina sequencing after rRNA depletion, which led to the generation of approximately 1.49 billion reads (average of 124.57 ± 0.29 million reads per sample). After quality control/trimming, 122.86 ± 0.21 million valid reads were obtained, accounting for 98.63% ± 0.19% of raw reads in each library. On alignment of all valid reads, we found that over 84.45% ± 1.46% clean reads could be successfully mapped to the porcine Sscrofa11.1 reference genome, including 78.88% ± 1.26% mapped reads with proper pair alignment (Table S2A). Transcript assemblies with StringTie revealed 138,278 isoforms across the 12 libraries, including approximately 24.98% identified candidates that completely matched Ensembl transcript regions (Table S2B). A comparison of known Ensembl transcripts revealed that 39,734, 39,909, 40,445, and 38,429 known transcripts were expressed in LT1D, LT90D, LW1D, and LW90D libraries, respectively; 42,081 transcripts existed in all libraries (Table S3A). Principal component analysis of globally expressed transcripts with fragments per kilobase of transcript per million mapped reads (FPKM) levels was performed, which indicated that muscle tissues in each group clustered together, as expected (Fig. 2A). We therefore applied the Ballgown algorithm to analyze differences in libraries between different breeds and growth stages (Fig. 2B). With normalized RPKM, there were 4,321 differentially expressed Ensembl transcripts between LT1D and LW1D libraries; 2,797 transcripts were upregulated and 1,524 were downregulated in LW1D libraries (Table S3B). Between LT1D and LT90D libraries, we detected 3,065 differentially expressed transcripts; 2,190 transcripts were upregulated and 875 were downregulated in LT90D libraries (Table S3C). Further, 4,292 differentially expressed transcripts were identified between LW1D and LW90D libraries; 1,335 transcripts were significantly upregulated and 2,957 were downregulated in LW90D libraries (Table S3D). In comparison with LW90D libraries, the expression levels of 1,365 transcripts were significantly different in LT90D libraries; 186 and 1,179 transcripts were up- and downregulated in LW90D libraries, respectively (Table S3E). In total, 6,962 unique differentially expressed transcripts were found on comparing LT1/90D, LW1/90D, LT/LW90D, and LT/LW1D, and only 498 transcripts were common (Fig. 2C). Gene ontology (GO) analysis revealed that these differentially expressed transcripts were significantly enriched (P < 0.05) in several biological processes associated with myogenesis, including skeletal muscle cell differentiation, muscle cell cellular homeostasis, positive regulation of smooth muscle cell proliferation, smooth muscle tissue development, muscle contraction, regulation of skeletal muscle satellite cell proliferation, and response to muscle stretch (Table S4A–D, Fig. 2D). Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that several differentially expressed transcripts were involved in muscle development and growth pathways, such as mTOR signaling pathway, Wnt signaling pathway, AMPK signaling pathway, and biosynthesis of amino acids (Table S4E–H). We randomly selected 10 dysregulated mRNAs (PFKM, ANKRD2, MSTN, MYOD1, SRF, IGF1, MYBPC2, LIMCH1, PFKFB1, and MEF2D; Fig. S2A) from these myogenesis-related GO terms and signaling pathways and validated their expression levels by performing quantitative real-time PCR (RT-qPCR). Between LT1D and LW1D, RT-qPCR data revealed that the expression levels of ANKRD2, MYOD1, LIMCH1, and MEF2D were significantly upregulated in LT1D, but those of PFKM, MSTN, MYBPC2, SRF, and PFKEB1 did not show a significant change. The RT-qPCR results of MSTN, ANKRD2, and SRF were inconsistent with those of RNA-seq. According to RNA-seq data, there was no significant difference in the expression level of ANKRD2, whereas the expression levels of SRF and MSTN were significantly upregulated in LT1D and LW1D, respectively. Further, in the comparison between LT1D and LT90D, the expression levels of ANKRD2, MYOD1, LIMCH1, and MEF2D were significantly upregulated in LT1D, and those of PFKM, MSTN, SRF, and MYBPC2 were significantly upregulated in LT90D; PFKFB1 was not significantly differentially expressed. RNA-seq did not reveal any significant differences in LIMCH1 expression between LT1D and LT90D. In the comparison between LW1D and LW90D, the expression levels of ANKRD2 and MEF2D were significantly upregulated in LT90D, and those of MSTN, SRF, MYBPC2, and PFKFB1 were significantly upregulated in LW90D; PFKM, MYOD1, and LIMCH1 expression showed no significant differences. According to RNA-seq data, the expression level of SRF was significantly upregulated in LW1D, which contradicted RT-qPCR results. In the comparison between LT90D and LW90D, the expression level of PFKM was significantly upregulated in LT90D and that of PFKFB1 was significantly upregulated in LW90D, but MSTN, MYOD1, MYBPC2, LIMCH1, ANKRD2, and SRF expression levels showed no significant differences. The results for MSTN, PFKFB1, and PFKM were inconsistent between RT-qPCR and RNA-seq. RNA-seq data indicated that the expression levels of PFKFB1 and PFKM did no show a significant change between LT90D and LW90D, while the expression level of MSTN was significantly upregulated in LW90D (Fig. S2B).
Identification of circRNAs
We characterized circRNA landscape and expression by performing deep RNA-seq experiments using the 12 aforementioned muscle tissue samples. In total, 52,133 circRNA candidates were identified using five different predicting algorithms (Fig. 3A); circRNA landscape differed quite radically depending on the algorithm used. To explain, 25,295, 33,283, 10,601, 38,292, and 4,751 circRNA candidates were detected by CIRCexplorer2, circRNA_Finder, CIRI, find_circ, and MapSplice algorithms, respectively (Table S5); find_circ and MapSplice exhibited the highest and lowest level of sensitivity, respectively. Only 3,352 circRNA candidates were commonly detected by all five algorithms, and these were subjected to further analyses. These circularization events were found to be produced from 1,745 hosting transcript loci, including 712 transcripts that generated multiple circRNA candidates (Table S6). With normalized back-splice junction reads, we analyzed significant differences in circRNA candidates across four comparisons: LT/LW1D, LT1/90D, LW1/90D, and LT/LW90D (Fig. 3B). Only three differentially expressed circRNA candidates were found between LT1D and LW1D libraries, and all three of them were significantly upregulated in LW1D library. Further, 39 and 38 circRNA candidates were differentially expressed between LT1D and LT90D libraries and between LW1D and LW90D libraries, respectively (Fig. 3C). Interestingly, 24 differentially expressed circRNA candidates were differentially expressed between LT90D and LW90D libraries, and all of them were downregulated in LW90D libraries.
Construction of circRNA-associated-ceRNA Networks
The expression of circRNAs potentially plays a key role in physiological and pathological conditions by regulating endogenous RNA targets [25]. We therefore performed Pearson correlation analysis to assess the association between differentially expressed circRNAs and mRNAs in each comparison (Fig. S3), which revealed 8, 187, 456, and 69 significant interactions in LT1D vs. LW1D, LT1D vs. LT90D, LW1D vs. LW90D, and LT90D vs. LW90D, respectively. Few circRNA candidates have been reported to directly modulate the transcription of their parent genes [26]. Herein we found that only circANKRD2, derived from exons 3 and 4 of ANKRD2, was positively correlated with its linear counterpart at the expression level between LW1D and LW90D, suggesting the involvement of circANKRD2 and ANKRD2 in myogenesis. In addition, it has been found that endogenous circRNAs serve as miRNA sponges to consequently repress the function of their targets [27]. This prompted us to predict shared miRNA-binding sites between differentially expressed circRNAs and mRNAs (Table S7A–G) and further analyze circRNA-miRNA-mRNA ceRNA networks. We identified 777, 855, and 22 convincing ceRNA interactions in LT1D vs. LT90D, LW1D vs. LW90D, and LT90D vs. LW90D, respectively (Table S7H–J); the number of putative interactions per miRNA markedly varied, ranging from 1 to 51 miRNA-associated ceRNA networks. We observed that the highly expressed circKANSL1L, circKANSL1L_2, circKANSL1L_3, circKANSL1L_4, and circKANSL1L_5 participated in 279 ceRNA transcriptional regulatory axes, including a total of 27 unique myo-miRNAs and 30 special myogenes. As evident from Fig. 3D, LT1D and LT90D comparison revealed 31 up–down–up regulation patterns: circKANSL1L, circKANSL1L_2, and circKANSL1L_3 were upregulated in LT90D and could sponge miR-128, miR-130a, miR-133b, miR-142-3p, miR-19a, miR-19b, miR-432-5p, miR-7142-3p, and miR-885-5p to significantly upregulate ATF3, CFL2, COPS2, DLG1, ITGA8, MSTN, MYOD1, PFKM, and TIPARP expression (Table S7H). On the contrary, 14 down–up–down regulation patterns were identified on comparing LW1D and LW90D: circKANSL1L_4 was downregulated in LW90D and could sponge miR-130a, miR-19a, miR-19b, miR-299, miR-376a-3p, miR-487b, and miR-493-5p to significantly downregulate ACTN2, COPS2, FBXO40, FOXN2, MYBPC1, SCN7A, TPM3, and TPM4 expression (Table S7I).
Characterization of Myogenesis-related circRNAs
To verify the circular structure of circRNAs, differentially and highly expressed circRNA candidates that were correlated with myogenes were selected for further analyses. cDNA was amplified using a pair of divergent primers, which led to the identification of nine circRNAs: circPFKFB1, circKANSL1L-3, circLIMCH1, circKANSL1L, circ4082, circKANSL1L-2, circMYBPC2, circMYBPC2-2, and circNR1H3 (Fig. S4A). Sanger sequencing further verified their head–tail junction structure (Fig. S4B).
Furthermore, on comparing the homology of the nine aforementioned circRNAs, we found that circKANSL1L sequence showed high homology between mice and pigs. To further verify their structure, we analyzed them in C2C12 cells using divergent and convergent primers (Fig. 4). Convergent primers could successfully amplify both cDNA and genomic DNA, but divergent primers could only amplify cDNA (Fig. 4A). On RNase R digestion and RT-qPCR of the three circRNAs and their host genes, we found that there was no significant change in circRNA expression levels between the RNase R treatment and control groups. However, the mRNA expression levels of the host genes and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were significantly different (P < 0.01, Fig. 4B). These findings suggested that the structure of circKANSL1L, circKANSL1L-2, and circKANSL1L-3 was indeed circular.
Inhibition of C2C12 Cell Proliferation by circKANSL1L
To elucidate the role of circKANSL1L in myogenesis, we constructed the overexpression plasmid OE-circKANSL1L and designed the knockdown gene si-circKANSL1L, and RT-qPCR was performed to verify their effects. We found that OE-circKANSL1L and si-circKANSL1L significantly increased and decreased the expression of circKANSL1L in C2C12 cells, but the expression of the host gene KANSL1L was unaffected (Fig. 5A).
Besides, RT-qPCR was performed to assess the relative expression levels of the cell proliferation-related genes PCNA, Cyclin D1, and Cyclin E. After circKANSL1L overexpression, the expression level of PCNA (P < 0.05) and Cyclin D1 (P < 0.05) significantly decreased in C2C12 cells and that of Cyclin E showed the same trend (P > 0.05). After circKANSL1L knockdown, the expression level of PCNA (P < 0.01) and Cyclin D1 (P < 0.05) significantly increased in C2C12 cells and that of Cyclin E also increased, but the change was not significant (P > 0.05; Fig. 5B). When circKANSL1L expression was upregulated, the protein expression levels of PCNA (P < 0.05), Cyclin D1 (P < 0.01), and Cyclin E (P < 0.05) significantly decreased, and when circKANSL1L expression was downregulated, the protein expression levels of PCNA (P < 0.05), Cyclin D1 (P < 0.05), and Cyclin E (P < 0.01) significantly increased (Fig. 5C); these results were consistent with RT-qPCR results. On transfecting C2C12 cells with empty vector and OE-circKANSL1L, cell proliferation was measured by Cell Counting Kit-8 (CCK-8) assay at 0, 12, 24, and 36 h. In comparison with the empty vector group, after circKANSL1L overexpression, absorbance (450 nm) significantly decreased at 24 h (P < 0.01); however, after circKANSL1L knockdown, absorbance (450 nm) significantly increased at 36 h (P < 0.01; Fig. 5D). In addition, our cell cycle analysis showed that when circKANSL1L was overexpressed, cells were arrested in the G1 phase, and the number of cells entering the S phase was significantly lower than that in the control group (P < 0.05). circKANSL1L knockdown promoted the progression of C2C12 cells to the S and G2 phases (P < 0.01, Fig. 5E). Altogether, these results indicated that circKANSL1L inhibited the proliferation of C2C12 cells.
Enhancement of C2C12 Cell Differentiation by circKANSL1L
RT-qPCR was performed to assess the relative expression levels of MYF5, MYOD1, Myogenin (MYOG), and MyHC (Fig. 6A). circKANSL1L overexpression significantly increased the expression levels of MYF5 (P = 0.06), MYOD1 (P < 0.05), MYOG (P < 0.05), and MyHC (P < 0.05), while circKANSL1L knockdown significantly decreased their expression levels (P < 0.05 for all). Western blotting was performed to detect MYOG and MyHC protein expression levels (Fig. 6B). The results showed the same trend as RT-qPCR results, but the data were not significant (P > 0.05). We also assessed RNA and protein expression levels of MyHC I and IIb. circKANSL1L overexpression promoted MyHC I and IIb expression, while circKANSL1L knockdown inhibited their expression, but the results were insignificant (Fig. 6C). Protein expression levels showed a consistent trend with RNA expression levels (Fig. 6D), indicating that circKANSL1L has a regulatory effect on muscle fiber type differentiation; further studies are nevertheless warranted.