Sheep skeletal muscle transcriptome analysis reveals muscle growth regulatory lncRNAs
- Published
- Accepted
- Received
- Academic Editor
- Juan Loor
- Subject Areas
- Agricultural Science, Genomics
- Keywords
- Transcriptome analysis, Mutton sheep, LncRNA, Skeletal muscle
- Copyright
- © 2018 Chao et al.
- Licence
- This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
- Cite this article
- 2018. Sheep skeletal muscle transcriptome analysis reveals muscle growth regulatory lncRNAs. PeerJ 6:e4619 https://doi.org/10.7717/peerj.4619
Abstract
As widely distributed domestic animals, sheep are an important species and the source of mutton. In this study, we aimed to evaluate the regulatory lncRNAs associated with muscle growth and development between high production mutton sheep (Dorper sheep and Qianhua Mutton Merino sheep) and low production mutton sheep (Small-tailed Han sheep). In total, 39 lncRNAs were found to be differentially expressed. Using co-expression analysis and functional annotation, 1,206 co-expression interactions were found between 32 lncRNAs and 369 genes, and 29 of these lncRNAs were found to be associated with muscle development, metabolism, cell proliferation and apoptosis. lncRNA–mRNA interactions revealed 6 lncRNAs as hub lncRNAs. Moreover, three lncRNAs and their associated co-expressed genes were demonstrated by cis-regulatory gene analyses, and we also found a potential regulatory relationship between the pseudogene lncRNA LOC101121401 and its parent gene FTH1. This study provides a genome-wide resolution of lncRNA and mRNA regulation in muscles from mutton sheep.
Introduction
As one of the most important domestic animals worldwide, sheep (Ovis aries) are raised mainly for meat and other agricultural products. To improve the meat production of existing mutton sheep breeds, research focused on the molecular mechanisms underlying sheep skeletal muscle development is of vital interest. To identify genes that affect muscle growth rates in mutton sheep, several studies on sheep skeletal muscle growth have been performed using transcriptome sequencing (Zhang et al., 2014; Sun et al., 2016a; Sun et al., 2016b; Bidwell et al., 2014). Many regulatory genes involved in sheep muscle growth and development have already been detected and reported. However, very few of these studies have focused on long noncoding RNAs (lncRNAs) in their transcriptome analyses.
In human and mouse, lncRNAs have already been reported as important regulatory factors of muscle growth and differentiation in multiple studies (Wang et al., 2015; Mueller et al., 2015; Han et al., 2015; Ballarino et al., 2015; Legnini et al., 2014). In goat and bovine, which are both closely related to sheep, lncRNAs have also been found to have crucial functions in skeletal muscle development (Zhan et al., 2016; Jin et al., 2017; Xu et al., 2017). Discovering differentially expressed lncRNAs in skeletal muscle from different sheep breeds would enable us to have a better understanding of the regulatory functions of lncRNAs in sheep muscle growth. Obviously, identification of relevant lncRNAs would improve our understanding of the regulatory mechanisms involved in sheep muscle growth and meat production. Recently, several studies have reported the detection and analysis of lncRNAs from multiple tissues in sheep (Bidwell et al., 2014; Yue et al., 2016; Bakhtiarizadeh et al., 2016), though similar studies on sheep skeletal muscle are still needed. In our previous report, we detected lncRNA transcripts expressed in sheep skeletal muscle (Chao et al., 2016). However, in the absence of samples and reference genome annotation, the biological functions of the lncRNAs remain unknown.
In this study, we used publicly available RNA-seq data from four different sheep skeletal muscle sequencing projects to identify and functionally predict lncRNAs. Among the four projects, two were designed to study the transcriptome differences between high production mutton sheep (Dorper sheep and Qianhua Mutton Merino sheep) and low production mutton sheep (Small-tailed Han sheep), which is consistent with our research purpose and was applied when detecting differential lncRNA expression. As a cross between Black headed Persian and Dorset, Dorper is a composite breed derived from South Africa (Milne, 2000). As a world-famous mutton sheep breed, Dorper is well known for its good muscle conformation for producing a desirable carcass. The Qianhua Mutton Merino sheep is a new breed reported by Sun et al. (2016a) that showed better meat performance than its parent breeds, the South Africa Mutton Merino and Northeast Fine-wool sheep. Comparison of the above two breeds showed that Small-tailed Han sheep had lower meat production in the two corresponding studies (Zhang et al., 2014; Sun et al., 2016a). However, with extremely high litter sizes (2.61) and good meat flavour performance, Small-tailed Han sheep showed high value in mutton sheep breeding (Tu, 1989). In this study, the top priority was to search for lncRNAs that may affect muscle growth, which would provide better knowledge on the mechanisms underlying muscle development and improve meat production in low production breeds such as Small-tailed Han sheep. This study significantly advances knowledge regarding sheep skeletal muscle lncRNA expression and will also provide the basis for further mutton sheep breeding.
Materials and Methods
Data acquisition and filtering
No new transcriptome sequencing datasets were generated in this study. RNA-seq datasets used for differential expression analysis were downloaded from the NCBI Sequence Read Archive database (SRA) under accession numbers SRP017799 and SRP080149. Detailed animal and gender information can be found under the BioProject accessions PRJNA185414 and PRJNA335752. Moreover, 2 additional RNA-seq datasets downloaded from the European Nucleotide Archive database (ENA) under accession numbers ERP005642 and SRP031629 were used for correlation analysis. Detailed animal and gender information could be found under BioProject accessions PRJEB6169 and PRJNA223213.
Before transcriptome mapping and assembling, we performed a filtering process with the raw sequencing reads. Using Trimmomatic-0.36 (Bolger, Lohse & Usadel, 2014), adapters, leading and trailing low-quality bases, and n bases (below quality 3) were removed. Then, all the reads were scanned with a 4-base wide sliding window and cut when the average quality dropped below 15. Finally, an average quality ≥ 20 and a minimum length of 36 bp were selected for the threshold.
Reads mapping and gene annotation
Using STAR (Dobin et al., 2013), we mapped our clean reads with the sheep reference genome Oar_v4.0 (ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Ovis_aries/all_assembly_versions/GCF_000298735.2_Oar_v4.0/GCF_000298735.2_Oar_v4.0_genomic.fna.gz). To detect the expression of the pseudogene loci, we created GTF files with genome-mapped BAM files in Cufflinks (Trapnell et al., 2010) and merged the GTF file with the sheep reference annotation GFF file (ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Ovis_aries/all_assembly_versions/GCF_000298735.2_Oar_v4.0/GCF_000298735.2_Oar_v4.0_genomic.gff.gz). All “gene_biotype=lncRNA” categories in the gene annotation GTF file were considered to be lncRNAs. Remapping, annotation and expression detection of the genes and lncRNAs were performed with STAR and RSEM (Li & Dewey, 2011).
Transcriptome data normalization and differentially expressed genes identification
For the eight samples from PRJNA185414 and PRJNA335752, expression data normalization was performed with the R package RUVSeq (Risso et al., 2014). All genes and lncRNAs were filtered by requiring more than 2 reads in at least 5 samples. Then, we normalized the expression matrix using upper-quartile (UQ) normalization from EDASeq (Bullard et al., 2010). Finally, according to the RUVSeq-recommended parameters, 5000 in-silico empirical negative control genes were used for unwanted variation factor estimation and expression data normalization. Relative log expression analysis (RLE) and principal components analysis (PCA) were performed according to the RUVSeq manual.
Normalized gene count data were used for differential expression analysis using DESeq2 (Love, Huber & Anders, 2014). DESeq2 used the Wald test for differential expression hypothesis testing (Love, Huber & Anders, 2014). The Wald test P-values were then independently filtered under the null hypothesis (Bourgon, Gentleman & Huber, 2010) and adjusted for multiple testing using the procedure of Benjamini & Hochberg (1995). The significant differentially expressed genes were declared at a log2-fold change ≥0.8 and a false discovery rate (FDR) <0.05.
Comparative sequence analysis
To identify differentially expressed lncRNAs that were already annotated in mode species (human and mouse), NCBI blast 2.7.1 was used for comparative sequence analysis. Using BLASTN, we compared our 39 lncRNAs to the Refseq RNA and non-Refseq RNA (human and mouse). With an e-value <0.001 as the threshold, we selected the following criteria for the comparative analysis: sequence mapping identity 75% in a covered region and 100 nt for hitting length.
Correlation analysis for lncRNAs with protein-coding genes and functional annotation
The correlation analysis between the lncRNAs and protein-coding genes was performed with expression data from all 26 samples from the four RNA-seq projects. The Pearson correlation test was used to estimate the co-expression relationships between lncRNAs and protein-coding genes. Moreover, the P value of the correlation coefficient was estimated. Using the MultiExperiment Viewer (MeV) (Saeed et al., 2003), hierarchical clustering was performed with the correlation r-values.
Before further functional annotation with the expression data, we calculated the Pearson correlation coefficient for the differentially expressed protein-coding genes with differentially expressed lncRNAs. Co-expressed genes from three clusters were identified by applying the correlation r-value >0.7 and p-value <0.001 as the threshold. Genes that achieved the threshold with at least one lncRNA were used for GO and KEGG enrichment analysis. The GO terms and KEGG pathway enrichment was performed using The Database for Annotation, Visualization and Integrated Discovery (DAVID v 6.8, https://david.ncifcrf.gov/) (Huang, Lempicki & Sherman, 2009).
LncRNA and protein-coding gene interaction network
We constructed an lncRNA-mRNA co-expression network with the correlation analysis results. Additionally, a Protein–protein interaction (PPI) network between the protein-coding genes was constructed based on information from STRING v.10.0 (Szklarczyk et al., 2015) and credible interactions (combined_score ≥ 0.4) were accepted for further network analysis. Using Cytoscape (version 3.5.1) (Shannon et al., 2003), the two networks were then merged as one, and the resulting network was defined as a lncRNA-gene interaction network. Interaction degree analysis was applied with Degree Sorted Layout, and all nodes were sorted with interaction degree values. Network module analysis identification was applied with the MCODE method using the MCODE plug-in in CytoScape; the node size was selected to be proportional with the interaction degree. For the MCODE analysis, the degree cut-off was selected as 2, while a node score cut-off = 0.2 and K-Core = 2 were used for Haircut Cluster finding (Bandettini et al., 2012).
Cis-regulatory gene analysis of differentially expressed lncRNAs
The differentially expressed genes were intersected to identify the genes 50 kbp upstream or downstream of the lncRNAs. A Pearson correlation p-value <0.001 was selected as the threshold for cis-regulatory gene prediction.
Results
Normalization analysis of sheep skeletal muscle transcriptome data
In this study, skeletal muscle RNA-seq data from eight sheep (four Small-tailed Han sheep, three Qianhua Mutton Merino sheep and one Dorper sheep) were downloaded for differential expression analysis. To distinguish the samples from the different sequencing projects, we named the three Qianhua Mutton Merino sheep from PRJNA335752 as M1, M2 and M3; the three Small-tailed Han sheep from PRJNA335752 as S1, S2 and S3; the Dorper sheep from PRJNA185414 as DP; and the Small-tailed Han sheep from PRJNA185414 as SH. The downloaded raw reads were filtered using the same thresholds and the clean reads were mapped to sheep reference genome Oar 4.0 (Table S1).
Because the data we used are derived from two sequencing projects, to check the batch effect between the samples, we performed relative log expression analysis (RLE) and principal components analysis (PCA) on the eight samples according to their gene counts data. As shown in Fig. 1, the RLE boxplots (Fig. 1A) and principal component plots (Fig. 1B) reveal a clear need for between-sample normalization. Using RUVseq, the gene count data were normalized with in-silico empirical negative control genes. After normalization, the eight samples showed consistency in the relative log expression analysis (Fig. 1C). The higher production group (Mutton Merino and Dorper) and lower production group (Small-tailed Han sheep) were divided into two groups by principal component 1 (Fig. 1D).
Figure 1: Relative log expression analysis (RLE) and principal components analysis (PCA) performed by R package RUVSeq of eight samples according to their gene counts data.
Green: high production group samples. Orange: low production group samples. (A) Boxplots of un-normalized sample RLE. (B) Plots of un-normalized sample PCA. (C) Boxplots of normalized sample RLE. (D) Plots of normalized sample PCA.
Identification of differentially expressed genes and lncRNA
To identify the differentially expressed genes, a comparison was performed between the higher production group (Mutton Merino and Dorper) and the lower production group (Small-tailed Han sheep). Using the DESeq2 algorithm (log2 Fold Change >0.8 and FDR <0.05), a total of 704 genes were identified as differentially expressed. In total, 386 genes were up-regulated in higher production group, while the remaining 318 genes were up-regulated in the lower production group (Table S2). Among them, 606 were known protein-coding genes, 39 were known lncRNAs, and 59 were novel gene loci. Interestingly, the top three log2-fold change value genes were all lncRNAs (LOC101101991: 5.23; LOC106991804: 4.09; and LOC105611977: −3.25).
Among the 39 differentially expressed lncRNAs, eight were detected from the pseudogene locus, while the potential features of the other 31 lncRNAs are all unknown. Performance of a comparative sequence analysis showed that only two of these 31 lncRNAs received acceptable matching results with known lncRNAs from human and mouse. LncRNA LOC105611269 showed high similarity with mouse lncRNA Nr6a1os, and, with the second highest fold change among all the differentially expressed transcripts, lncRNA LOC106991804 (log2 Fold Change = 4.09) showed high similarity with the human lncRNA LOC285847.
Co-expression analysis between lncRNA and protein-coding genes
To reveal the genes and functions that the lncRNAs may be related to, we performed co-expression analysis between the differentially expressed lncRNAs and differentially expressed protein-coding genes. The correlation r-values and p-values are shown in Table S3.
After the Pearson correlation analysis, hierarchical clustering was performed with the correlation r-values. As shown in Fig. 2, differentially expressed lncRNAs were clustered into 4 groups based on their correlation relationship with coding genes. With r > 0.7 and p-value <0.001 as the threshold, a total of 369 genes were found to be co-expressed with 32 differentially expressed lncRNAs from clusters 1 (three lncRNAs), 2 (16 lncRNAs) and 3 (13 lncRNAs), while the other seven lncRNAs (including the two lncRNAs from cluster 4) were filtered out from further analysis (Table S4).
Figure 2: Correlation r-value clustering heat map perform by MultiExperiment Viewer showing the co-expression patterns of differentially expressed lncRNAs.
The functions of the differentially expressed lncRNAs were predicted using GO and KEGG enrichment analyses of their co-expression genes. For cluster 1, we failed to get any significant enrichment for the co-expression genes. For cluster 2, the correlated genes were significantly enriched into 14 GO terms and 13 KEGG pathways, among which, the top enriched terms and pathways were associated with metabolism (Fig. 3A). For cluster 3, the correlated genes were significantly enriched into 19 GO terms and eight KEGG pathways, among which, the top enriched terms and pathways were mostly associated with cell proliferation and apoptosis (Fig. 3B).
Figure 3: GO term and KEGG pathways achieved significant enrichment levels.
Green: biological process GO term; blue: cellular component GO term. Yellow: molecular function Go term. Red: KEGG pathway. (A) Cluster 2 lncRNA co-expression genes GO term and KEGG pathway enrichment result. (B) Cluster 3 lncRNA co-expression genes GO term and KEGG pathway enrichment result.
Construction of the lncRNA–gene interaction network
Based on the expression correlation analysis, we constructed a lncRNA-gene co-expression network with lncRNAs (from clusters 1, 2 and 3) and their correlated genes. A total of 1,206 interactions between 32 lncRNAs and 369 genes were observed. Then, we performed a Protein–protein interaction analysis with the 369 genes that were significantly correlated with lncRNAs. A total of 526 interactions between those genes were obtained from STRING. These two networks were merged as a single lncRNA–gene interaction network (Fig. 4). The merged network was constructed by 391 nodes (185 expressed higher in the higher production group and 206 expressed higher in the lower production group) and 1,726 interactions. Among the 391 nodes, 30 were lncRNAs and 361 were gene nodes. All the gene nodes and lncRNA nodes were sorted with interaction degrees (Table S5), and the top 12 interaction degree nodes are all lncRNAs, with all five pseudogene lncRNAs included.
Figure 4: Merged lncRNA–gene interaction network created with Cytoscape.
The node size was decided on the basis of the interaction degree value. Square node, Protein coding genes; Hexagon node, LncRNA nodes; Red node, Higher expressed genes in the higher production group; Blue node, higher expressed genes in the lower production group; orange node, higher expressed lncRNAs in the higher production group; purple node, higher expressed lncRNAs in the lower production group; gray edge, lncRNA-gene interaction; green edge, protein–protein interaction.
Further MCODE analysis revealed that there are 8 modules in the network, among which only the top 2 modules showed scores ≥5. We named these two modules module A (Fig. 5A) and module B (Fig. 5B). In module A, all 14 nodes, except for the GAMT gene, showed higher expression in the lower production group (Small-tailed Han sheep). Protein-coding genes from module A could be significantly enriched into 2 GO terms, poly(A) RNA binding (FDR = 0.00024) and nucleolus (FDR = 0.0023). As three of the highest interaction degree nodes, lncRNAs LOC101121401 (first highest degree node), LOC105616222 (second highest degree node) and LOC105616344 (fourth highest degree node) were all included in module A. In module B, 10 nodes showed higher expression in higher production group (Mutton Merino and Dorper), and the other five nodes were expressed higher in lower production group (Small-tailed Han sheep). Protein-coding genes from module B could be significantly enriched into two GO terms, transition between fast and slow fibre (FDR = 0.0075) and troponin complex (FDR = 0.011). We also found three lncRNAs in module B, including LOC106991804 (third highest degree node), LOC101101991 (seventh highest degree node) and LOC106991092 (ninth highest degree node). As a result, these six lncRNAs were identified as hub lncRNAs.
Figure 5: Module analysis performed by plugin MCODE of CytoScape on merged lncRNA–gene interaction network.
Square node, protein coding genes; hexagon node, lncRNA nodes; red node, higher expressed genes in the higher production group; blue node, higher expressed genes in the lower production group; orange node, higher expressed lncRNAs in the higher production group; purple node, higher expressed lncRNAs in the lower production group; gray edge, lncRNA-gene interaction; green edge, protein–protein interaction. (A) LncRNA–gene interaction network Module A. (B) LncRNA–gene interaction network Module B.
Regulatory analysis of differentially expressed lncRNA
To detect potential lncRNA cis-regulatory genes, we identified the chromosomal co-expression of genes 50 kbp upstream and downstream of the 32 lncRNAs. As a result, 3 cis-regulatory genes of 3 lncRNAs were detected (LOC106990587 with TRIM7, LOC105603392 with KLHL40, and LOC106991804 with ARMC12). These genes were considered potential lncRNA cis-regulatory genes. Among these three lncRNAs, only LOC106990587 showed higher expression in the higher production group (Mutton Merino and Dorper).
Of the 32 lncRNAs, five were pseudogene transcripts. Among them, we detected one pseudogene exhibiting co-expression with its parent gene. With a Pearson r value = 0.88, pseudogene LOC101121401 and its parent gene FTH1 were both expressed higher in the lower production group (Small-tailed Han sheep). It is worth noting that FTH1 showed the highest expression of all the differentially expressed genes, while LOC101121401 had the third highest expression among the differentially expressed lncRNAs.
Discussion
Sheep skeletal muscle gene profiling studies have had enormous impacts on our understanding of muscle growth, providing the identification of novel regulators of skeletal muscle gene expression and function and defining the pathways that interplay to promote mutton sheep production. In this study, we performed a transcriptome level analysis of lncRNAs in skeletal muscle tissues from high production mutton sheep (Dorper sheep and Qianhua Mutton Merino sheep) and low production mutton sheep (Small-tailed Han sheep). With the goal of identifying genes affecting muscle production in mutton sheep, similar experimental designs were used by both of the studies (Zhang et al., 2014; Sun et al., 2016a) that were selected as our data source. However, because the experiments were performed with different muscles (the Longissimus, a back muscle with high commercial value, and the Biceps brachii, a shoulder muscle) and read types (single-end vs. paired-end), the statistical power of our design might be weakened, though the exact degree of weakening is still unknown. To solve this problem, a normalization analysis was performed to eliminate the existing batch effect. However, the unbalanced data design and normalization might still be a weakness of our study. Furthermore, similar problems were also present among the samples used for the co-expression analysis, where six samples were from Longissimus dorsi and biceps from three Texel sheeps (PRJEB6169) and 12 samples were from undefined muscles from crossbred sheep harbouring a callipyge mutation or not (PRJNA223213). Further verification experiments on lncRNA expression and lncRNA-gene regulatory relationships will be performed to confirm the discoveries from our in silico analysis.
Compared with the two data source studies (Zhang et al., 2014; Sun et al., 2016a), the most important innovation in this research is the detection of differentially expressed lncRNAs. A number of studies (Wang et al., 2015; Mueller et al., 2015; Han et al., 2015; Ballarino et al., 2015; Legnini et al., 2014) have proven that lncRNAs can regulate muscle growth and differentiation through cis-regulatory, trans-regulatory or competing endogenous pathways, which indicates that lncRNAs could be important muscle growth regulatory factors and potential valuable molecular marker regions for mutton sheep breeding. In our previous report (Chao et al., 2016), we performed a novel lncRNA detection analysis with sequencing data from Zhang et al. (2014), though the functions of the identified lncRNAs remain unknown. It is worth noting that the number of differentially expressed lncRNAs detected in this study is far lower than our previous report (Chao et al., 2016), which may be caused by the different methods and the statistical power. Unlike protein-coding genes, the functions of the lncRNAs could not be annotated or predicted with functional enrichment analysis. In this study, to confirm the functions of the differentially expressed lncRNAs, we applied co-expression analysis to seek potential related genes. GO and KEGG pathway enrichment revealed that the co-expression genes primarily include those known to be related to muscle development, metabolism, cell proliferation and apoptosis. These findings suggest a role for lncRNAs in the growth of sheep skeletal muscle.
Although some pseudogenes have been reported as nonfunctional, ‘dying’ genes (Zhu et al., 2007), others have been found to be transcribed as noncoding RNAs (Emadi-Baygi et al., 2017). Pseudogene RNAs can sequester microRNAs, RNA-binding proteins or translation machinery (Poliseno, 2012), as well as produce natural antisense transcripts, which fine-tune the expression of the corresponding sense transcripts at the epigenetic, transcriptional or posttranscriptional levels (Emadi-Baygi et al., 2017; Morris, 2009). In this study, we specifically examined the correlation between pseudogene lncRNA expression and their corresponding genes. Among them, only LOC101121401 was co-expressed with its corresponding gene, FTH1. A major function of FTH1 is the storage of iron in a soluble and nontoxic state (Arosio, Ingrassia & Cavadini, 2009); FTH1 was also found to be associated with several diseases (Xu et al., 2014; Chekhun et al., 2014; Li et al., 2015). Excessive iron accumulation may cause the nitration of tyrosine residues, resulting in extensive protein damage or iron-mediated nucleic acid damage, thereby leading to muscle damage (Bian et al., 2003; Xu et al., 2012). Therefore, although there is still no direct evidence, FTH1 may regulate muscle growth through its effect on iron homeostasis. Furthermore, several studies have proven that the FTH1 gene can negatively regulate cell proliferation in human (Cozzi et al., 2000; Feng et al., 2012). In this study, FTH1 showed the highest expression among all the differentially expressed genes, which is consistent with another report on the high expression of FTH1 in skeletal muscle (Polonifi et al., 2010). This gene has also been detected with multiple pseudogenes, though their functions are still unknown. Displaying a similar expression pattern, LOC101121401 may regulate the expression of FTH1 and affect muscle growth. However, since the sequence characteristics of pseudogenes are highly similar to protein-coding genes, we cannot completely rule out possibility of an incorrect mapping between them. More experiments and evidence will be required to confirm their accuracy and regulatory relationships.
For the 31 differentially expressed non-pseudogene lncRNAs, only two lncRNAs (LOC105611269 and LOC106991804) showed high similarity with human or mouse lncRNAs. However, to our knowledge, there are no reports on the functions of these lncRNAs to date. LOC285847 has been reported as a sensitivity molecular signature for proliferative diabetic retinopathy (Pan et al., 2016). In our research, LOC106991804 showed higher expression in Small-tailed Han sheep muscle and co-expressed with 124 genes, and its co-expression genes were enriched for AMPK signalling, FoxO signalling, HIF-1 signalling, cell cycle and p53 signalling pathways. It is worth noting that LOC106991804 has been detected as a hub lncRNA in network analysis, and it has also shown a potential cis-regulatory relationship with the ARMC12 gene. These results indicate that LOC106991804 could play important roles regulating apoptosis in skeletal muscle, though its high expression may be negatively correlated with muscle growth.
Further study of the merged network with MCODE revealed eight modules, among which two top score modules were selected and analysed based on GO enrichment. For the first module, module A, most protein-coding genes were significantly enriched into two RNA editing-related GO terms, while their related lncRNAs were expressed higher in Small-tailed Han sheep. In the next module, module B, only three genes (TNNT1, TNNI1 and TNNC1) were significantly enriched into two muscle contraction-related GO terms. As slow skeletal muscle troponin genes, TNNT1, TNNI1 and TNNC1 were reported with consistent expression patterns in sheep (Sun et al., 2016b); they have also been reported as important genes for maintaining slow myofibres (Pierzchala et al., 2014). These results showed that genes in module B appear to be correlated with muscle type, and their only related lncRNA in this module, LOC106991092, could be related to skeletal muscle type determination or maintenance.
An important feature of noncoding RNAs is the capacity for cis-regulation (Guttman et al., 2009). It has been reported that lncRNAs function through cis-regulation of nearby protein-coding genes as a common mechanism (Nagel et al., 2014). In the present study, three lncRNAs showed potential cis-regulation with nearby coding genes. Among them, only LOC106990587 and its cis-regulation gene TRIM7 showed higher expression in the high production mutton sheep. TRIM7 was first identified as glycogenin interacting protein (GNIP) and reported to have a high expression in skeletal muscle (Skurat et al., 2002). TRIM7 also mediates c-Jun/AP-1 activation through Ras signalling and showed ubiquitin ligase activity towards RACO-1 (Chakraborty et al., 2015). It can be inferred that the cis-regulation between LOC106990587 and TRIM7 may increase cell proliferation and promote the growth of skeletal muscle. The other two lncRNAs, LOC105603392 and LOC106991804, showed higher expression in low production mutton sheep. For LOC106991804 and its cis-regulation gene ARMC12, we could not find any relevant studies. Their differential expression may indicate that they play a role in skeletal muscle growth, though more tests are required to reveal their detailed functions. The KLHL40 gene, which is predicted to be a cis-regulation gene of LOC105603392 in this study, was reported as a muscle-specific transcript gene locus (Garg et al., 2014). A number of studies have verified that KLHL40 is a key regulatory control gene in sarcomere thin filament growth (Chen et al., 2016; Seferian et al., 2016; Winter et al., 2016), though its impact on muscle production is still unknown. Our results indicate that LOC105603392 may affect muscle growth through the cis-regulation of KLHL40.
Overall, our computational analysis revealed a number of potential muscle growth/ production-related lncRNAs in mutton sheep. These findings will provide novel insights into transcriptional level research on mutton sheep skeletal muscle systems.
Conclusions
In summary, a total of 39 differentially expressed lncRNAs were detected in this study. Subsequent bioinformatics analyses revealed that 29 of these lncRNAs were associated with muscle development, metabolism, cell proliferation and apoptosis. Six lncRNAs were detected as hub lncRNAs, and four lncRNAs showed potential regulatory relationships with specific genes. Our research therefore provides the basis for further functional studies focused on the roles of lncRNAs in sheep skeletal muscle growth and mutton sheep production.
Supplemental Information
Statistic summary results of 26 samples from 4 sequencing projects
Sequencing project and Biosample information were collected from NCBI SRA database. Data file filtering was performed with Trimmomatic-0.36.
Differentially expressed genes and lncRNAs
Differential expression analysis was performed with DESeq2.
Correlation r-value and p-value between differentially expressed lncRNA and differentially expressed protein coding genes
Sheet 1: Pearson correlation r-value matrix. Sheet 2: Pearson correlation p-value matrix.
Information of lncRNAs and co-expressed genes
Differentially expressed lncRNAs from Cluster 1, 2 and 3 were contained in. Cis-regulatory genes were identified with 50 k bp upstream or downstream of the lncRNAs and a Pearson correlation p-value <0.001.
Network node interaction degree list
Interaction degree and MCODE score of 391 nodes were obtained with Cytoscape 3.5.1. All nodes were ranked with interaction degree.