- 1Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province, College of Animal Science and Technology, Sichuan Agricultural University, Chengdu, China
- 2College of Life Science, Sichuan Agricultural University, Ya’an, China
The liver is the major organ of lipid biosynthesis in the chicken. In laying hens, the liver synthesizes most of the yolk precursors and transports them to developing follicles to produce eggs. However, a systematic investigation of the long non-coding RNA (lncRNA) and mRNA transcriptome in liver across developmental stages is needed. Here, we constructed 12 RNA libraries from liver tissue during four developmental stages: juvenile (day 60), sexual maturity (day 133), peak laying (day 220), and broodiness (day 400). A total of 16,930 putative lncRNAs and 18,260 mRNAs were identified. More than half (53.70%) of the lncRNAs were intergenic lncRNAs. The temporal expression pattern showed that lncRNAs were more restricted than mRNAs. We identified numerous differentially expressed lncRNAs and mRNAs by pairwise comparison between the four developmental stages and found that VTG2, RBP, and a novel protein-coding gene were differentially expressed in all stages. Time-series analysis showed that the modules with upregulated genes were involved in lipid metabolism processes. Co-expression networks suggested functional relatedness between mRNAs and lncRNAs; the DE-lncRNAs were mainly involved in lipid biosynthesis and metabolism processes. We showed that the liver transcriptome varies across different developmental stages. Our results improve our understanding of the molecular mechanisms underlying liver development in chickens.
Introduction
The liver is the largest metabolic organ, and it is involved in the metabolism of numerous biologically important compounds, synthesis of essential proteins, detoxification of exogenous materials, and host defense against invading pathogens (Grant, 1991; Brockmöller and Ivar, 1994; Calne, 2000; Rui, 2014). The liver produces a large number of fatty acids and play a crucial role in lipid transportation through lipoprotein synthesis, depending on species (Nguyen et al., 2008; Mashek, 2013). In chickens, over 70% of the de novo lipogenesis occurs in the liver (Leveille et al., 1968; O’Hea and Leveille, 1969), especially during the egg production period. The production of egg yolk proteins in chicken liver is mainly regulated by several steroid hormones, such as estrogen. As hens mature sexually, the secretion of hydrophobic lipids including free fatty acids, triacylglycerols, and cholesteryl esters begins in the liver and these are assembled together to form two major components of yolk protein precursors including very low density apolipoprotein II (ApoVLDL II) and vitellogenin II (VTG II) under estrogen stimulation (Wiskocil et al., 1980; Denslow et al., 1999; Walzem et al., 1999). These particles are subsequently transported to blood and circulated to the developing oocyte for embryo growth and development (Nimpf and Schneider, 1991; Schneider, 1995).
Large-scale transcriptome analyses reveal that about 75% of the human genome can be transcribed into RNAs, and most of these transcripts are defined as different types of non-coding RNAs (ncRNAs) (Djebali et al., 2012; Liu et al., 2013). Of these, long non-coding RNAs (lncRNAs) are the most abundant and play versatile role in gene regulation of expression (Derrien et al., 2012a). The lncRNAs are generally >200 nucleotides (nt) in length with little or no coding potential (Esteller, 2011). LncRNAs were initially considered to be “noise” of genome transcription and to lack biological function. However, studies in different species have fully demonstrated that lncRNAs are an important regulatory factor in a series of biological processes, such as controlled recruitment of transcription factors (Martianov et al., 2007; Maamar et al., 2013; Ng et al., 2013), regulation of mRNA decay and translation rate (Gong and Maquat, 2011; Carrieri et al., 2012), and formation and function of subnuclear structures (Mao et al., 2010; Yang et al., 2011).
The chicken is an established model organism for studying vertebrate development (Stern, 2005). Compared with mammals, the liver is the accessory organ of lipid biosynthesis in the chicken (Leveille et al., 1975). Thus, liver is responsible for lipid metabolism and homeostasis that play a critical effect on chicken growth and laying performance. It is also important for birds to adapt environmental changes through the regulation of the lipid metabolism in liver (Bedu et al., 2002). Although the genes and gene products related to hepatic lipid metabolism in chicken liver are very similar to those in mammals, they differed in specialized functions (Kirchgessner et al., 1987). For example, the microsomal triglyceride transfer protein (MTTP) is known to assist lipoprotein assembly to form very low density lipoprotein (VLDL) in mammals (Hussain et al., 2008, 2012); however, a previous study indicated that upregulation of MTTP in the liver does not require increased VLDL assembly during laying process in chicken (Ivessa et al., 2013). In addition, a battery of genes associated to lipid metabolism were lost from the chicken genome during the evolutionary process (Ðaković et al., 2014). Most RNA-seq studies in chickens have focused on protein-coding genes and ignored lncRNAs because of their poor conservation and low expression abundance in tissues compared with mRNAs. Therefore, the repertoire of potential genes and lncRNAs involved in hepatic lipid metabolism remains to be completely elucidated in laying chickens.
There, in this study, we comprehensively investigated the lncRNAs and mRNAs in chicken liver at four developmental stages using high-throughput RNA sequencing technology. We identified many differentially expressed transcripts and explored these enrichment pathways across the four stages. In addition, by constructing a co-expression network of lncRNA and mRNAs, we predicted the potential functions of lncRNAs. These results would contribute to better comprehension of transcriptome changes and comparative understanding of molecular mechanisms of liver development in chickens.
Materials and Methods
Animals and Liver Tissue Collection
In this study, the healthy hen chickens of Luhua (a Chinese indigenous breed) were obtained from the Experimental Chicken Farm of Sichuan Agricultural University. Chickens were examined at the four developmental stages: juvenile period (JP, 60 days of age), sexual maturity period (SM, 133 days of age), peak laying period (PL, 220 days of age), and broodiness period (BP, 400 days of age). Birds were euthanized by intravenous injection containing 2% pentobarbital sodium (25 mg/kg of BW). Liver tissue was rapidly collected from three biological replicates at each stage and promptly frozen in liquid nitrogen. All samples were stored at −80°C until total RNA extraction. All experimental methods and sample collection of this study were conducted in accordance with guidelines provided by the Institutional Animal Care and Use Committee (IACUC) of Sichuan Agricultural University, under permit No. DKY-B20171910.
RNA Isolation, Library Preparation, and Sequencing
Total RNA of each sample was extracted using RNAiso Plus reagent (TaKaRa, Otsu, Shiga, Japan) following the manufacturer’s instruction. We estimated the integrity and quality of total RNA using Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, United States) with an RNA 6000 Nano kit. Subsequently, 12 strand-specific libraries were generated after depleting rRNA by the Ribo-ZeroTM Gold Kit (Illumina, San Diego, CA, United States) and then sequenced on the Illumina HiSeq X Ten platform (Illumina Inc., San Diego, CA, United States) with a paired-end sequencing length of 150 bp (PE150) at Novogene Corporation (Beijing, China).
Identification of lncRNAs
We removed the low-quality reads, adaptor sequences, empty reads, and ribosomal (r)RNA reads from the raw data in order to obtain high-quality lncRNAs. The processed reads were mapped to the chicken reference genome (Gallus gallus-5.0) using HISAT2 2.1.0 (Kim et al., 2015), and the mapped reads were assembled by Stringtie v1.3.3 (Pertea et al., 2015). Then, we merged the assembled transcripts into consensus transcripts using custom Python scripts after filtering reads with length less than 200 nt. Coffcompare v2.2.1 was applied to removed transcripts annotated as “C” and “=” (“C” for partial match, and “=” for full match) (Trapnell et al., 2010). Next, the remaining known protein-coding transcripts were removed by BLASTX and Hmmscan (Eddy, 2009). The Coding Potential Calculator (CPC) (Kong et al., 2007) was used to assess the coding potential of the remaining transcripts; transcripts with FPKM > 0 (where FPKM = fragments per kilobase of transcript per million mapped reads) in at least one biological replicate were annotated as lncRNAs. The pipeline of lncRNA identification is depicted in Supplementary Figure S1.
Classification of lncRNAs
The identified lncRNAs were classified by FEELnc (Wucher et al., 2017), a tool for lncRNA prediction and annotation. According to the location with corresponding gene, lncRNAs were classified into two categories: intergenic lncRNAs and intragenic lncRNAs. The intragenic lncRNAs can be further subclassified into four categories: (1) sense intronic lncRNAs, (2) antisense intronic lncRNAs, (3) sense exonic lncRNAs, and (4) antisense exonic lncRNAs.
Transcript Expression Level Analysis
Transcript abundance was evaluated by FPKM values of mRNAs and lncRNAs using StringTie v1.3.3 (Pertea et al., 2015). The transcripts that exhibited FPKM > 0.1 were used to filter the expressed protein-coding genes. Subsequently, log2-transformed values of FPKM + 1 were used for further analysis. For reliability of the experimental data, Spearman correlations were checked across the four developmental stages. t-Distributed stochastic neighbor embedding (t-SNE) was performed using R software and hierarchical clustering was carried out using MultiExperiment Viewer (MeV v4.9.0) (Howe et al., 2010). The differentially expressed genes (DEGs) and lncRNAs transcripts (DE-lncRNAs) were identified using edgeR package1 based on the read count data (Robinson et al., 2009). The significant DEGs and DE-lncRNAs were screened with a false discovery rate <0.05 and | fold change (FC)| ≥ 1.5 as cutoff.
Time-Series Analysis
We performed the four time-series analysis using the Short Time-series Expression Miner (STEM) (Ernst and Bar-Joseph, 2006). Different colors represented significantly enriched model profiles (Bonferroni-adjusted P-values ≤ 0.05), where model profiles with the same color were considered the same cluster of expression profiles.
LncRNA-mRNA Co-expression Network Analysis
To predict the function of lncRNAs, a co-expression network was constructed by Weighted Gene Co-expression Network Analysis (WGCNA) in the R package (Langfelder and Horvath, 2008). We calculated the co-expression network based on the DEGs and DE-lncRNAs that were screened by edgeR. The “co-expression modules” detected by average linkage hierarchical clustering represented genes with coordinated expression patterns.
Functional Enrichment Analysis for Genes
The DEGs were subjected to functional enrichment analysis for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway categories using Metascape (Zhou et al., 2019) online tool2. The GO terms and KEGG pathways with P-values < 0.01 were considered significant.
Quantitative Real-Time PCR Validation
To verify the reliability of the expression profiles of RNA-sequencing data, we randomly evaluated expression of five each for DEGs and DE-lncRNAs by quantitative real-time PCR (qPCR). The primers used for qPCR are shown in Supplementary Table S1. The cDNA was synthesized using the EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech, Beijing, China) following the manufacturer’s recommendation. We performed qPCR on a CFX96 Real-Time PCR Detection system (Bio-Rad Laboratories, Hercules, CA, United States) using TransStart Top Green qPCR SuperMix (Transgen Biotech, Beijing, China). The reaction volume was 20 μl, containing 1 μl of template cDNA, 0.4 μl of forward and reverse primers, 10 μl of Top Green qPCR SuperMix, and 8.2 μl of RNase-free water. The qPCR measurements were evaluated in triplicate followed the qPCR system: 94°C for 30 s, followed by 40 cycles of 94°C for 5 s and 30 s at the Tm indicated in Supplementary Table S1. The melting curve analysis was carried out from 65 to 95°C with increments of 0.5°C for 5 s. The expression levels of each transcript were normalized to that of ACTB used as endogenous control. The fold change in expression levels was determined using the 2–ΔΔCt method.
Results
Data Summary of lncRNAs and mRNAs Profile in Chicken Liver
To systematically explore the hepatic transcriptome during development in chicken, we constructed 12 cDNA libraries from 12 female Luhua chickens in four crucial stages during development: JP (day 60), SM (day 133), PL (day 220), and BP (day 400). The strand-specific libraries were sequenced on Illumina’s HiSeq platform. A total of 195.53 Gb of raw data with 150-bp paired-end sequences were generated. After filtering adaptors and low-quality reads, 192.31 Gb clean reads were retained corresponding to an average of ∼16.03 Gb of high-quality reads per sample. The high-quality reads were mapped to the chicken genome (G. gallus-5.0) with a mapping rate of 92.20 to 96.26% using HISAT2.1.0 (Supplementary Table S2).
Genomic Characterization and Classification of lncRNAs
To identify putative lncRNAs in chicken liver, we analyzed the homology of transcripts with known genes, inclusion of a known protein-coding region, and their coding potential. As a result, we obtained 16,930 putative lncRNAs that were expressed in at least one biological replicate (FPKM > 0), and 18,260 mRNAs that were substantially expressed (FPKM > 0.1 in at least one replicate) (Supplementary Tables S3, S4). We found that 14,959 (88.36%) lncRNAs and 15,983 (87.53%) mRNAs were assembled on chromosomes of chicken (Figure 1A and Supplementary Table S5), and most of these transcripts tended to be distributed in the large chromosomes (nos. 1–5 and the Z chromosome). Additionally, the numbers of mRNAs and lncRNAs were highly correlated with chromosome size (r = 0.9712 and 0.9371, respectively) (Supplementary Figure S2). Further comparative analysis showed that lncRNAs were shorter with fewer exons and were expressed at a lower level than that of mRNAs (Figures 1B–D).
Figure 1. Genomic characterization of lncRNAs and mRNAs. (A) Chromosome distribution of lncRNAs and mRNAs identified in chicken liver. (B) The distribution of transcript length of lncRNAs and mRNAs. (C) The proportion of transcripts with different number of exons. (D) The expression level of mRNAs and lncRNAs in four developmental stages. (E) The classification of lncRNAs.
According to their genomic location, we used FEELnc to separate the 14,730 identified lncRNAs into five groups: 9091 intergenic lncRNAs (lincRNAs, 53.7%), 2489 sense exonic overlapping lncRNAs (14.7%), 1247 antisense exonic overlapping lncRNAs (7.4%), 734 sense intronic overlapping lncRNAs (4.3%), and 1169 antisense intronic overlapping lncRNAs (6.9%) (Figure 1E). The results of functional enrichment showed that genes whose exons and introns overlapped with lncRNAs with the most significantly enriched GO categories were annotated with terms that involved development and cell activities, such as regulation of system process (GO:0044057), developmental growth (GO:0048589), regulation of cell adhesion (GO:0030155), and cell surface receptor signaling pathway involved in cell–cell signaling (GO:1905114). KEGG pathways were enriched with terms such as MAPK signaling pathway (hsa04010), endocytosis (hsa04144), and prolactin signaling pathway (hsa04917) (Supplementary Table S6). These results indicated that the potential functions of exon- and intron-overlapping lncRNAs may be closely involved in these enriched categories and pathways.
Expression Profiles of mRNAs and lncRNAs
Based on the expression change of mRNAs and lncRNAs, we calculated Spearman correlations between each pair of samples. The results showed that both the mRNAs and lncRNAs of 12 samples were grouped into four clusters based on their developmental stages: samples at the PL and JP stages were separated from the SM and BP stages, whereas the SM and BP stages were more convergent than the PL and JP stages (Figures 2A,B). Similarly, t-SNE plots of 12 samples indicated that mRNAs and lncRNAs were expressed in a stage-specific manner (Figures 2C,D). In addition, we found that the correlation between two consecutive stages of lncRNAs was weaker than that for the mRNAs (Figures 2A,B). This observation indicates that lncRNAs have a narrower time windows on expression profiles and more restricted temporal expression than mRNAs. To further validate the inference, we used the Shannon entropy (H) value to calculate the specificity of expression profile across developmental stages. Compared with mRNAs, all three types of lncRNAs, intergenic, exon-overlapping, and intron-overlapping, showed increased temporal specificity (Figure 2E).
Figure 2. The expression profiles of mRNAs and lncRNAs. (A,B) Dynamic changes in expression profile of mRNAs and lncRNAs. The top and left panel is the sample tree and the value represents the pairwise Spearman correlation. (C,D) Two-way t-SNE plot of mRNAs and lncRNAs based on expression profiles. (E) Distributions of Shannon entropy-based temporal specificity scores were calculated for distinct types of lncRNAs and mRNAs.
Identification of Differentially Expressed mRNAs and lncRNAs
To identify the DEGs and DE-lncRNAs in chicken liver, we performed pairwise comparisons of expression level between the four developmental stages with a threshold of |FC| > 1.5 and a false discovery rate < 0.05. Overall, by comparing six pairs, we have identified 4405 unique mRNAs and 752 unique lncRNAs, which were significantly differentially expressed over four developmental stages. A greater number of DEGs were obtained from JP versus SM, followed by JP versus BP and SM versus PL, whereas most DE lncRNAs were detected at SM versus PL, followed by JP versus SM and JP versus BP (Supplementary Tables S7, S8). The results of functional enrichment showed that most of these DEGs were enriched in metabolic processes, hormone stimulus, and developmental processes, such as fatty acid metabolic process (GO:0006631), oxidoreductase activity (GO:0016491), regulation of hormone levels (GO:0010817), and reproductive structure development (GO:0048608) (Supplementary Table S9).
Subsequently, we analyzed the common DEGs and DE lncRNAs for the four developmental stages of chicken liver. The Venn diagram showed that only three genes were differentially expressed in all four developmental stages, and no common DE lncRNAs were identified in these stages (Figures 3A,B). The three common DEGs were vitellogenin 2 (VTG2), riboflavin binding protein (RBP), and a novel gene (ENSEMBL_Id: ENSGALG00000034504). Interestingly, we found that all three protein-coding genes had the same pattern of expression; that is, expression increased during the first three developmental stages, peaked at the PL stage, and then declined at the BP stage (Figure 3C), indicating that these genes contributed to lipoprotein synthesis and egg yolk formation for meeting demand for egg production at the PL stage.
Figure 3. Differentially expressed mRNAs and lncRNAs. The Venn diagrams of common differentially expressed mRNAs (A) and lncRNAs (B) in four developmental stages. (C) Dynamic expression profiles of VTG2, RBP, and a novel gene (ENSGALG00000034504).
Time-Series and Co-expression Network Analysis of mRNAs and lncRNAs
To further investigate the global expression pattern of mRNAs and lncRNAs during four hepatic developmental stages in chicken, we conducted a time-series analysis using STEM. A total of 11,749 mRNAs and 994 lncRNAs were partitioned into 6 and 5 module profiles, respectively, comprising 10 and 8 significantly enriched (Bonferroni-adjusted P-value ≤ 0.05) model profiles for mRNAs and lncRNAs, respectively (Figure 4A and Supplementary Tables S10, S11). In these model profiles, we found that mRNAs and lncRNAs existed in some model profiles with similar dynamic expression patterns, suggesting that their functions were highly correlated. Subsequently, we used the mRNAs in the yellow module (monotonic trajectory with increased expression pattern across time) to analyze the functional distribution. Interestingly, the genes in the yellow module were mainly enriched in sulfur compound metabolic process (GO:0006790), lipid biosynthetic process (GO:0008610), cellular response to hormone stimulus (GO:0032870), and glucose homeostasis (GO:0042593) (Supplementary Table S12). This finding was consistent with the biological function of chicken liver in different periods.
Figure 4. Time-series model and co-expression networks of lncRNAs and mRNAs. (A) Time-series modules of mRNAs and lncRNAs. The top panel shows mRNAs and the second panel shows lncRNAs. Numbers in the top left corner indicate module number. Numbers in the lower left corner indicate number of mRNAs or lncRNAs in each module. The same color was used to represent each cluster. (B) Co-expression network constructed by the WGCNA method. The top panel shows cluster dendrogram obtained by average linkage hierarchical clustering. The low panel shows co-expression modules of lncRNAs and mRNAs with different colors. (C) Heat map showing the mRNAs corresponding to the largest three co-expression network modules of lncRNAs. Values in heat map represent log2(FPKM + 1) of each mRNA in each sample minus the mean value of each mRNA across all samples. (D) Functional categories of the mRNAs corresponding to the largest three co-expression network modules of lncRNAs.
Due to the characteristics of lncRNA and the lack of annotation libraries, it remains difficult to predict the function of lncRNAs. Here, we explored functional relatedness between mRNAs and lncRNAs using co-expression analysis. We built the co-expression network based on a non-redundant list of DEGs and DE-lncRNAs by WGCNA. Ultimately, 14 co-expression modules of lncRNAs and mRNAs were identified (Figure 4B). To predict the function of the DE-lncRNAs, we only considered the top three modules of lncRNA and these overlapping genes (Figure 4C and Supplementary Table S13), which accounted for 63.22% of the total genes in the 14 modules. We submitted these genes to Metascape for functional enrichment and the result showed that the co-expression genes in the top three modules were mainly enriched in terms related to lipids (Figure 4D and Supplementary Table S14), such as lipid biosynthetic process (GO:0008610), neutral lipid metabolic process (GO:0006638), fatty acid metabolic process (GO:0006631), and lipid catabolic process (GO:0016042). The KEGG analysis of the peroxisome proliferator-activated receptor (PPAR) signaling pathway (hsa03320) was also closely related to intracellular synthesis of lipid. These results indicated that most of the DE-lncRNAs may be involved in lipid synthesis and metabolism in the liver of chicken.
Quantitative Real-Time PCR Validation of Gene Expression
We selected five DEGs (VTG2, ApoVLDL2, APOB, FGB, and FASN) and five DE-lncRNAs (TCONS_00220501, TCONS_00179354, TCONS_00083873, TCONS_00155619, and TCONS_00037924) for validation by qPCR (Figure 5). The relative expression of the five DEGs and five DE lncRNA of each sample by qPCR was compared with the transformed log2(FPKM + 1) values of RNA-seq. All genes except FGB showed the same expression pattern in both qPCR and RNA-seq, having the highest expression level at the PL stage. The altered expression pattern of DEGs and DE-lncRNAs at most stages were also consistent with the RNA-seq data. These results illustrated the reliability of our RNA-seq data and confirmed the accuracy of the identified transcripts.
Figure 5. Validation of five expressed DEGs and DE lncRNAs by qPCR. The x-axis represents the four developmental stages. The y-axis indicates the relative expression of each gene; green lines are log2(FPKM + 1) values of RNA-seq and the blue lines show relative expression by qPCR.
Discussion
It is well established that the liver is the major organ of lipid biosynthesis in the chicken (Leveille et al., 1968; O’Hea and Leveille, 1968, 1969). In laying hens, the liver synthesizes a series of lipoproteins and most of the additional minor yolk precursors. This biological process represents a complex and highly orchestrated program of gene expression and regulation of various regulatory molecules, such as microRNAs (miRNAs) and lncRNAs. Although the mRNA and miRNA transcriptome in chicken liver has been partially explored in previous studies (Wang et al., 2014; Li et al., 2015, 2016), a systematic evaluation of the lncRNA repertoire across chicken liver developmental stages is lacking. To further investigate change in the transcriptome expression profile in different developmental stages, we performed RNA-seq and identified lncRNAs and mRNAs at days 60, 133, 220, and 400, representing four developmental stages of the chicken. In total, we obtained 192.31 Gb of clean data from 12 libraries and identified 16,930 lncRNAs and 18,260 mRNAs. We found that lncRNAs had several specific genomic characteristics compared with mRNAs (shorter sequences, fewer exons, lower expression level, and more temporal specificity) consistent with results in other studies (Derrien et al., 2012b; Peng et al., 2014; Muret et al., 2017; Jin et al., 2018).
According to our results, lncRNAs and mRNAs were expressed in a stage-dependent manner of expression profiles. Spearman correlation and t-SNE showed that the SM stage is more closed to BP stage, which means lncRNA and mRNA had similar expression profiles at the two stages. Indeed, the capacity of hepatic lipogenesis weakens because of the absence of egg production in the SM and BP stages, which is in line with the lipogenic enzyme activities involved in this process (Pearce, 1971). We found that lncRNAs were expressed for shorter duration compared with mRNAs, suggesting that lncRNAs have more restricted temporal expression than mRNAs (Pauli et al., 2012).
Pairwise comparisons of lncRNAs and mRNAs were performed to identify DEGs and DE-lncRNAs. We found that most of these DEGs were enriched in various lipid metabolic and developmental processes, especially between JP and PL stage, which were consistent with a study of Li et al. (2015). In addition, only three genes, VTG2, RBP, and a novel gene (ENSEMBL_Id: ENSGALG00000034504), were found to be differentially expressed in all four stages. These three DEGs had similar expression profiles, reaching peak expression at the PL stage, indicating their important roles involved in chicken liver development and egg production. It is clear that VTG2, a marker gene in yolk precursor formation, has a vital role in follicle development and maturation (Schneider, 1996). The VTG2 gene encodes vitellogenin, which is dependent on estrogen stimulation (Denslow et al., 1999; Li et al., 2014). In laying hens, three subtypes of VTG are identified in liver; the ratio of VTG I:II:III RNA is 4:10:1 (Evans et al., 1988). As a glycophospholipoprotein, VTG is synthesized in the liver and then deposited in the rapidly developing oocytes. It has a critical effect on reproduction performance by providing nutrients and functional substances for embryonic development (Richards, 1997). The protein encoded by RBP, which is stored in eggs, is also essential for maintaining growth and developing embryos. It binds riboflavin (vitamin B2) and is involved in its transport from the serum to the oocyte (White and Merrill, 1988). Changes in the expression levels of VTG2 and RBP in our study were consistent with the biological functions of chicken liver in different periods. Surprisingly, the two egg-expressed genes, VTG2 and RBP, were lost in mammals during evolution (Tian et al., 2010). Unlike the oviparous vertebrates that depend on nutritional reserves, mammals have evolved placenta to nourish their embryos and mammary gland as milk resources for early offspring (Brawand et al., 2008; Tian et al., 2010). In addition, we identified a novel protein-coding gene that was differentially expressed in liver at all four developmental stages, indicating that it may participate in lipid metabolism in chicken liver.
Time-series analysis suggested that both lncRNAs and mRNAs were dynamically expressed in liver during four developmental stages and have similar expression patterns. Genes in the yellow module with upregulated expression, such as ACACA (Cronan and Waldrop, 2002), CYP3A4 (Wang et al., 2008), and INSIG1 (Gong et al., 2006), were involved in metabolic activity and lipid homeostasis. At the PL stage, the liver is an important organ for the synthesis of egg yolk substances, and many genes related to lipid metabolism and synthesis of important yolk proteins are highly expressed to meet the demand for egg production. Besides, we found that some genes maintained a high expression level at the BP stage in the yellow module, such as PRLR, which was enriched in hormone stimulus-related process. The PRLR gene, which encodes prolactin receptor and modulates prolactin activity in an endocrine and autocrine manner, is generally accepted as a candidate gene for the onset and maintenance of broodiness in birds (Wilkanowska et al., 2014). Studies also have illustrated that this gene polymorphisms are associated with egg production, body weight at hatch, and age at sexual maturity (Rashidi et al., 2012; Zhang et al., 2012). The function of these genes can explain the changes in expression levels during corresponding stages. To determine the potential function of lncRNAs, we constructed co-expression networks between DE-lncRNAs and DEGs using WGCNA. The co-expression networks can be applied to associate lncRNAs with functionally annotated mRNAs (van Dam et al., 2017). Transcripts with similar expression patterns are often biologically related and can therefore be divided into the same module by WGCNA. Our results showed that five clusters of GO and KEGG terms, including lipid biosynthetic process, neutral lipid metabolic process, fatty acid metabolic process, lipid catabolic process, and PPAR signaling pathway were enriched in the co-expressed genes. For example, APOA1 encodes key components modulating lipid metabolism (Sorci-Thomas et al., 1989). Apolipoprotein A1 is the major protein of high-density lipoprotein (HDL), and it functions to regulate the reverse cholesterol transport (Segrest et al., 2000). The SCD gene encodes stearoyl CoA desaturase-1, which is located in endoplasmic reticulum membrane and transforms different saturated fatty acids into monounsaturated fatty acids (Ntambi and Miyazaki, 2004). These co-expressed genes play crucial roles in lipid biosynthesis and metabolism, suggesting the biological function of DE-lncRNAs by regulating hepatic lipogenesis and liver development in chicken.
Conclusion
In summary, we profiled a systematic lncRNA repertoire in liver across different developmental stages of the chicken. Our results indicated the potential roles of lncRNA in liver development based on co-expression network, but further research is need for function validation. Our study provides a high-quality resource for future transcriptomic studies and improves the comparative understanding of molecular mechanisms of liver development in chickens.
Data Availability Statement
The datasets generated for this study can be found in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) with accession number GSE138152.
Ethics Statement
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Sichuan Agricultural University.
Author Contributions
CN, YY, QZ, MY, and DL designed the study and drafted the manuscript. CN, PZ, ZX, XF, BZ, and DY carried out the feeding trial and collected the liver tissue samples. TM, SH, XZ, YW, HY, and YH extracted the liver tissue RNA and contributed to the data analysis. QN, YL, and MZ contributed to the qPCR validation. HX, QZ, and DL reviewed and improved this manuscript. All authors contributed to the manuscript and approved the submitted version.
Funding
This research was funded by the Sichuan Provincial Department of Science and Technology Program (2019JDTD0009), the Fok Ying-Tong Education Foundation for Yong Teachers in the High Education Institutions of China (161026), and the Independent Research Project of Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province (SNDK-ZZ201801).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Louise Adam, ELS(D), from Liwen Bianji, Edanz Editing, China (www.liwenbianji.cn/ac) for editing the English text of a draft of this manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00574/full#supplementary-material
FIGURE S1 | Identification pipeline of lncRNAs.
FIGURE S2 | Correlation between size of chromosome and number of identified transcripts: (A) mRNA and (B) lncRNAs.
FIGURE S3 | The number of differentially expressed genes (DEGs, A) and lncRNAs (DE lncRNAs, B) identified by pairwise comparison.
Footnotes
References
Bedu, E., Chainier, F., Sibille, B., Meister, R., Dallevet, G., Garin, D., et al. (2002). Increased lipogenesis in isolated hepatocytes from cold-acclimated ducklings. Am. J. Physiol. Regul. Integr. Comp. Physiol. 283, R1245–R1253. doi: 10.1152/ajpregu.00681.2001
Brawand, D., Wahli, W., and Kaessmann, H. (2008). Loss of egg yolk genes in mammals and the origin of lactation and placentation. PLoS Biol. 6:e63. doi: 10.1371/journal.pbio.0060063
Brockmöller, J., and Ivar, R. (1994). Assessment of liver metabolic function. Clin. Pharmacokinet. 27, 216–248. doi: 10.2165/00003088-199427030-00005
Calne, R. Y. (2000). Immunological tolerance – the liver effect. Immunol. Rev. 174, 280–282. doi: 10.1034/j.1600-0528.2002.017419.x
Carrieri, C., Cimatti, L., Biagioli, M., Beugnet, A., Zucchelli, S., Fedele, S., et al. (2012). Long non-coding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat. Nature 491, 454–457. doi: 10.1038/nature11508
Cronan, J. E., and Waldrop, G. L. (2002). Multi-subunit acetyl-CoA carboxylases. Prog. Lipid Res. 41, 407–435. doi: 10.1016/s0163-7827(02)00007-3
Ðaković, N., Térézol, M., Pitel, F., Maillard, V., Elis, S., Leroux, S., et al. (2014). The loss of adipokine genes in the chicken genome and implications for insulin metabolism. Mol. Biol. Evol. 31, 2637–2646. doi: 10.1093/molbev/msu208
Denslow, N. D., Chow, M. C., Kroll, K. J., and Green, L. (1999). Vitellogenin as a biomarker of exposure for estrogen or estrogen mimics. Ecotoxicology 8, 385–398. doi: 10.1023/a:1008986522208
Derrien, T., Guigó, R., and Johnson, R. (2012a). The long non-coding RNAs: a new (P)layer in the “dark matter”. Front. Genet. 2:107. doi: 10.3389/fgene.2011.00107
Derrien, T., Johnson, R., Bussotti, G., Tanzer, A., Djebali, S., Tilgner, H., et al. (2012b). The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 22, 1775–1789. doi: 10.1101/gr.132159.111
Djebali, S., Davis, C. A., Merkel, A., Dobin, A., Lassmann, T., Mortazavi, A., et al. (2012). Landscape of transcription in human cells. Nature 489, 101–108. doi: 10.1038/nature11233
Eddy, S. R. (2009). A new generation of homology search tools based on probabilistic inference. Genome Inform. 23, 205–211.
Ernst, J., and Bar-Joseph, Z. (2006). STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics 7:191. doi: 10.1186/1471-2105-7-191
Esteller, M. (2011). Non-coding RNAs in human disease. Nat. Rev. Genet. 12, 861–874. doi: 10.1038/nrg3074
Evans, M. I., Silva, R., and Burch, J. B. (1988). Isolation of chicken vitellogenin I and III cDNAs and the developmental regulation of five estrogen-responsive genes in the embryonic liver. Genes Dev. 2, 116–124. doi: 10.1101/gad.2.1.116
Gong, C., and Maquat, L. E. (2011). lncRNAs transactivate STAU1-mediated mRNA decay by duplexing with 3’ UTRs via Alu elements. Nature 470, 284–288. doi: 10.1038/nature09701
Gong, Y., Lee, J. N., Lee, P. C. W., Goldstein, J. L., Brown, M. S., and Ye, J. (2006). Sterol-regulated ubiquitination and degradation of Insig-1 creates a convergent mechanism for feedback control of cholesterol synthesis and uptake. Cell Metab. 3, 15–24. doi: 10.1016/j.cmet.2005.11.014
Grant, D. M. (1991). “Detoxification pathways in the liver,” in Journal of Inherited Metabolic Disease, eds R. A. Harkness, R. J. Pollitt, and G. M. Addison (Dordrecht: Springer), 421–430. doi: 10.1007/978-94-011-9749-6_2
Howe, E., Holton, K., Nair, S., Schlauch, D., Sinha, R., and Quackenbush, J. (2010). “MeV: multiexperiment viewer,” in Biomedical Informatics for Cancer Research, eds M. F. Ochs, J. T. Casagrande, and R. V. Davuluri (Boston, MA: Springer), 267–277. doi: 10.1007/978-1-4419-5714-6_15
Hussain, M. M., Rava, P., Pan, X., Dai, K., Dougan, S. K., Iqbal, J., et al. (2008). Microsomal triglyceride transfer protein in plasma and cellular lipid metabolism. Curr. Opin. Lipidol. 19, 277–284. doi: 10.1097/mol.0b013e3282feea85
Hussain, M. M., Rava, P., Walsh, M., Rana, M., and Iqbal, J. (2012). Multiple functions of microsomal triglyceride transfer protein. Nutr. Metab. 9:14. doi: 10.1186/1743-7075-9-14
Ivessa, N. E., Rehberg, E., Kienzle, B., Seif, F., Hermann, R., Hermann, M., et al. (2013). Molecular cloning, expression, and hormonal regulation of the chicken microsomal triglyceride transfer protein. Gene 523, 1–9. doi: 10.1016/j.gene.2013.03.102
Jin, L., Hu, S., Tu, T., Huang, Z., Tang, Q., Ma, J., et al. (2018). Global long noncoding RNA and mRNA expression changes between prenatal and neonatal lung tissue in pigs. Genes 9:443. doi: 10.3390/genes9090443
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317
Kirchgessner, T. G., Heinzmann, C., Svenson, K. L., Gordon, D. A., Nicosia, M., Lebherz, H. G., et al. (1987). Regulation of chicken apolipoprotein B: cloning, tissue distribution, and estrogen induction of mRNA. Gene 59, 241–251. doi: 10.1016/0378-1119(87)90332-5
Kong, L., Zhang, Y., Ye, Z.-Q., Liu, X.-Q., Zhao, S.-Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35(Suppl._2), W345–W349. doi: 10.1093/nar/gkm391
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559
Leveille, G. A., O’Hea, E. K., and Chakrabarty, K. (1968). In vivo lipogenesis in the domestic chicken. Proc. Soc. Exp. Biol. Med. 128, 398–401. doi: 10.3181/00379727-128-33022
Leveille, G. A., Romsos, D. R., Yeh, Y.-Y., and O’Hea, E. K. (1975). Lipid biosynthesis in the chick. A consideration of site of synthesis, influence of diet and possible regulatory mechanisms1. Poult. Sci. 54, 1075–1093. doi: 10.3382/ps.0541075
Li, H., Ma, Z., Jia, L., Li, Y., Xu, C., Wang, T., et al. (2016). Systematic analysis of the regulatory functions of microRNAs in chicken hepatic lipid metabolism. Sci. Rep. 6:31766. doi: 10.1038/srep31766
Li, H., Wang, T., Xu, C., Wang, D., Ren, J., Li, Y., et al. (2015). Transcriptome profile of liver at different physiological stages reveals potential mode for lipid metabolism in laying hens. BMC Genomics 16:763. doi: 10.1186/s12864-015-1943-0
Li, J., Leghari, I. H., He, B., Zeng, W., Mi, Y., and Zhang, C. (2014). Estrogen stimulates expression of chicken hepatic vitellogenin II and very low-density apolipoprotein II through ER-α. Theriogenology 82, 517–524. doi: 10.1016/j.theriogenology.2014.05.003
Liu, G., Mattick, J., and Taft, R. J. (2013). A meta-analysis of the genomic and transcriptomic composition of complex life. Cell Cycle 12, 2061–2072. doi: 10.4161/cc.25134
Maamar, H., Cabili, M. N., Rinn, J., and Raj, A. (2013). linc-HOXA1 is a noncoding RNA that represses Hoxa1 transcription in cis. Genes Dev. 27, 1260–1271. doi: 10.1101/gad.217018.113
Mao, Y. S., Sunwoo, H., Zhang, B., and Spector, D. L. (2010). Direct visualization of the co-transcriptional assembly of a nuclear body by noncoding RNAs. Nat. Cell Biol. 13, 95–101. doi: 10.1038/ncb2140
Martianov, I., Ramadass, A., Serra Barros, A., Chow, N., and Akoulitchev, A. (2007). Repression of the human dihydrofolate reductase gene by a non-coding interfering transcript. Nature 445, 666–670. doi: 10.1038/nature05519
Mashek, D. G. (2013). Hepatic fatty acid trafficking: multiple forks in the road. Adv. Nutr. 4, 697–710. doi: 10.3945/an.113.004648
Muret, K., Klopp, C., Wucher, V., Esquerré, D., Legeai, F., Lecerf, F., et al. (2017). Long noncoding RNA repertoire in chicken liver and adipose tissue. Genet. Sel. Evol. 49:6.
Ng, S.-Y., Bogu, G. K., Soh, B. S., and Stanton, L. W. (2013). The long noncoding RNA RMST interacts with SOX2 to regulate neurogenesis. Mol. Cell 51, 349–359. doi: 10.1016/j.molcel.2013.07.017
Nguyen, P., Leray, V., Diez, M., Serisier, S., Bloc’h, J. L., Siliart, B., et al. (2008). Liver lipid metabolism. J. Anim. Physiol. Anim. Nutr. 92, 272–283. doi: 10.1111/j.1439-0396.2007.00752.x
Nimpf, J., and Schneider, W. J. (1991). Receptor-mediated lipoprotein transport in laying hens. J. Nutr. 121, 1471–1474. doi: 10.1093/jn/121.9.1471
Ntambi, J. M., and Miyazaki, M. (2004). Regulation of stearoyl-CoA desaturases and role in metabolism. Prog. Lipid Res. 43, 91–104. doi: 10.1016/s0163-7827(03)00039-0
O’Hea, E. K., and Leveille, G. A. (1968). Lipogenesis in isolated adipose tissue of the domestic chick (Gallus domesticus). Comp. Biochem. Physiol. 26, 111–120. doi: 10.1016/0010-406x(68)90317-4
O’Hea, E. K., and Leveille, G. A. (1969). Lipid biosynthesis and transport in the domestic chick (Gallu domesticus). Comp. Biochem. Physiol. 30, 149–159. doi: 10.1016/0010-406x(69)91309-7
Pauli, A., Valen, E., Lin, M. F., Garber, M., Vastenhouw, N. L., Levin, J. Z., et al. (2012). Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res. 22, 577–591. doi: 10.1101/gr.133009.111
Pearce, J. (1971). An investigation of lipogenic and glycolytic enzyme activity in the liver of sexually immature and mature domestic fowl. Biochem. J. 123, 717–719. doi: 10.1042/bj1230717
Peng, L., Paulson, A., Li, H., Piekos, S., He, X., Li, L., et al. (2014). Developmental programming of long non-coding RNAs during postnatal liver maturation in mice. PLoS One 9:e114917. doi: 10.1371/journal.pone.0114917
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Rashidi, H., Rahimi-Mianji, G., Farhadi, A., and Gholizadeh, M. (2012). Association of prolactin and prolactin receptor gene polymorphisms with economic traits in breeder hens of indigenous chickens of Mazandaran province. Iran. J. Biotechnol. 10, 129–135.
Richards, M. P. (1997). Trace mineral metabolism in the avian embryo. Poult. Sci. 76, 152–164. doi: 10.1093/ps/76.1.152
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Rui, L. (2014). Energy metabolism in the liver. Compr. Physiol. 4, 177–197. doi: 10.1002/cphy.c130024
Schneider, W. J. (1995). Yolk precursor transport in the laying hen. Curr. Opin. Lipidol. 6, 92–96. doi: 10.1097/00041433-199504000-00006
Schneider, W. J. (1996). “Vitellogenin receptors: oocyte-specific members of the low-density lipoprotein receptor supergene family,” in International Review of Cytology, ed. K. W. Jeon (Cambridge, MA: Academic Press), 103–137. doi: 10.1016/s0074-7696(08)62507-3
Segrest, J. P., Li, L., Anantharamaiah, G. M., Harvey, S. C., Liadaki, K. N., and Zannis, V. (2000). Structure and function of apolipoprotein A-I and high-density lipoprotein. Curr. Opin. Lipidol. 11, 105–115.
Sorci-Thomas, M., Prack, M. M., Dashti, N., Johnson, F., Rudel, L. L., and Williams, D. L. (1989). Differential effects of dietary fat on the tissue-specific expression of the apolipoprotein A-I gene: relationship to plasma concentration of high density lipoproteins. J. Lipid Res. 30, 1397–1403.
Stern, C. D. (2005). The chick: a great model system becomes even greater. Dev. Cell 8, 9–17. doi: 10.1016/s1534-5807(04)00425-3
Tian, X., Gautron, J., Monget, P., and Pascal, G. (2010). What makes an egg unique? Clues from evolutionary scenarios of egg-specific genes. Biol. Reprod. 83, 893–900. doi: 10.1095/biolreprod.110.085019
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
van Dam, S., Võsa, U., van der Graaf, A., Franke, L., and de Magalhães, J. P. (2017). Gene co-expression analysis for functional classification and gene–disease predictions. Brief. Bioinform. 19, 575–592. doi: 10.1093/bib/bbw139
Walzem, R. L., Hansen, R. J., Williams, D. L., and Hamilton, R. L. (1999). Estrogen induction of VLDLy assembly in egg-laying hens. J. Nutr. 129, 467S–472S. doi: 10.1093/jn/129.2.467S
Wang, K., Chen, S., Xie, W., and Wan, Y.-J. Y. (2008). Retinoids induce cytochrome P450 3A4 through RXR/VDR-mediated pathway. Biochem. Pharmacol. 75, 2204–2213. doi: 10.1016/j.bcp.2008.02.030
Wang, X., Yang, L., Wang, H., Shao, F., Yu, J., Jiang, H., et al. (2014). Growth hormone-regulated mRNAs and miRNAs in chicken hepatocytes. PLoS One 9:e112896. doi: 10.1371/journal.pone.0112896
White, H. B., and Merrill, A. H. (1988). Riboflavin-binding proteins. Annu. Rev. Nutr. 8, 279–299. doi: 10.1146/annurev.nu.08.070188.001431
Wilkanowska, A., Mazurowski, A., Mroczkowski, S., and Kokoszyński, D. (2014). Prolactin (PRL) and prolactin receptor (PRLR) genes and their role in poultry production traits. Folia Biol. 62, 1–8. doi: 10.3409/fb62_1.1
Wiskocil, R., Bensky, P., Dower, W., Goldberger, R. F., Gordon, J. I., and Deeley, R. G. (1980). Coordinate regulation of two estrogen-dependent genes in avian liver. Proc. Natl. Acad. Sci. U.S.A. 77, 4474–4478. doi: 10.1073/pnas.77.8.4474
Wucher, V., Legeai, F., Hédan, B., Rizk, G., Lagoutte, L., Leeb, T., et al. (2017). FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 45:e57. doi: 10.1093/nar/gkw1306
Yang, L., Lin, C., Liu, W., Zhang, J., Ohgi, K. A., Grinstein, J. D., et al. (2011). ncRNA- and Pc2 methylation-dependent gene relocation between nuclear structures mediates gene activation programs. Cell 147, 773–788. doi: 10.1016/j.cell.2011.08.054
Zhang, L., Li, D., Liu, Y., Wang, Y., Zhao, X., and Zhu, Q. (2012). Genetic effect of the prolactin receptor gene on egg production traits in chickens. Genet. Mol. Res. 11, 4307–4315. doi: 10.4238/2012.October.2.1
Keywords: long non-coding RNA, mRNA, liver, chicken, development
Citation: Ning C, Ma T, Hu S, Xu Z, Zhang P, Zhao X, Wang Y, Yin H, Hu Y, Fan X, Zeng B, Yang M, Yang D, Ni Q, Li Y, Zhang M, Xu H, Yao Y, Zhu Q and Li D (2020) Long Non-coding RNA and mRNA Profile of Liver Tissue During Four Developmental Stages in the Chicken. Front. Genet. 11:574. doi: 10.3389/fgene.2020.00574
Received: 09 December 2019; Accepted: 11 May 2020;
Published: 16 June 2020.
Edited by:
Jiuzhou Song, University of Maryland, College Park, United StatesReviewed by:
Guirong Sun, Henan Agricultural University, ChinaQinghua Nie, South China Agricultural University, China
Carl Joseph Schmidt, University of Delaware, United States
Copyright © 2020 Ning, Ma, Hu, Xu, Zhang, Zhao, Wang, Yin, Hu, Fan, Zeng, Yang, Yang, Ni, Li, Zhang, Xu, Yao, Zhu and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qing Zhu, zhuqingsicau@163.com; Diyan Li, diyanli@sicau.edu.cn
†These authors have contributed equally to this work