http://orcid.org/0000-0002-2141-913X
Figures
Abstract
Long noncoding RNAs (lncRNAs) play important roles in transcriptional and posttranscriptional regulation. However, the effects of lncRNAs on the meat quality of chicken hasn’t been elucidated clearly yet. Gushi chickens are popular in China because of their superior meat quality, particularly the tender flesh, and unique flavor. Gushi chickens are popular in China because of their superior meat quality, delicate flesh, and unique flavor. We performed RNA-Seq analysis of breast muscle from Gushi chicken at two physiological stages, including juvenile (G20W) and laying (G55W). In total, 186 lncRNAs and 881 mRNAs were differentially expressed between G20W and G55W (fold change ≥ 2.0, P < 0.05). Among them, 131 lncRNAs presented upregulated and 55 were downregulated. We identified the cis and trans target genes of the differentially expressed lncRNAs, and constructed lncRNA-mRNA interaction networks. The results showed that differentially expressed mRNAs and lncRNAs were mainly involved in ECM-receptor interaction, glycerophospholipid metabolism, ubiquitin-mediated proteolysis, and the biosynthesis of amino acids. In summary, our study utilized RNA-seq analysis to predict the functions of lncRNA on chicken meat quality. Furthermore, comprehensive analysis identified lncRNAs and their target genes, which may contribute to a better understanding of the molecular mechanisms underlying in poultry meat quality and provide a theoretical basis for further research.
Citation: Li D, Li F, Jiang K, Zhang M, Han R, Jiang R, et al. (2019) Integrative analysis of long noncoding RNA and mRNA reveals candidate lncRNAs responsible for meat quality at different physiological stages in Gushi chicken. PLoS ONE 14(4): e0215006. https://doi.org/10.1371/journal.pone.0215006
Editor: Juan J. Loor, University of Illinois, UNITED STATES
Received: September 25, 2018; Accepted: March 25, 2019; Published: April 9, 2019
Copyright: © 2019 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The RNA-Seq dataset supporting the conclusions of this article is available in Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) under accession numbers SRR5367021, SRR5367022, SRR5367023, SRR5367024, SRR5367025 and SRR5367026.
Funding: This work was supported by the NSFC-Henan joint grant (U1804107), the China Agriculture Research System (CARS-40-K04), the Program for Innovation Research Team of Ministry of Education (IRT16R23), and the Key Science and Technology Research Project of Henan Province (151100110800).
Competing interests: The authors have declared that no competing interests exist.
Introduction
Along with the improvement in life quality in modern society, more and more consumers paid attention to the food safety and food quality, such as the meat quality, which includes color, flavor, juiciness, fat content and tenderness. The market demands guided the direction of breeding program in poultry industry, which has made a conversion to balance the growth rate and meat quality to some extent. Meat quality is significantly affected by a series of complex factors, including breed, age, genetics, nutrition and feeding environment [1–5]. The Gushi chicken was an eminent representative with high meat quality of indigenous chicken breeds, which origins from the Gushi county of Henan province in China. The high meat quality performs in many facets. For example, high intramuscular fat (IMF) content, low shear force (SF), low drip loss, low pH decreasing, nice color and flavor. In our previous study, the laying-stage Gushi hens (55 weeks, G55W) exhibited higher serum lipid levels, IMF deposition capacity, water holding capacity (WHC), SF, intermuscular fat width (IFW) and unsaturated fatty acid (UFA) levels than the juvenile Gushi hens (20 weeks, G20W) showed, while G20W showed higher drip losses and muscle tenderness. Global DNA methylation profiles indicated that G55W exhibited higher global DNA methylation levels than G20W [2]. But phenotypic characteristics of meat quality are not only controlled by coding genes but also by noncoding genes at transcriptional regulatory level. Although many studies have investigated meat quality differences among different breeds at the mRNA level [6–8], little has been demonstrated between the different physiological stages from the same breed.
Long noncoding RNAs (lncRNAs), or mRNA-like long noncoding RNAs, are both important members of the noncoding RNA families, with lengths ranging from 200 bp to > 100 kb [9–13]. In the past, lncRNAs were thought to have no coding ability due to their lack of typical open reading frames, but recent studies have found that lncRNAs have small open reading frames (sORF) that could encode small peptides and these small peptides have important biological functions [14–16]. These studies have expanded the traditional understanding of the potential for genetic coding within the genome and have elucidated the rich variety and diversity of lncRNAs. Recent studies have shown that lncRNAs play specific roles in the development of different organs and tissues. For example, the lncRNA Bvht can affect the development of cardiac tissue in mice, as the RNAi knockout of Bvht can change the expression of cardiac specific genes, and inhibit the normal development of neonatal cardiomyocytes into mature cardiomyocytes [17]. Linc-MD1, Lnc-mg and Linc-RAM were observed to be involved in myogenesis [18–20]. LncRNA-Six1 could promote the cell proliferation, which participated in muscle growth [21]. However, little is known what role of lncRNAs played in regulating the meat quality of breast muscle from chickens at different physiological stages.
Herein, we conducted a comprehensive transcriptome profile on breast muscle tissue in juvenile and laying-stage Gushi hens. The results revealed the expression patterns of mRNAs and lncRNAs, which were important during different developmental stages in the two groups. An integrated analysis of differentially expressed (DE) lncRNAs and mRNAs was performed to elucidate the regulatory patterns of lncRNAs and their interactive network with putative target genes. Our study will help to identify the predicted interplay between lncRNAs and protein-coding genes in breast muscle, and the information generated from these predictions can be utilized in further studies of lncRNA function in chicken breast meat quality. This study may also provide insights into gene regulation and lay a foundation for genetic breeding strategies to improve meat quality.
Materials and methods
Experimental animals and sample collection
All animal experiments were performed in accordance with a protocol approved by the Institutional Animal Care and Use Committee (IACUC) of Henan Agricultural University [22]. The animals involved in this study were humanely sacrificed as necessary to ameliorate suffering.
Animals were raised in the same environment with the same ad libitum diet. Six healthy Gushi hens were selected randomly at 20 weeks old (G20W) and 55 weeks old (G55W), representing juvenile and laying-stage hens, respectively. Three individuals per stage were regarded as biological replicates. Samples were harvested from the breast muscle, snap frozen in liquid nitrogen and subsequently stored at -80°C.
RNA extraction and Solexa sequencing
Total RNA was extracted from chicken breast tissues using TRIzol reagent (Invitrogen, USA) following the manufacturer’s protocol. RNA quality was determined using the RNA Nano6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), which analyzes the integrity of the RNA. The cDNA libraries for RNA-Seq were constructed using a TruSeq RNA Sample Preparation Kit (San Diego, CA, USA) according to its instructions. First, purification of mRNA was isolated by polyA selection with magnetic oligo (dT) beads, and the purified products were randomly sheared in Fragmentation Buffer. Second, cDNA was synthesized with random primers. Then, the double-stranded cDNA was end-repaired, A-bases were added, and adaptors were ligated, and the ligation products were purified. At last, cDNA libraries were obtained by PCR enrichment. Finally, six libraries were sequenced at Novogene Co. Ltd. (Beijing, China) on an Illumina HiSeq 2000 with 2×150 paired-end (PE) reads. The RNA-Seq dataset supporting the conclusions of this article is available in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) under accession numbers SRR5367021, SRR5367022, SRR5367023, SRR5367024, SRR5367025 and SRR5367026.
Data analysis
After sequencing, the raw data were stored in fastq format (http://en.wikipedia.org/wiki/Fastq). First, clean data were obtained by removing the low-quality reads, reads containing an unknown nucleotide “N” and duplicate sequences and trimming the first 15 nt at the5’-end. The Q20, Q30 and GC content of the clean data were calculated [23]. The resulting clean, high-quality data were used for all downstream analyses. The clean, paired-end reads were aligned to the chicken genome sequence assembly using TopHat [24], and the mapped reads for each sample were assembled using both Scripture (beta2) [25] and Cufflinks [26].
LncRNA identification and prediction of target genes
To ensure the quality of the obtained lncRNAs, four criteria were used to identify the desired lncRNAs in the transcriptome assemblies: 1) transcripts with length ≤ 200 nt, exon number ≤ 2 and open reading frame (ORF) length ≥300 nt were removed; 2) Cufflinks was used to calculate the read coverage of every transcript [27], and transcripts with an FPKM value of more than 0.01 were removed; 3) the remaining transcripts were BLASTed against known chicken lncRNAs in ALDB, and only lncRNA transcripts, for which the splice sites were completely congruent between our results and those in ALDB, were immediately labeled in our results as known lncRNAs; and 4) the coding potential of the transcripts was predicted using the coding potential calculator (CPC < 1) [28], Coding-Potential Assessment Tool (CPAT) [29], Coding Noncoding Index (CNCI) [30], and Pfam [31] protein domain families to further remove coding genes, Subsequent steps were performed based on the intersection of the results obtained using the four databases. Gene transcripts within 100 kb upstream or downstream of the lncRNA were considered cis target genes. The correlation between lncRNA and target genes was analyzed using Pearson’s correlation coefficient. The lncRNA-gene pairs with absolute values of the coefficients greater than 0.95 were retained.
GO and KEGG analysis
Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed assessed using DAVID (http://david.abcc.ncifcrf.gov/). We retained the “GO terms” and “KEGG pathways” with P-values < 0.05.
LncRNA-mRNA network analysis
Pearson’s correlation was used to determine whether the expression levels of DE lncRNAs were correlated with DE mRNAs between G20W and G55W. The lncRNA-mRNA correlations with P-values <0.05 were retained. We used Cytoscape 3.4.0 to visualize lncRNA-mRNA networks.
Quantitative RT-PCR and data analysis
To verify the accuracy of the data obtained by high-throughput sequencing, nine lncRNAs were randomly selected, and the results were confirmed by qRT-PCR for nine randomly selected lncRNAs. Total RNA was isolated from six Gushi chicken breast muscle samples at 20weeks and 55weeks of age, which was the same as for the RNA-seq samples. All gene-specific qRT-PCR primers were designed using Primer 6.0 software and synthesized by Shenggong Biological Engineering (Wuhan) Limited by Share Ltd. All qPCR was conducted using the LightCycler 96 Real-Time PCR system (Roche Applied Science). In a 10 μL reaction mixture, 1.0 μL of cDNA was used for amplification with 5 μL of SYBR Premix Ex Taq (TRNaseH Plus) (TaKaRa, Dalian, China), 0.5 μL of forward primer (10 μM), 0.5 μL of universal primer (10 μM), and 3 μL of deionized water. The reactions were incubated at 95°C for 1 min, followed by 35 cycles of 95°C for 15 s, 60°C for 45 s, and 72°C for 40 s and an extension for 10 min at 72°C. All reactions were performed in triplicate. The threshold cycle (CT) was determined using the default threshold settings, and the data were analyzed using the 2–ΔΔCt program [32]. Statistical analyses of the data were performed using SPSS 22.0 (SPSS Inc., Chicago, IL, USA). The data are presented as the means ± standard error (SE).
Results
Overview of lncRNA expression
To identify DE lncRNAs, we established six cDNA libraries that represented two different physiological stages: G20W-1, G20W-2, and G20W-3 from the breast muscle tissue of 20-week-old juvenile hens and G55W-1, G55W-2, and G55W-3 from the breast muscle tissue of 55-week-old laying hens. RNA sequencing obtained a total of 43.51 Gb of data, each stage averaging 7.3 Gb of data. More than 61.05% of the clean reads were perfectly mapped to the chicken reference genome, and 59.69 to 63.33% uniquely mapped reads were obtained from the total reads from the six samples. The Q30 results in each sample were > 90%, and the GC percentage was less than 60%, as listed in Table 1. After mapping the reads to each chromosome, no significant difference between the different physiological stages was observed, and the number of reads mapping to different chromosomes was consistent with the chromosome length (Fig 1). The same results were obtained in a previous study [33].
The outer ring represents the chicken genome labeled with chromosome number and position. The green circle shows the distribution of identified lncRNAs in G20W, and the orange circle shows the distribution in G55W.
After assembly, we obtained 4, 404 lncRNAs from the two groups (Fig 2A). Among the total lncRNAs, 4, 029 and 3, 452 lncRNAs were identified in the G20W and G55W samples, respectively, with 3, 074 lncRNAs observed to be expressed in both groups. We found that 952 lncRNAs were specifically expressed in G20W and that 378 lncRNAs were expressed only in G55W (Fig 2B). Furthermore, we examined the lengths and chromosomal locations of the lncRNAs. The length of the lncRNAs ranged from 200 to 10, 000 nucleotides (nt). 60% of the total number of lncRNAs were approximately 200 to 2, 000 nt in length (Fig 2C). Interestingly, the lncRNAs in chicken breast tissue were similar in length to the lncRNAs in mouse and Atlantic salmon [34, 35]. To understand the distribution of the chicken lncRNA reference genome, we analyzed the distribution of lncRNAs in Gushi chicken breast samples. The results showed that lncRNAs were distributed across the chicken chromosomes (chromosomes 1–28, W), with chromosomes 1–4 and Z containing many lncRNAs (Fig 2D). that is, the larger chromosomes contained more lncRNAs. Widespread transcription along all chromosomes with some bias in transcriptional activity on specific chromosomes indicates that the presence of lncRNAs is not due to transcriptional noise.
(A) Computational pipeline for the systematic identification of lncRNAs. (B) Venn diagram illustrating the overlap of lncRNAs in breast muscle tissue libraries from G20W and G55W chickens. G20W: 20 weeks of age; G55W: 55 weeks of age. Length (C) and chromosomal (D) distribution of Gushi chicken lncRNAs.
Differentially expressed lncRNAs and mRNAs in chicken breast muscle
The correlation coefficient can represent the degree of similarity between samples. We found that G20W and G55W were highly correlated (Fig 3A). We used Heatmap3 [36] to perform an unsupervised cluster analysis of the differences in lncRNA expression in breast muscle tissue at different developmental stages. The results of the cluster analysis showed that distinct lncRNA expression patterns are associated with G20W and G55W (Fig 3B). We identified 234 DE lncRNAs between the G20W and G55W breast muscle tissues (|Log2 Fold-Change| ≥ 1.0, P < 0.05). Among the total 234 lncRNAs, 131 were upregulated in G20W compared with G55W, and 55 were downregulated (Fig 3C). Among the top 10 expressed lncRNAs in the libraries, the first 7 lncRNAs were present in both samples. None of the top 10 lncRNAs were specifically expressed in G20W or G55W. The fold change in the expression of the DE lncRNAs varied from -1.8582 to 2.3798 (Table 2). These findings may help to identify DE lncRNAs at specific stages in chickens.
(A) Representative of the degree of similarity between samples. The correlation coefficient is represented by color; deeper color represents a stronger correlation. (B) Clustering analysis of differentially expressed lncRNAs in G20W and G55W. (C) Volcano plot of differentially expressed lncRNAs in G20W and G55W. Green and red represent downregulated and upregulated expression, respectively.
To further compare transcriptional profiles and networks in breast tissue, gene expression profiling was performed using the Illumina HiSeq 2000. Finally, we identified a total of 881 DE genes in breast tissue samples, comprising 451 upregulated and 430 downregulated genes in G20W compared with G55W (|Log2 fold-change| ≥ 1.0, P < 0.05).
GO enrichment analysis
We predicted the cis and trans targets of DE lncRNAs to further investigate how lncRNAs interact with the mRNA of their target genes to regulate chicken breast meat quality and to identify key molecular participants in this process. LncRNA can not only regulate the expression of neighboring protein-coding genes through a cis mechanism [37, 38], but also regulate the expression of genes located on other chromosomes via a trans mechanism [21, 39]. In this study, bioinformatic analysis was used to predict the potential target genes of DE lncRNAs via cis and trans. A total of 315 cis-mRNA pairs (S1 Table) and 15, 722 random trans-mRNA pairs (S2 Table) were identified in the present study. Furthermore, we identified the genes common to both categories and were therefore able to filter key lncRNAs and their associated target genes that could affect the quality of breast muscle tissue.
We used GOseq to compare the GO classifications of the groups of upregulated and downregulated genes (adjusted P < 0.05) [40]. Overall, 203 DE lncRNAs and coexpressed mRNAs were enriched for 955 GO terms in the cis targets of lncRNAs. Among them, 689 of the enriched GO terms were categorized as biological processes (BP), 153 as for molecular functions (MF), and 113 ascellular components (CC). The top 30 significantly enriched GO terms are shown in Fig 4A. The most significantly enriched GO term between the G20W and G55W groups was cAMP-dependent protein kinase inhibitor activity. For trans targets, 774 DE lncRNAs and coexpressed mRNAs were enriched for 6, 398 GO terms. Among them, 4, 852 GO terms were enriched for BP, 961 were enriched for MF, and 585 were enriched for CC. The top 30 significantly enriched GO terms are shown in Fig 4B. The top three GO terms were membrane-bound protein binding, protein metabolic process, and cellular component organization or biogenesis.
(A) Categories of biological processes, cellular components and molecular functions of the target genes of differentially expressed lncRNAs and mRNAs (cis). (B) Categories of the biological processes, cellular components and molecular functions of the target genes of differentially expressed lncRNAs and mRNAs (trans).
KEGG pathway enrichment analysis
We enriched the differentially coexpressed lncRNAs and mRNAs in the G20W and G55W groups through KEGG pathway analysis to further identify the biological pathways in which these genes participate. The pathways of the DE lncRNA target genes were significantly enriched for tight junctions, ABC transporters, retinol metabolism, fatty acid degradation and MAPK signaling pathways in the cis targets of lncRNAs and mRNAs (Fig 5A). For trans targets of lncRNAs and mRNAs, the proteasome, ECM-receptor interaction, focal adhesion, adherens junctions, ubiquitin-mediated proteolysis, endocytosis, the Wnt-signaling pathway and the calcium-signaling pathway were identified (Fig 5B). We observed that several genes were involved in skeletal muscle development and lipid metabolism, such as COL6A3, COL4A1, COL5A1, COL4A2, LAMA4, LAMB4, LAMC1, Lpin1, PGK1, UBE2B, UBE2Q2, and UBE2D1 (S3 Table). Some of these pathways may be associated with muscle development and meat quality.
(A) Scatter plot of the top 20 pathways enriched for differentially expressed lncRNAs and mRNAs in breast tissue from G20W and G55W chickens (cis). The abscissa represents the richness factor, and the ordinate represents the enriched pathway terms. The Q-value represents the corrected P (B) Scatter plot of the top 20 pathways enriched for differentially expressed lncRNAs and mRNAs in breast tissue from G20W and G55W chickens (trans). The abscissa represents the richness factor, and the ordinate represents the enriched pathway terms. The Q-value represents the corrected P.
LncRNA-mRNA network analysis
An lncRNA-mRNA network was constructed based on the correlation analysis between the differentially coexpressed lncRNAs and mRNAs, which were determined to be enriched in several pathways. We constructed gene coexpression networks to identify interactions among DE mRNAs and lncRNAs (Fig 6). The majority of DE lncRNAs and target genes (G55W vs G20W) were identified in glycerophospholipid metabolism, ECM-receptor interactions, and amino acid pathway biosynthesis. Interestingly, the expression of PGK1, AGPAT9, and PHOSPHO1 was upregulated, and their corresponding lncRNAs, ALDBGALG000000005749, ALDBGALG0000005602, and ALDBGALG0000004644, respectively, were downregulated. COL2A1, COL6A3, COL5A1, LAMB1, LAMB4, and LAMA4 were downregulated, and their corresponding lncRNAs were upregulated.
Circles represent coding genes, and arrows represent lncRNAs. Red nodes represent upregulated nodes, and green nodes represent downregulated nodes. Yellow nodes represent pathway terms. The node size represents -logP; the smaller the P value is, the greater is the node size. Lines between lncRNA-mRNA represent interactions between them.
Verification of lncRNA expression profiles using qRT-PCR
To verify the accuracy of our RNA-seq results, we randomly selected 9 candidate lncRNAs from common DEGs in G20W and G55W. Total RNA was obtained by the same method as for RNA-seq, and β-actin was used as a reference gene. The primer sequences are presented in S4 Table. The results showed that the expression of ALDBGALG (5114, 1940, 5729, 1736 and 4424) was higher at 20 weeks than at 55 weeks, whereas the expression of ALDBGALG (2353, 0605, 0119 and 1104) was much higher at 55 weeks than at 20 weeks. All the lncRNAs expression levels were consistent with the RNA-Seq findings (Fig 7), suggesting that the detection and expression abundance of lncRNAs in our RNA-seq was highly accurate.
Data were analyzed by the 2-ΔΔCt method using β-actin as a reference gene. Each column represents the mean ±SE. Different letters indicate significant differences in expression levels between the two stages (P ≤ 0.05). Black bars represent read from RNA-Seq. Gray bars represent the results of qRT-PCR.
Discussion
Meat quality traits are comprehensive economic characteristics controlled by polygenic traits, which are easily affected by many factors, such as castration, gender, age, and species. Species, as a major genetic factor, plays the leading role in defining the characteristics of meat quality. IMF is one of the important factors influencing meat quality and is mainly distributed in muscle and muscle fibers. IMF contains many phospholipids, and phospholipids containing large amounts of polyunsaturated fatty acids can be produced via lipid degradation reactions. IMF deposition can loosen muscle fibers, fat and connective tissue in the overlapping structure, thereby decreasing physical strength and promoting the separation of muscle fiber bundles, which especially improves the tenderness of meat [41, 42].
LncRNAs are widely distributed among mammals. Their expression is spatiotemporally and tissue specific and plays important roles in many aspects, including epigenetic regulation, transcriptional regulation and posttranscriptional regulation [43]. With the rapid development of next generation sequencing technologies, lncRNAs have been examined as novel regulatory players in cellular and biological processes using transcriptome analyses [44–49]. Some previous studies have investigated lncRNA in the muscle tissue of chickens [50, 51]. Even so, the diversity of lncRNA expression and the biological lncRNA functions at two different chicken physiological stages remain unclear. To the best of our knowledge, the lncRNA profiles in breast tissue from juvenile and laying-stage chickens are not well known. In this study, we comprehensively investigated lncRNAs in breast tissue from Gushi chickens at two physiological stages, juvenile and laying.
We used RNA-Seq technology to sequence the Gushi chicken breast tissue and identified approximately 4, 404 putative noncoding lncRNAs (4, 026 in G20W and 3, 452 in G55W). Among these lncRNAs, 131 were upregulated and 55 were downregulated in G20W compared with G55W. One difficulty in lncRNA study is that, in contrast to mRNA functions, the functions of most lncRNAs have not been determined, and there is no existing database that can be used to identify their functional annotations. Increasing evidence suggests that lncRNAs can regulate the expression of neighboring mRNAs and that their expression is highly correlated with that expression of neighboring mRNAs [9, 52]. Thus, we searched for trans target genes of lncRNAs to help predict the function of the lncRNAs.
The regulation of muscle development and composition is likely a function of complex networks. Therefore, examining the regulatory network is the preferred method of analysis [53]. In the present study, we constructed a coexpression network of lncRNAs and coding gene transcripts to predict the potential biological functions of DE lncRNAs. DE lncRNAs and mRNAs were enriched in several pathways in common, including protein processing in endoplasmic reticulum, ECM-receptor interaction, glycerophospholipid metabolism, ubiquitin-mediated proteolysis and the biosynthesis of amino acids, which play crucial roles in muscle development, glycerophospholipid metabolism, protein catabolism, energy metabolism and lipid metabolism processes.
The deposition of IMF is affected by a variety of genetic factors, and changes in the IMF content are closely related to different breeds and developmental stages [54]. A previous study demonstrated that ECM-receptor interaction might form a network with pathways related to lipid metabolism to influence the deposition of IMF [55]. The extracellular matrix (ECM) is a part of the connective tissue layers surrounding muscle fibers [56] and is a dynamic network structure composed of fibrous and nonfibrous proteins including collagen, proteoglycans and glycoproteins and other large molecules. The ECM provides support for the survival and activity of cells, and its effects include signal transduction that affects cell shape, function, metabolism, migration, proliferation and differentiation. Overall, there might exist a relationship between the expression of genes related to encoding collagen and ECM and meat quality traits [57]. Our results suggest that some DE genes encode collagen and epidermal growth factor-like domains, including collagen type VI alpha 3 chain (COL6A3), collagen type IV alpha 1 chain (COL4A1), collagen type V alpha 1 chain (COL5A1), collagen type IV alpha 2 chain (COL4A2), collagen type II alpha 1 chain (COL2A1), laminin subunit alpha 4 (LAMA4), laminin subunit beta 4 (LAMB4), laminin subunit gamma 1 (LAMC1), laminin subunit beta 1 (LAMB1), and reelin (RELN). It was suggested that these genes have an important relationship to the ECM [58]. Among them, LAMA4, which is encoded by the LAMA4 gene, mainly functions outside neuronal membranes, in muscle fibers during the development of endothelial cells and surrounding the basement membrane. Another study suggested that LAMA4 likely functions as a key protein in muscle development during the embryonic period. In our study, COL2A1, COL6A3, COL5A1, LAMB1, LAMB4, and LAMA4 were downregulated, and their corresponding DE lncRNAs (ALDBGALG 2611, 5463, 5689, 3544, 4424 and 5595) were upregulated. The upregulation of a given lncRNA may inhibit its target gene and affect meat quality, a possibility that requires further study.
Lipid phosphate phosphohydrolase 1 (Lpin1) is a member of the LPIN family and is a key gene in adipocyte differentiation and lipid metabolism. In a previous study using RNA-Seq in broilers of two groups with different shear forces of breast tissue, the LPIN1 gene was found to be downregulated and involved in the regulation of meat tenderness [59]. The corresponding lncRNA, ALDBGALG599, may have a similar function, and both may influence meat quality. Interestingly, we found that PGK1 (phosphoglycerate kinase 1), a key glycolytic enzyme involved in the biosynthesis of amino acids, was upregulated at two different chicken physiological stages. A lack of this enzyme in organisms can cause metabolic disorders. Another study suggested that PGK1 was involved in acid synthesis and influenced pig IMF deposition in skeletal muscle, affecting the flavor of meat. In our study, the corresponding lncRNA of PGK1, ALDBGALG5749 was downregulated, and its function in chickens should be examined in future studies.
Ubiquitin-mediated proteolysis may play a crucial role in protein degradation during the aging process. Protein degradation also affects muscle shear and thus alters meat tenderness [2, 60, 61]. Our results suggest that some differentially coexpressed genes are associated with ubiquitin modification, including ubiquitin-conjugating enzyme E2 B (UBE2B), ubiquitin conjugating enzyme E2 Q2 (UBE2Q2), ubiquitin conjugating enzyme E2 D1 (UBE2D1), ubiquitin conjugating enzyme E2 D3(UBE2D3) and ubiquitin conjugating enzyme E2 K (UBE2K). The results indicate that co-expressed genes might affect the expression of these ubiquitin-related genes and thus participate in protein degradation in breast muscle.
Conclusion
In conclusion, we first generated the expression profiles of lncRNAs in breast tissue from two developmental stages of Gushi chicken (G20W and G55W) based on RNA-Seq. Our results suggest that those lncRNAs might play important roles in meat quality. Overall, this study takes the first step toward understanding the molecular mechanisms underlying variations in poultry meat quality and helps to provide a theoretical basis for further research. However, the regulatory mechanisms of lncRNAs require further verification.
Supporting information
S3 Table. KEGG pathways associated with lncRNA-mRNA interaction network analysis of differentially expressed genes.
https://doi.org/10.1371/journal.pone.0215006.s003
(XLSX)
S4 Table. Primers used for qRT-PCR in this study.
https://doi.org/10.1371/journal.pone.0215006.s004
(XLSX)
References
- 1. Meluzzi A, Sirri F, Castellini C, Roncarati A, Melotti P, Franchini A, et al. Influence of genotype and feeding on chemical composition of organic chicken meat. Italian Journal of Animal Science. 2010; 8(2s): 766–768.
- 2. Zhang M, Yan FB, Li F, Jiang KR, Li DH, Han RL, et al. Genome-wide DNA methylation profiles reveal novel candidate genes associated with meat quality at different age stages in hens. Scientific Reports. 2017; 7: 45564. pmid:28378745
- 3. Zhang CY, Wang Z, Bruce HL, Janz J, Goddard E, Moore S, et al. Associations between single nucleotide polymorphisms in 33 candidate genes and meat quality traits in commercial pigs. Animal Genetics. 2014; 45(4): 508–516. pmid:24707962
- 4. Lopes LS, Martins SR, Chizzotti ML, Busato KC, Oliveira IM, Machado Neto OR, et al. Meat quality and fatty acid profile of Brazilian goats subjected to different nutritional treatments. Meat Science. 2014; 97(4): 602–608. pmid:24795167
- 5. Gregory NG, Sant'Ana ADS. How climatic changes could affect meat quality. Food Research International. 2010; 43(7): 1866–1873.
- 6. Fanatico AC, Pillai PB, Emmert JL, Owens CM. Meat quality of slow- and fast-growing chicken genotypes fed low-nutrient or standard diets and raised indoors or with outdoor access. Poult Sci. 2007; 86(10): 2245–2255. pmid:17878457
- 7. Hou X, Yang Y, Zhu S, Hua C, Zhou R, Mu Y, et al. Comparison of skeletal muscle miRNA and mRNA profiles among three pig breeds. Molecular Genetics & Genomics. 2016; 291(2): 559–573.
- 8. Tellam RL, Cockett NE, Vuocolo T, Bidwell CA. Genes Contributing to Genetic Variation of Muscling in Sheep. Frontiers in Genetics. 2012; 3: 164. pmid:22952470
- 9. Yang L, Yi K, Wang H, Zhao Y, Xi M. Comprehensive analysis of lncRNAs microarray profile and mRNA-lncRNA co-expression in oncogenic HPV-positive cervical cancer cell lines. Oncotarget. 2016; 7(31): 49917–49929. pmid:27363024
- 10. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009; 458(7235): 223–227. pmid:19182780.
- 11. Kedra D, Bussotti G, Prieto P, Sørensen P, Bagnato A, Mckay SD, et al., editors. Identification and conservation of novel long noncoding RNAs in cattle using RNASeq data. Meeting of the European Federation of Animal Science; 2013.
- 12. Li XQ, Ren ZX, Li K, Huang JJ, Huang ZT, Zhou TR, et al. Key Anti-Fibrosis Associated Long Noncoding RNAs Identified in Human Hepatic Stellate Cell via Transcriptome Sequencing Analysis. International Journal of Molecular Sciences. 2018; 19(3): 675.
- 13. Novikova IV, Hennelly SP, Tung CS, Sanbonmatsu KY. Rise of the RNA machines: exploring the structure of long non-coding RNAs. Journal of Molecular Biology. 2013; 425(19): 3731–3746. pmid:23467124
- 14. Nelson BR, Makarewich CA, Anderson DM, Winders BR, Troupes CD, Wu F, et al. A peptide encoded by a transcript annotated as long noncoding RNA enhances SERCA activity in muscle. Science. 2016; 351(6270): 271–275. pmid:26816378
- 15. Anderson D, Anderson K, Chang CL, Makarewich C, Nelson B, Mcanally J, et al. A Micropeptide Encoded by a Putative Long Noncoding RNA Regulates Muscle Performance. Cell. 2015; 160(4): 595. pmid:25640239
- 16. Anderson DM, Makarewich CA, Anderson KM, Shelton JM, Bezprozvannaya S, Basselduby R, et al. Widespread control of calcium signaling by a family of SERCA-inhibiting micropeptides. Science Signaling. 2016; 9(457): ra119. pmid:27923914
- 17. Klattenhoff CA, Scheuermann JC, Surface LE, Bradley RK, Fields PA, Steinhauser ML, et al. Braveheart, a long noncoding RNA required for cardiovascular lineage commitment. Cell. 2013; 152(3): 570–583. pmid:23352431
- 18. Shen H, Mcelhinny AS, Cao Y, Gao P, Liu J, Bronson R, et al. The Notch coactivator, MAML1, functions as a novel coactivator for MEF2C-mediated transcription and is required for normal myogenesis. Genes Dev. 2006; 20(6): 675–688. pmid:16510869
- 19. Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, et al. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell. 2011; 147(2): 358–369. pmid:22000014
- 20. Zhu M, Liu J, Xiao J, Yang L, Cai M, Shen H, et al. Lnc-mg is a long non-coding RNA that promotes myogenesis. Nature Communications. 2017; 8: 14718. pmid:28281528
- 21. Cai B, Li Z, Ma M, Wang Z, Han P, Abdalla BA, et al. LncRNA-Six1 Encodes a Micropeptide to Activate Six1 in Cis and Is Involved in Cell Proliferation and Muscle Growth. Frontiers in Physiology. 2017; 8.
- 22. Li H, Ma Z, Jia L, Li Y, Xu C, Wang T, et al. Systematic analysis of the regulatory functions of microRNAs in chicken hepatic lipid metabolism. Sci Rep. 2016; 6: 31766. pmid:27535581
- 23. Bai B, Wu J, Sheng WT, Zhou B, Zhou LJ, Zhuang W, et al. Comparative Analysis of Anther Transcriptome Profiles of Two Different Rice Male Sterile Lines Genotypes under Cold Stress. International Journal of Molecular Sciences. 2015; 16(5): 11398. pmid:25993302
- 24. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009; 25(9): 1105–1111. pmid:19289445
- 25. Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Adiconis X, et al. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs. Nature Biotechnology. 2010; 28(5): 503. pmid:20436462
- 26. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology. 2010; 28(5): 511–515. pmid:20436464
- 27. Li A, Zhang J, Zhou Z, Wang L, Liu Y, Liu Y. ALDB: a domestic-animal long noncoding RNA database. PloS one. 2015; 10(4): e0124003. pmid:25853886
- 28. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research. 2007; 35(Web Server issue): W345. pmid:17631615
- 29. Wang L, Park HJ, Dasari S, Wang S, Kocher JP, Li W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Research. 2013; 41(6): e74–e74. pmid:23335781
- 30. Kaplan RS, Roll R. Investor Evaluation of Accounting Information: Some Empirical Evidence. Journal of Business. 1972; 45(2): 225–257.
- 31. Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Research. 2016; 44(Database issue): D279–D285. pmid:26673716
- 32. Costa V, Angelini C, De FI, Ciccodicola A. Uncovering the Complexity of Transcriptomes with RNA-Seq. Journal of Biomedicine & Biotechnology. 2010; 2010(5757): 853916.
- 33. Zhang Y, Li D, Han R, Wang Y, Li G, Liu X, et al. Transcriptome analysis of the pectoral muscles of local chickens and commercial broilers using Ribo-Zero ribonucleic acid sequencing. Plos One. 2017; 12(9): e0184115. pmid:28863190
- 34. Okazaki Y, Furuno M, Kasukawa T, Adachi J, Bono H, Kondo S, et al. Okazaki Y. Analysis of the mouse transcriptome based on functional annotation of 60,770 full-length cDNAs. Nature 420, 563–573. Nature. 2003; 420(6915): 563–573. pmid:12466851
- 35. Boltaña S, Valenzuelamiranda D, Aguilar A, Mackenzie S, Gallardoescárate C. Long noncoding RNAs (lncRNAs) dynamics evidence immunomodulation during ISAV-Infected Atlantic salmon (Salmo salar). Scientific Reports. 2016; 6: 22698. pmid:26939752
- 36. Zhao S, Guo Y, Sheng Q, Shyr Y. Advanced Heat Map and Clustering Analysis Using Heatmap3. BioMed Research International. 2014; 2014: 1–6. https://doi.org/10.1155/2014/986048.
- 37. Bu Q, Hu Z, Chen F, Zhu R, Deng Y, Shao X, et al. Transcriptome analysis of long non‐coding RNAs of the nucleus accumbens in cocaine‐conditioned mice. Journal of Neurochemistry. 2012; 123(5): 790–799. pmid:22957495
- 38. Han L, Zhang K, Shi Z, Zhang J, Zhu J, Zhu S, et al. LncRNA profile of glioblastoma reveals the potential role of lncRNAs in contributing to glioblastoma pathogenesis. International Journal of Oncology. 2012; 40(6): 2004. pmid:22446686
- 39. Yang L, Yi K, Wang H, Zhao Y, Xi M. Comprehensive analysis of lncRNAs microarray profile and mRNA–lncRNA co-expression in oncogenic HPV-positive cervical cancer cell lines. Oncotarget. 2016; 7(31): 49917–49929. pmid:27363024
- 40. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology. 2010; 11(2): R14. pmid:20132535
- 41. Cannata S, Engle TE, Moeller SJ, Zerby HN, Radunz AE, Green MD, et al. Effect of visual marbling on sensory properties and quality traits of pork loin. Meat Sci. 2010; 85(3): 428–434. pmid:20416803.
- 42. Okeudo NJ, Moss BW. Interrelationships amongst carcass and meat quality characteristics of sheep. Meat Science. 2005; 69(69): 1–8.
- 43. Quan M, Chen J, Zhang D. Exploring the Secrets of Long Noncoding RNAs. International Journal of Molecular Sciences. 2015; 16(3): 5467–5496. pmid:25764159
- 44. Fu LL, Li CJ, Xu Y, Li LY, Zhou X, Li DD, et al. Role of lncRNAs as Novel Biomarkers and Therapeutic Targets in Ovarian Cancer. Crit Rev Eukaryot Gene Expr. 2017; 27(2): 183–195. pmid:28845767
- 45. Yuwaiman C, Owen N, Lees J, Tagalakis A, Hart S, Webster A, et al. Genome-wide RNA-Sequencing analysis identifies a distinct fibrosis gene signature in the conjunctiva after glaucoma surgery. Scientific Reports. 2017; 7(1).
- 46. Guo CJ, Xiao X, Sheng L, Chen L, Zhong W, Li H, et al. RNA Sequencing and Bioinformatics Analysis Implicate the Regulatory Role of a Long Noncoding RNA-mRNA Network in Hepatic Stellate Cell Activation. Cellular Physiology & Biochemistry. 2017; 42(5): 2030–2042.
- 47. Zhang K, Han X, Zhang Z, Zheng L, Hu Z, Yao Q, et al. The liver-enriched lnc-LFAR1 promotes liver fibrosis by activating TGFβ and Notch pathways. Nature Communications. 2017; 8(1): 144. pmid:28747678
- 48. Sun Y, Zeng C, Gan S, Li H, Cheng Y, Chen D, et al. LncRNA HOTTIP-Mediated HOXA11 Expression Promotes Cell Growth, Migration and Inhibits Cell Apoptosis in Breast Cancer. International Journal of Molecular Sciences. 2018; 19(2): 472.
- 49. Li Y, Chen X, Sun H, Wang H. Long non-coding RNAs in the regulation of skeletal myogenesis and muscle diseases. Cancer Letters. 2017; 417: 58–64. pmid:29253523
- 50. Cabianca DS, Casa V, Bodega B, Xynos A, Ginelli E, Tanaka Y, et al. A long ncRNA links copy number variation to a polycomb/trithorax epigenetic switch in FSHD muscular dystrophy. Cell. 2012; 149(4): 819–831. pmid:22541069.
- 51. Sun L, Zhang L, Liu H. Prediction of Long Non-Coding RNAs Based on RNA-Seq*. Progress in Biochemistry and Biophysics. 2013; 39(12): 1156–1166. https://doi.org/10.3724/sp.j.1206.2012.00287.
- 52. Wang Y, Xue S, Liu X, Liu H, Hu T, Qiu X, et al. Analyses of Long Non-Coding RNA and mRNA profiling using RNA sequencing during the pre-implantation phases in pig endometrium. Scientific Reports. 2016; 6: 20238. pmid:26822553
- 53. Du YF, Ding QL, Li YM, Fang WR. Identification of Differentially Expressed Genes and Pathways for Myofiber Characteristics in Soleus Muscles between Chicken Breeds Differing in Meat Quality. Anim Biotechnol. 2016: 1–11. pmid:27623936.
- 54. Gui-Juan QU, Yang LY, Qin GX, Dong XQ. Developmental changes of LPL gene expression of longissimusdorsi and their effects on IMF in different crossbred cattle breeds. Chinese Journal of Veterinary Science. 2011.
- 55. Cui HX, Liu RR, Zhao GP, Zheng MQ, Chen JL, Wen J. Identification of differentially expressed genes and pathways for intramuscular fat deposition in pectoralis major tissues of fast-and slow-growing chickens. BMC Genomics. 2012; 13: 213. pmid:22646994.
- 56. Mariman ECM, Ping W. Adipocyte extracellular matrix composition, dynamics and role in obesity. Cellular & Molecular Life Sciences Cmls. 2010; 67(8): 1277–1292.
- 57. Smith SH, Judge MD. Relationship between pyridinoline concentration and thermal stability of bovine intramuscular collagen. Journal of Animal Science. 1991; 69(5): 1989–1993. pmid:2066308
- 58. Nakajima I, Muroya S, Tanabe R, Chikuni K. Extracellular matrix development during differentiation into adipocytes with a unique increase in type V and VI collagen. Biology of the Cell. 2002; 94(3): 197–203. pmid:12206658
- 59. Piorkowska K, Zukowski K, Nowak J, Poltowicz K, Ropka-Molik K, Gurgul A. Genome-wide RNA-Seq analysis of breast muscles of two broiler chicken groups differing in shear force. Anim Genet. 2016; 47(1): 68–80. pmid:26592359.
- 60. Jin L, Jiang Z, Xia Y, Lou P, Chen L, Wang H, et al. Genome-wide DNA methylation changes in skeletal muscle between young and middle-aged pigs. BMC Genomics. 2014; 15(1): 1–12.
- 61. Zahn JM, Sonu R, Vogel H, Crane E, Mazanmamczarz K, Rabkin R, et al. Transcriptional Profiling of Aging in Human Muscle Reveals a Common Aging Signature. Plos Genetics. 2006; 2(7): e115. pmid:16789832