Background The difference between the skeletal muscle growth rates of Western and domestic breeds is remarkable, but the potential regulatory mechanism involved is still unclear. Numerous studies have pointed out that long intergenic noncoding RNA (lincRNA) plays a key role in skeletal muscle development. This study used published Yorkshire (LW) and Tibetan pig (TP) transcriptome data to explore the possible role of lincRNA in the difference in skeletal muscle development between the two breeds.
Results Through differential expression analysis, 138 differentially expressed lincRNAs (DELs) were obtained between the two breeds, and their potential target genes (PTGs) were predicted. The results of Gene Ontology and pathway analysis revealed that PTGs are involved in multiple biological processes and pathways related to muscle development. The quantitative trait loci (QTLs) of DELs were predicted, and the results showed that most QTLs are related to muscle development. Finally, we constructed a co-expression network between muscle development related PTGs (MDRPTGs) and their corresponding DELs on the basis of their expression levels. The expression of DELs was significantly correlated with the corresponding MDRPTGs. Also, multiple MDRPTGs are involved in the key regulatory pathway of muscle fiber hypertrophy, which is the IGF-1-AKT-mTOR pathway.
Conclusions In summary, multiple lincRNAs that may cause differences in skeletal muscle development between the two breeds were identified, and their possible regulatory roles were explored. The findings of this study may provide a valuable reference for further research on the role of lincRNA in skeletal muscle development.
In the past few decades, pig breeding has been focused on increasing lean meat ratio and growth rate, resulting in a decline in the quality and flavor of pork. However, the current increase in pork quality requirements has prompted breeders to increase the growth rates of pigs while ensuring the quality and flavor of pork.
Tibetan pigs, also known as ginseng pigs, have less subcutaneous fat, more lean meat, higher amounts of protein and amino acids, and better taste than Yorkshire pigs. Tibetan pigs are particularly popular in the high-end market. The price of Tibetan pigs is at least five times the prices of other commercial pork varieties in the Chinese market. However, compared with traditional commercial pig breeds, such as Yorkshire, Tibetan pigs have slower growth rates and lower reproductive rates. At 12 months old, Tibetan pigs usually weigh approximately 25 kg. As a typical Western breed, Yorkshire pigs have a high lean rate, rapid muscle growth, and heavy body weight. Therefore, understanding the difference in skeletal muscle growth and development between the two breeds is beneficial to the genetic improvement of pigs in the future.
In pig breeding, the growth and development of skeletal muscles directly affect the quantity and quality of animal meat. Related studies have shown that skeletal muscle growth is affected by the number, size, and type of muscle fibers. and the number and size of muscle fibers are closely related to the tenderness of pork[5, 6].
Long noncoding RNA (lncRNA) was once widely regarded as a transcriptional “noise.” However, its potential function has attracted attention in recent years. Numerous studies have shown that lncRNA plays an important regulatory role in many key biological processes, such as metabolism, immune response and apparent modification. Moreover, lncRNA plays a key role in muscle development; for instance, the mouse lncRNA H19, promotes muscle differentiation by participating in miRNA coding. Tens of thousands of lncRNAs have been identified in the pig genome by high-throughput sequencing technology, but little is known about their functions. lncRNA can affect the expression of nearby genes through cis-acting mechanisms and performs a series of functions in the whole cell through trans-acting mechanisms.
Here, DELs were obtained through pipeline analysis and differential expression analysis, and potential target genes in cis were predicted. The role of DELs in the regulation potential target genes (PTGs) was explored by performing Gene Ontology and pathway analysis on the PTGs. Meanwhile, the quantitative trait loci (QTL) of DELs were analyzed for the prediction of the potential function of DELs from other perspectives. The co-expression network of MDRPTGs and their corresponding long intergenic noncoding RNAs (lincRNAs) was constructed and used in exploring the potential role of DELs in skeletal muscle development. Finally, a comprehensive analysis was performed, and the possible role of DELs in the difference in skeletal muscle development between the two breeds was discussed. The aim of this work is to provide a valuable reference for the study of lncRNA related to porcine skeletal muscle development.
Six published RNA-seq data (Figure.S1,Table 1) sets of two pig breeds (Tibetan and Yorkshire) were used in identifying lincRNAs that may cause phenotypic differences between the breeds. As shown in Fig. 1A, potential lincRNAs were obtained on the basis of this pipeline. Approximately 245.05 of 258.43 million reads were mapped to the pig reference genome (Sus scrofa 11.1) by HISAT2. The mapped reads of each data were assembled into one set of transcripts with StringTie, and all of the transcripts from six data were merged into a nonredundant transcript set. A total of 828 transcripts, which were > 200 bp intergenic transcripts with more than two exons, were obtained from this pipeline. Finally, three different methods namely, CPC, HMMER, and BLAST were used for the assessment of the coding potential of the transcripts. A total of 361 potential lincRNAs were obtained. In addition, 53 of the potential lincRNAs had no overlaps with currently annotated coding or noncoding transcripts (Figure.1B). These lincRNAs are distributed in all chromosomes except the Y chromosome (Figure.1C).
Previous studies have shown by comparing structural features that pig lincRNA is identical to the lincRNAs of other mammals (human and mouse). The pig lincRNA has fewer and longer exons than the coding gene; the lincRNA transcript, owing to their small number of exons, is shorter than the coding gene. Thus, the present study compared the difference in exon number (Figure.2A), length (Figure.2B), and length of transcription (Figure.2C) between lincRNA (known lincRNAs and novel lincRNAs) and protein coding genes in this data, consistent with previous reports. The accuracy of lincRNAs obtained from the pipeline of this study was confirmed.
The mammalian genome is universally transcribed and encodes thousands of lincRNAs distributed throughout the genome, which are less conserved and have low expression levels[14, 15]. The present study compared the average expression levels of 361 lincRNAs (known lincRNAs and novel lincRNAs) and protein coding genes from six samples to investigate whether this expression pattern is also present in pigs. The results showed that the average expression level of lincRNAs (known lincRNAs and novel lincRNAs) in pigs is generally lower than that of genes encoding proteins (Figure.3A). In order to study the lincRNA that may cause phenotypic differences between the two breeds of pigs (Yorkshire pigs and Tibetan pigs). DeSeq2 in the R package was used to perform differential expression analysis on the two breeds of pig samples on the basis of expression levels. Between the two breeds, 66 of the 138 DELs of Yorkshire pigs were upregulated and 72 were downregulated (Figure.3B). Between the two pig breeds, 326 of the 682 differentially expressed coding genes identified were downregulated and 356 were upregulated (Figure.S2).
Given that lncRNA can silence or activate cis-gene expression, it can act on neighboring genes at lncRNA sites.Target genes (Table.S1, the methods of target genes prediction refer to methods) in the range of 100 kb upstream and downstream of the DELs position were searched. Conducting the online GO and pathway analysis (Table.S2) through Metascape to explore the functions of target genes that may be regulated by DELs. The results showed that 409 PTGs were significantly (P < 0.05) involved in 155 biological processes and 21 pathways. Many biological processes and pathways involved in muscle development (Figure.4A). In addition, we found that the genes in the pathways related to muscle development differently expressed between the two species (Figure.4B). It is speculated that the differential expression of target genes may be related to the differential expression of linRNA between the two species.
The functions of DELs were further explored by performing QTL mapping analysis (Table.S3) after the prediction of the target genes of DELs. The results indicated that approximately 37% of QTLs are associated with muscle growth and development (Figure.5A). At the same time, we calculated the chromosome distribution of QTLs associated with muscle development were mainly distributed on chromosomes 4 and 6 (Figure.5B). Interestingly, the top 10 QTLs were associated with muscle development, and were mainly concentrated on the average back fat thickness QTL, waist muscle area QTL and body weight QTL (Figure.5C). This result confirms the potential mechanism of DELs associated with muscle growth and development.
To further explore the potential role of DELs in muscle development. We collected PTGs from biological processes and pathways involved in muscle development (Table.S4). Based on the expression levels of these PTGs and the corresponding lincRNA, we construction the DELs- MDRPTGs co-expression network by using Cytoscape_3.6.1  (Fig. 6). It was found that 34 of the 138 DELs may regulate PTGs which associated with muscle development, and we found that 25 of 34 DELs upregulate their target genes. In addition, there are eight DELs corresponding to more than two MDRPTGs, two MDRPTGs correspond to multiple DELs. Therefore, we speculate that there are some regulatory mechanisms related to muscle development between MDRPTG and DEL.
In the PTG prediction section, we predicted 409 PTGs corresponding to 138 DELs. To confirm this result, we randomly selected 5 lincRNA genes with significant positive correlation based on their expression levels. The correlation coefficients were all greater than 0.80, and the p values were less than 0.05. We performed RT-qPCR experiments on nine samples, and the results were analyzed using linear regression. The expression levels suggested that the five pairs of lincRNA genes and their PTGs are significantly positively correlated, with a correlation coefficient greater than 0.80 and p value of less than 0.05. The experimental results of RT-qPCR showed that the results of the two datasets are in good agreement, further improving the reliability of the present study (Fig. 7).
The skeletal muscle is the largest organ in animals. In pigs, skeletal muscles have important economic significance for production, and understanding the development of skeletal muscles is important for improving productivity and meat quality. Previous studies have detected thousands of lincRNAs in skeletal muscle, but only a few number of lincRNAs have been characterized. In this study, 361 potential lincRNAs were identified on the basis of the designed pipeline, of which 53 lincRNAs were novel lincRNAs. Thereafter, 138 DELs were statistically obtained using the DeSeq2 software, of which 66 DELs were upregulated in the LW group and 72 were downregulated. The transcription length, exon length, number and expression characteristics of potential lincRNA and coding genes were compared and found to be consistent with the characteristics of other mammalian lincRNAs. The accuracy of our pipeline identification lincRNA was verified from the other hand.
The PTGs of DELs were predicted, and gene ontology analysis was used to understand the biological processes that PTGs may participate in to predict the possible role of corresponding DELs. Multiple biological processes are associated with skeletal muscle development. Such processes include muscle structure development, muscle organ development, muscle cell differentiation, and skeletal muscle tissue development. DELs may indirectly affect the development of skeletal muscle by regulating its target genes. In general, multiple genes need to cooperate with one another to play the ultimate biological function; thus, Kyoto Encyclopedia of Genes and Genomes pathway analysis was performed on PTGs. Path-based analysis is used to understand the possible role of DELs. Multiple genes were speculated to participate in the cAMP signaling pathway, the mTOR signaling pathway, and the Rap1 signaling pathway. QTL mapping analysis was performed to improve the credibility of DEL potential functional prediction. The results show that 1137 of the predicted 3018 QTLs are associated with skeletal muscle development. The proportion of loin muscle area QTL and body weight QTL is the highest. The muscle development-related PTGs were studied by generating the MDRPTG–DEL co-expression network. We then focused on muscle development related PTGs. Generating the MDRPTG–DEL co-expression network.
Based on the above analysis results, we comprehensive analysis found that multiple DELs may participate in the IGF-1-Akt-mTOR signaling pathway by regulating their PTGs. The growth of skeletal muscle depends on muscle fiber hypertrophy, and the size of muscle fibers is increased when the rate of protein synthesis is higher than the rate of degradation. Under normal physiological conditions, the IGF-1-Akt-mTOR pathway plays a key regulatory role in skeletal muscle protein synthesis[19, 20]. Interestingly, multiple PTGs are related to this pathway.
Mitochondrial calcium unidirectional transporter (MCU) is a highly selective channel for Ca2+ transport into the mitochondria. Mammucari et al. reported that MCU participates in IGF-1-Akt-mTOR signaling by increasing Ca2+ level in the mitochondria, activating the PGC-1α4, which is a transcriptional coactivator; the IGF-1 gene is activated through the PGC-1α4, leading to muscle hypertrophy[21, 22]. In the present study, our analysis showed that compared with the TP group, the LW group had higher MCU expression level (Figure.S3), which may be associated with the growth characteristics of Yorkshire pig breeds. More importantly, we found that DEL-MSTRG.8035 is positively related to the expression of MCU and highly expressed in the LW group. Insulin-like growth factor 2 (IGF-2) is a maternal blotting growth factor that regulates prenatal skeletal muscle development. It can be involved in the IGF1-Akt-mTOR signaling pathway by activating the IGF1 receptor. A significant positive correlation exists between the DEL-MSTRG.12010 and IGF-2, which were upregulated in the LW group. Interestingly, MSTRG.12010 was significantly negatively correlated with the troponin T-3 (TNNT3) gene. Wang et al. It is predicted that TNNT3 can regulate muscle growth and muscle fibers. TNNT3 is an important part of pig skeletal muscle filaments, that can affect the taste and tenderness of pork[26, 27]. Its expression level was low in the LW group. This may, on the other hand, find the reason for the decrease in meat quality as the growth rate increases. In addition, a potential target gene, serum response factor (SRF), plays an important role in controlling muscle fiber hypertrophy[28, 29]. SRF can control the transcription of miR-486, which as a potential regulator of PI3K/Akt signal transduction in muscle cells, can phosphorylate Akt and activate the IGF-1-Akt-mTOR signaling pathway, leading to muscle fiber hypertrophy. The co-expression network, suggests that DEL-MSTRG.21771 is significantly positively correlated with SRF expression and highly expressed in the LW group. It is speculated that MSTRG.21771 regulates the high expression of its potential target gene SRF, in the LW group, which possibly useful in maintaining the fast skeletal muscle rate in Yorkshire pig. PLD1 is an isoform of phospholipase D (PLD), which can stimulate phosphatidylcholine (PC) to produce phosphatidic acid (PA), which can bind to mTOR and participate in the IGF-1-Akt-mTOR signaling pathway. The substrate S6K1 of mTORC1 is phosphorylated to enhance protein translation, resulting in muscle fiber hypertrophy. Furthermore, DEL-MSTRG.6293 was positively correlated with its expression, and highly expressed in the LW group. QTL results indicated that MSTRG.8035, MSTRG12010, MSTRG21771, and MSTRG.6293 were all mapped to the QTL loci, such as Loin weight QTL, Loin muscle area QTL, and backfat above muscle dorsi QTL, which are related to muscle development. It is speculated that these DELs may be related to skeletal muscle development, and participate in the IGF-1-Akt-mTOR signaling pathway by regulating the expression of its potential target genes, thereby affecting the muscle fibrous hypertrophy process. However, the specific molecular regulation of this phenomenon remains unclear,and further studies are needed. Figure 8 shows that DELs may affect muscle protein synthesis by regulating their PTGs to participate in the IGF-1-Akt-mTOR pathway.
In conclusion, Given that native livestock breeds are adapted to local conditions and are not as strongly selected as commercial livestock breeds, they are important reservoirs of genetic variation. Studying lincRNAs that cause phenotypic differences between local and Western pig breeds is important to breed improvement. In this study, linRNAs that may cause skeletal muscle growth differences between Tibetan and Yorkshire pigs were identified, and a batch of novel linRNAs were discovered. Furthermore, we speculate that multiple DELs may participate in the IGF-1-Akt-mTOR pathway that affects muscle protein synthesis by regulating their target genes, and affect the size of muscle fibers. Our finding serve as valuable reference and new ideas for lincRNA research.
In this study, the sows used for RNA-seq were reared under similar environmental and feeding conditions(12). Tibet Agricultural University, under the same dietary and drinking water standards, raises Tibetan pigs (TP) and imported pig breeds, Yorkshire pigs (LW). After 60 days of fertilization, 9 embryos were taken from each group, and the longissimus dorsi muscle (LD) muscle tissue of the 12th rib was taken and frozen in liquid nitrogen to extract total RNA. We obtained six published RNA-seq data for two pig breeds (TP&LW) from the GEO database. All RNA-seq datasets were divided into two groups with 3 replicates in each group(12). The pig gene annotations were downloaded from ftp://ftp.ensembl.org/pub/release-91/gtf/sus_scrofa, and the non-redundant reference sequence database was downloaded from https://ftp.ncbi.nih. gov/blast/db/.
The experimental Yorkshire pigs were provided by the National Livestock Engineering Research Center of Huazhong Agricultural University. All Yorkshire pigs were raised under the same temperature, humidity, ventilation conditions and feeding standards. After fasting for 12 hours, nine sows with 55 days of pregnancy were randomly selected and euthanized by electric shock and rapid bleeding. Then take the longissimus dorsi muscle of embryo and save it in liquid nitrogen for later use.
We use FastQc to evaluate the quality of sequencing reads in the data. Low quality reads were removed using Trimmomatic (version 0.3.2) with the default parameters. The clean data obtained were aligned to the Sus scrofa genome (SusScrofa11.1) from University of California Santa Cruz (UCSC) using the HISAT2 (version 2.0.1) default parameters. Sorting mapped reads and remove duplicates via SAMtools(version 0.1.19). In addition, we assembled the read map using the default parameters of StringTie (version 1.2.2). At the same time, we set the StringTie's "-G" option for the novel transcript assembly. Finally, we used the merge function of StringTie software to combine the transcription files from six samples (GTF format) into a non-redundant transcriptome file.
361 lincRNAs were screened based on established pipelines in our laboratory . The complete pipeline follows the one shown in Fig. 1A. Step 1, Transcripts representing intergenic transcripts classified as "U" were screened using the gffcompare program StringTie (version1.2.2). Step 2, based on the transcript characteristics of lincRNA, transcripts with transcript lengths greater than 200 bp and exon numbers greater than 2 were screened. Step 3, the coding potential of lincRNAs were verified using a coding potential calculator (CPC) tool, and the lincRNAs with CPC value < 0 was entered into the downstream screening program. Step 4, to improve the accuracy of lincRNA screening, we predicted the coding potential of transcripts from multiple perspectives. We translated the transcript into six possible protein sequences, which were then transcribed and compared to the Pfam database. Finally, no Pfam hit (E value < 1e-5) transcripts are retained. Step 5, Transcripts were aligned with NCBI NR and UniRef90 databases using the BLASTX program, and transcripts of similar proteins in known proteins were filtered (E value < 1e-5). Step 6, to minimize the chance of false positives, transcripts that were not expressed (FPKM values) in all samples were removed.
We selected 45788 "gene_biotype = protein_coding" transcripts from the pig's genome annotation file (SusScrofa 11.1) to define them as transcripts of protein coding genes. In addition, we use the "blastn" instruction to divide lincRNA into known lincRNA and novel lincRNA. We then identified and compared the transcript lengths, exon lengths, exon numbers, and FPKM averages for these three categories.
We used the counting software HTSeq to count the number of reads in the six samples, and then divided the six samples into two groups according to the variety, namely TP and LW, and compared them in R using the "Deseq2" package. The gene of |log2FoldChange|>1, padj < 0.05 is a differentially expressed gene. Then, it is calculated by taking the intersection of the potential lincRNA obtained from the pipeline to obtain the differentially expressed lincRNA(DELincRNA), at the same time, it is intersected with the protein coding gene expressed in the sample to obtain a differentially expressed protein encoding gene.
In accordance with previous studies. Using bedtools (version 2.17.0) to search for mRNAs in the 100 kb upstream and downstream of the lincRNA locus, and the Pearson correlation coefficient of DELs and mRNAs was calculated in R. Finally, mRNA with a correlation coefficient greater than 0.9 was identified as a cis potential target of DELs.
Due to the limitations of pig genome annotation, this study included background human orthologous genes. After transforming the pig gene into a human gene in the Ensembl website, gene ontology (GO) and the Kyoto Gene and Genomic Encyclopedia (KEGG) pathway enrichment analysis were performed in Metascape. Subsequent selection of p value less than 0.05 is a valid result.
In this study, the pig QTL annotation file was downloaded from the animal QTL database, and the location information of DEL was proposed in the non-redundant transcription file according to the ID of DEL. After that, we performed QTL mapping on DEL using the BEDTools (version 2.17.0).
We used the longissimus dorsi muscle from nine 55-day-old embryos of Yorkshire pigs and performed RT-qPCR to verify the expression correlation between DELs and PTGs. For quantitative verification, in the first step, total RNA was extracted using Trizol reagent (Invitrogen, Life Technologies, CA, USA) and performed according to the manufacturer's instructions. The purity and concentration of total RNA were measured by a micro-spectrophotometer (Thermo, NanoDrop 2000, United States) at 260 and 280 nm. Next, cDNA synthesis was performed using the RevertAid First Strand cDNA Synthesis Kit (Thermo, Wuhan, Cat#k1622).According to the manufacturer's instructions, qPCR for DELs and PTGs detection in Roche LightCyler 480 system (Roche, Mannheinm, Germany) was performed using SYBR Green (CWBIO, Beijing, China, CW0957). Ten pairs of RT-qPCR primers were designed using the Primer 5 program (Table.S5, Table.S6). 18S rRNA is used as an endogenous control gene. The RT-qPCR data were analyzed using the 2-△△CT method and R scripts were used to perform related linear regression analysis.
Ethics approval and consent to participate
The experimental Yorkshire pigs were provided by the National Livestock Engineering Research Center of Huazhong Agricultural University. Animal breeding and slaughtering were carried out in accordance with the pre-approval guidelines of the Standing Committee of the Hubei Provincial People's Congress No. 5. And all the experimental schemes were approved by the Scientific Ethics Committee of Huazhong Agricultural University.
Availability of data and material
The sequence datasets of this paper from the GEO database. The GEO data set ID is GSE65983 (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP108727&o=acc_s%3Aa).
CCL conceived and designed the experiments and explained the data. ZYH analyzed main content of the data with the help of QQL. MXL performed the experiment. ZYH wrote the paper with the help of CCL.
We are grateful for the generous help of CCL and QQL in data analysis. We are grateful to MXL for the help of the experimental verification.
The work was supported by National Natural Science Foundation of China (NSFC, 31872322), the Fundamental Research Funds for the Central Universities (2662017PY030), National R&D Project of Transgenic Animals of Ministry of Science and Technology of China (2016ZX08006-003) and National high technology research and development plan (863, 2011AA100302).
No potential conflicts of interest were disclosed.
Consent for publication
long intergenic noncoding RNA
differentially expressed lincRNAs
potential target genes
quantitative trait loci
development related PTGs
insulin-like growth factor-1
insulin receptor substrate
phosphatidylinositol 3 kinase
mechanistic target of rapamycin in complex 1
Insulin-like growth factor 2
Serum response factor
Mitochondrial calcium unidirectional transporter
Additional file 1: Figure.S1. PCA analysis results of RNA-seq data.
Additional file 2: Figure.S2. Differential expression analysis of differential protein coding genes. the bar code represents the color scale of the log10(FPKM) .
Additional file 3: Table.S1. The statistics table of DELs target gene.
Additional file 4: Table.S2. Gene ontology and pathway analysis of adjacent genes of lincRNAs.
Additional file 5: Table.S3. Statistical table for QTLs analysis of DELs.
Additional file 6: Table.S4. The expression regulation relationship between DEL genes and their PTGs.
Additional file 7: Figure.S3. Heatmap of muscle development-related differentially expressed genes. the bar code represents the color scale of the log10(FPKM).
Additional file 8: Table.S5. The information of ten pairs of RT-qPCR primers.
Additional file 9: Table.S6. Sequence information of lincRNAs used in RT-qPCR.
Due to technical limitations, table 1 is only available as a download in the supplemental files section.