Transcriptome profiling of the liver among the prenatal and postnatal stages in chickens.

The liver is an important organ that has pivotal functions in the synthesis of several vital proteins, the metabolism of various biologically useful materials, the detoxification of toxic substances, and immune defense. Most liver functions are not mature at a young age and many changes happen during postnatal liver development, which lead to differential functions of the liver at different developmental stages. However, the transcriptome details of what changes occur in the liver after birth and the molecular mechanisms for the regulation of the developmental process are not clearly known in chickens. Here, we used RNA-sequencing to analyze the transcriptome of chicken liver from the prenatal (at an embryonic day of 13) to the postnatal stages (at 5 wk and 42 wk of age). A total of approximately 161.17 Gb of raw data were obtained, with 4,127 putative and 539 differentially expressed lncRNAs, and with 13,949 putative and 6,370 differentially expressed mRNAs. Coexpression of lncRNAs-mRNAs in hepatic transcriptome analysis showed that the liver plays important roles in providing energy for organisms through the mitochondrial respiratory chain in chickens, meanwhile, acting as a crucial part of antioxidant stress. The developmental transcriptome date revealed that antioxidant defenses are likely to act on chicken embryo development and that significant functional changes during postnatal liver development are associated with the liver maturation of chickens. These results provide a timeline for the functional transcriptome transition from the prenatal to adult stages in chickens and will be helpful to reveal the underlying molecular mechanisms of liver development.

INTRODUCTION

The chicken (Gallus gallus domesticus), a subspecies of the red jungle fowl, is one of the most common and widespread domestic animals. Humans keep chickens primarily as a source of food (consuming both their meat and eggs) (Burt, 2005; Burt, 2007). The liver is found in all vertebrates and is typically the largest visceral organ. Its principal functions in metabolism include the decomposition of red blood cells, the regulation of glycogen storage, and the production of hormones (Brockmoller and Roots, 1994). In addition, the liver participates in protein synthesis, production of biochemical necessary for digestion, and is responsible for detoxification and immunological effects (Calne, 2000; Chen et al., 2009). The liver is a mature organ in the adult with important, well-defined roles in the maintenance of bile acid synthesis, nutrient homeostasis, the metabolism of xenobiotics and endogenous hormones, and the detoxification of exogenous compounds (Grijalva and Vakili, 2013). However, in the developing embryonic period and in the newborn stage, the liver acts as a hematopoietic organ with a key function in generating blood cells (Guo et al.,2009).

Embryo completion occurs from day 8 to day 14 after early embryogenesis in chicken embryonic development, and the liver produced the major plasma proteins during fetal life, such as α-fetoprotein (Réhault-Godbert et al., 2014). During the postnatal developmental periods, a transition occurs in the liver in order for it to gain mature metabolic activity (Peng et al., 2017). Significant changes in gene expression take place during the prenatal and postnatal stages in order to facilitate the switch from a hepatic microenvironment supportive of hematopoiesis to one that promotes a broad range of metabolic and detoxifying functions (Lee et al., 2012). However, the genetic background of changes in the liver physiology for each developmental period needs to be further studied, especially in embryonic development stage, rapid growth period, and laying peak of chickens.

High-throughput sequencing confirmed that long noncoding RNAs (lncRNAs) lncRNA plays an important role in the regulation of physiological metabolism. New, large-scale advances in livestock have been achieved by employing RNA-Seq in recent lncRNA studies, such as for the trout (Altobasei et al., 2016), cow (Koufariotis et al., 2015), pig (Ramayocaldas et al., 2012), sheep (Yue et al., 2016), chicken (Bourin et al., 2012), and duck (Ren et al., 2017). In chickens, most of the RNA-Seq studies have been focused on the protein coding RNAs and ignored the non-protein coding RNAs, such as miRNAs and lncRNAs. Due to the complexity of liver metabolic development, here, we performed rRNA-depleted RNA-Seq to explore the hepatic transcriptome profile of hens in 3 crucial liver developmental stages including embryos under day 13 of incubation (E13), pre-laying pullets at the age of 5 wk (W5) and laying hens with peak production at the age of 42 wk (W42). We also intended to identify functional transcripts (mRNAs) and regulatory transcripts (lncRNAs), and explored the differentially expressed transcripts across different liver developmental stages. In addition, we predicted the potential functions of lncRNAs by construction of the lncRNA-mRNA coexpression network. These findings would be of vital significance in contributing to a comprehensive understanding of transcriptome changes in chicken liver developmental stages.

MATERIALS AND METHODS

Chicken Embryo Incubation, Rearing, and Sampling

A total of 30 fertilized eggs laid by Roman hens were hatched in a microcomputer-regulated automatic incubator at 38°C from embryonic day 1 to 18 and at 37.5–37.8°C from day 19 to 21, with the humidity was maintained at 50–65%. The liver samples were collected from 3 well developed eggs at embryonic day 13 (E13). Hatched chickens were raised at the experimental farm of Zhejiang Academy of Agricultural Sciences. Three healthy chickens were sacrificed and liver tissues were collected at the age of 5 wk (W5) and 42 wk (W42), respectively. All samples were snap-frozen in liquid nitrogen and stored at −80°C until further analysis.

All protocols in this study were approved by the Institutional Animal Care and Use Committee of Zhejiang Academy of Agricultural Sciences.

RNA Extraction, Library Construction, and Sequencing

Total RNA was isolated from 9 liver tissues using TRIzol reagent (Invitrogen) following the manufacturer's protocols. We used a 1% agarose gel electrophoresis and Agilent 2100 Bioanalyzer (Agilent Technologies) to examine the RNA integrity and quality, respectively. The quality of RNA with an RNA Integrity Number (RIN) ≥7.5 and 28S/18S ≥1.0 were supposed to meet the high-throughput sequencing requirement. A Ribo-Zero™ Gold Kit (Illumina, San Diego, CA) was used to deplete rRNA for strand-specific library construction, and libraries were sequenced on Illumina's HiSeq platform with PE (paired-end) 150 bp sequencing.

Mapping and Assembly of lncRNA

Raw reads of quality were controlled using FastQC, and the high quality reads remaining after removing the low-quality reads and trimming the adaptor sequences were mapped to the chicken genome (Gallus_gallus-5.0) using Hisat (Kim et al., 2015). The transcriptome for each sample was assembled from the Hisat mapped reads separately by Stringtie v1.3.3 (Pertea et al., 2015). After filtering reads with length less than 200 nt, custom Python scripts (Iyer et al., 2015) were used to merge stringent transcripts into a consensus. Then, the transcripts annotated as “c” and “=” were removed by Cuffcompare v2.2.1 (Trapnell et al., 2010). Finally, the read coverage and fragments per kilobase of transcript per million mapped reads (FPKM) values were estimated by Cufflinks v2.2.1 (Trapnell et al., 2010) and Stringtie v1.3.3 (Pertea et al., 2015).

lncRNA Identification and Classification

First, we removed transcripts contained a known protein-coding domain using Hmmscan (Eddy, 2011) and BLASTX (Altschul et al., 1990). Second, the Coding Potential Calculator (CPC) (Kong et al., 2007) was used to assess the coding potential of the remaining transcripts and transcripts with CPC >0 were filtered. Finally, the remaining transcripts with FPKM >0 in more than 1 biological replicate were annotated as lncRNAs. The identified lncRNAs were classified based on their genomic location, with corresponding mRNAs through the FEELnc program (Muret et al., 2017); this methodology allowed us to distinguish intragenic and intergenic lncRNAs, where the intragenic lncRNAs were subclassified into 4 categories according to their orientation with respect to a reference set of genes.

Differential Expression Analysis

The expression levels of mRNAs were quantified with Cufflinks version 2.2.1 (Trapnell et al., 2010), and those of lncRNA transcripts were evaluated with Stringtie version 1.3.3. Genes with FPKM values less than 0.1 were removed (Lonsdale et al., 2015), and log2 transformed values of (FPKM+1) were used in subsequent analyses. Pearson correlations were calculated for all samples, and hierarchical clustering was performed using MultiExperiment Viewer (MeV version 4.9.0) (Saeed et al., 2003). Differential expression analyses were performed to detect the differentially expressed mRNAs and lncRNAs by MultiExperiment Viewer based on log2-transformed values of (FPKM+1). Benjamini adjusted P -values ≤0.05 were identified as differentially expressed genes (DEGs). Principal component analysis (PCA) and Pearson correlations were also performed as expression analysis.

Alternative Splicing Level Estimation

We identified 1:1 orthologous exons and estimated the percent spliced in (PSI) values following the methods described by Barbosamorais (Barbosamorais et al., 2005). We also examined the expression of all alternatively spliced transcripts.

Coexpression Analysis and Functions Prediction

To predict the putative function of lncRNAs, DEGs screened out by DESeq2 were further used to perform coexpression analysis through WGCNA (Langfelder and Horvath, 2008), which has been implemented in R and recently applied to the RNA-Seq data analysis; it was developed to analyze transcriptomic profiling experiments by quantifying the correlations between genes and transcripts that share similar expression patterns. The actual connectivity of features in the network was indicated by their position in a dendrogram or other network diagram, and genes with similar features were clustered into coexpressed “modules”. Functional annotation enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were conducted with the DAVID (Dennis et al., 2003). GO terms and KEGG pathways with P -values ≤0.05 were considered significantly enriched.

Validation of the Gene Expression Profile by qPCR

Five DEGs were randomly selected to confirm the expression profiles of RNA-Seq using qPCR. cDNA was synthesized using the EasyScript First-Strand cDNA Synthesis SuperMix (Transgen Biotech, Beijing) following the manufacturer's instructions. Primers were designed by Primer-BLAST, and the primer sequences were listed in Supplementary Table S1, and the β-actin gene was used as an internal control as previously described (Liu et al., 2019). The qPCR reactions were performed on the Bio-Rad iQ5 Real-Time PCR Detection system to detect RNA expression using TransStart Tip Green qPCR SuperMix (Transgen Biotech, Beijing). The amplification efficiency of qPCR ranged from 95 to 105%. Each sample was conducted with 3 technique replicates. The relative expression levels were calculated using the 2−△△Ct method (Livak and Schmittgen, 2001).

RESULTS

Data Summary of Hepatic Transcriptome in Chickens

To explore the transcriptome repertoire of chicken liver, we constructed 9 cDNA libraries focused on 3 crucial stages during development (embryonic day 13, at the age of 5 wk, and at 42 wk). Libraries were sequenced using the Illumina HiSeq platform and 150 bp paired-end reads were generated. All sequencing data have been submitted to the NCBI Gene Expression Omnibus with the accession number GSE121019. As a result, we obtained a total of approximately 161.17 Gb of raw data, corresponding to an average of approximately 17.56 Gb of clean data per sample. After removing adaptors and low quality reads, we finally obtained 158.00 Gb of clean data, with the percentage of clean reads in each library ranging from 97.54 to 98.50%. The proportion of reads with a Phred quality value of more than 30 (Q30) among the clean data ranged from 91.93 to 94.71%, and the overall average alignment rate of each sample was as high as 95.91% (Table 1).

Table 1

A summary of raw data, clean data, clean reads rate, and clean Q30 bases rate of sequenced samples.

Sample ID	Raw data (Gb)	Clean data (Gb)	Clean reads rate (%)	Clean Q30 bases rate (%)	Mapping rate (%)
E13–1	19.65	19.29	98.19	91.93	95.17
E13–2	19.68	19.39	98.50	93.43	95.90
E13–3	19.42	19.11	98.42	93.14	95.84
W5–1	18.72	18.37	98.11	93.75	96.48
W5–2	19.89	19.40	97.54	93.25	95.64
W5–3	17.94	17.59	98.08	94.07	96.41
W42–1	15.00	14.64	97.64	94.36	96.19
W42–2	17.48	17.11	97.89	94.71	95.84
W42–3	13.40	13.10	97.78	94.50	95.69

Characterization and Classification of Expressed Transcripts

The putative lncRNAs were identified and a total of 4,127 putative lncRNAs were obtained after discarding transcripts and all the lncRNAs were expressed in at least one biological replicate (FPKM >0). The Illumina RNA-Seq also produced 13,949 mRNAs. To depict the genomic characterizations of these transcripts, we analyzed their exon number, transcript length and the expression level based on comparisons between lncRNAs and mRNAs. The results showed that an average of 2.19 and 8.74 exons existed in lncRNAs and mRNAs, respectively. On the other hand, the number of lncRNAs with 1 or 2 exons was much higher than for the mRNAs, and the exon size in lncRNAs were larger than that of mRNAs (Figure 1 A-B). Furthermore, lncRNAs were shorter in length (Figure 1 C) and were expressed at lower levels compared with mRNAs (Figure 1 D).

Figure 1

Genomic characterizations of lncRNAs and mRNAs. (A) Exon number distribution of lncRNAs and mRNAs. (B) The percent of lncRNAs and mRNAs with different transcript length. (C) The distribution of transcript length of lncRNAs and mRNAs. (D) The expression levels of lncRNAs and mRNAs in 3 stages. (E) The classification of lncRNAs.

We used FEELnc to distinguish intergenic lncRNAs which would not overlap with any genes from the intragenic lncRNAs, and classified the identified lncRNAs into 3 groups based on their genomic location: 4844 intergenic lncRNAs (55.22%), 1564 intronic (17.83%) and 2364 exonic (26.95%) overlapping lncRNAs. According to the orientations, we subdivided the intronic overlapping lncRNAs to 500 sense and 1064 antisense lncRNAs, and the exonic overlapping lncRNAs to 1650 sense and 714 antisense lncRNAs (Figure 1 E).

Expression Profile of Hepatic Transcriptome in Chickens

Of the 13,949 protein coding genes, 9,749 (69.89%) were simultaneously expressed in the 3 stages (E13, W5, and W42). More specifically, a total of 10,003, 10,004, and 10,005 genes were expressed in E13, W5 and W42, respectively. Among them, only 57 genes expressed at E13 were stage-specific (Supplementary Figure S1 A-B). In contrast, there were no stage specific transcripts expressed among the 4127 putative lncRNAs. The expression profiles for the lncRNAs and mRNAs are illustrated in Figure 2, where the 9 samples lncRNAs and mRNAs were separated into 3 clusters based on their developmental stages; samples at E13 were differentiated first, while W5 and W42 were convergent (Figure 2 A-B). According to the PCA, both lncRNAs and mRNAs were distinguished with the first two components explained as having a variance of 83.72% and 96.28%, respectively (Supplementary Figure S1 C-D), which indicated that lncRNAs and mRNAs expressed in stage-specific patterns. Simultaneously, we calculated Pearson correlations between each two samples randomly to verify the expression dynamics of lncRNAs and mRNAs with chicken liver development. Similarly, the lncRNAs and mRNAs were grouped into 3 clusters by stage, and the correlation between two adjacent time points of the lncRNAs was weaker than for the mRNAs (Figure 2 C). Meanwhile, the Shannon entropy (H) value was used to further validate the specificity of gene expression during developmental stages, and the results showed that the temporal specificity of lncRNAs, including intergenic, extron- and intron-overlapping lncRNAs, was higher than that of mRNAs (Figure 2 D).

Figure 2

The expression profile and temporal expression profiles of mRNAs and lncRNAs. The heat-map showing the expression profile of mRNAs (A) and lncRNAs (B). The top panel is the tree constructed by Pearson correlation. The value represents the log2 transformed values of (FPKM+1). (C) Correlation between every two samples calculated by log2 (FPKM+1) each sample. The upper panel shows mRNAs and the lower panel shows lncRNAs and values represent the pairwise Pearson correlation. (D) The distributions of Shannon entropy-based temporal specificity scores that were calculated for distinct classes of lncRNAs and mRNAs.

Differentially Expressed Genes Identification

A total of 539 differentially expressed lncRNAs were identified through pairwise comparisons (Figure 3). Among them, approximately equivalent numbers of differentially expressed lncRNAs were obtained from E13 vs. W5 (228) and E13 vs. W42 (222), and the least amount of differentially expressed lncRNAs was in the comparison between W5 and W42 (89). Notably, 7 differentially expressed lncRNAs were common among the three comparisons. We identified 6,370 differentially expressed mRNAs (DEGs), and most DEGs were identified from E13 vs. W5, followed by E13 vs. W42 and W5 vs. W42. Importantly, 110 DEGs were common among the 3 comparisons.

Figure 3

Differentially expressed mRNAs and lncRNAs and functional enrichment. Venn diagram of common differential expression genes (A) and lncRNAs (B) in 3 liver developmental stages. (C) Pathway enrichment of 110 DEGs in 3 stages. (D) PCA analysis of alternative splicing levels in mRNAs.

To gain insight into the similarities and differences of the chicken liver across the 3 crucial stages, we analyzed the functions of DEGs by submitting them to DAVID. Functional enrichment of the 3,185 DEGs between E13 and W5 showed that the most significantly enriched GO-BP categories were annotated with terms that involved cell mitosis such as “DNA replication,” “Mitotic nuclear division,” “Cell division,” “DNA replication initiation,” and “DNA duplex unwinding”. The KEGG pathways were enriched with terms such as “DNA replication,” “Carbon metabolism,” “Fatty acid metabolism,” “PPAR signaling pathway,” “Glutathione metabolism,” and “Biosynthesis of unsaturated fatty acids” (Supplementary Table S2). While comparing E13 vs. W42, except for terms involved in cell division and development, GO-BP categories were related to damage repair (DEGs were enriched in “DNA repair,” “Cellular response to DNA damage stimulus,” “Mitotic G2 DNA damage checkpoint,” “Positive regulation of telomerase RNA localization to Cajal body”, “Double-strand break repair via non-homologous end joining”) and immune response (“T cell homeostasis,” “Response to UV”). KEGG pathway annotations were enriched with “Cell cycle,” “Mismatch repair,” “Oxidative phosphorylation,” ect ( Supplementary Table S3). In the comparison of W5 vs. W42, “ATP metabolic process,” “ATP synthesis coupled proton transport,” “Positive regulation of protein secretion,” and “Positive regulation of peptide hormone secretion” were significantly enriched in GO-BP (Supplementary Table S4). These results indicated that chickens at the embryonic stage (E13) were primarily undergoing cell differentiation and that basal metabolisms were active in the developmental stage (W5), while at the age of 42 wk, chickens were highly reproductive and appeared to have signs of aging, with their repair and immune mechanisms activated.

In addition, we analyzed the pathway enrichment of 110 DEGs for the 3 stages, and “Oxidative phosphorylation,” “Mitochondrial respiratory chain complex I,” and “Regulation of T cell differentiation” were significantly enriched in these DEGs (Figure 3 C).

Alternative Splicing Expression of mRNAs

To evaluate the alternative splicing expression of mRNAs, we further estimated the levels of alternative splicing mRNAs in the 3 developmental stages by calculating the PSI value of exons. We observed a stage-dominated clustering pattern of alternative splicing in PC1 (Figure 3 D), which was consistent with the results of the mRNA expression profile, although the power of alternative splicing levels were weaker than the mRNA.

Coexpression of lncRNAs-mRNAs in Hepatic Transcriptome

We used WGCNA to predict the function of lncRNAs by exploring the functional correlation based on the coexpression network of mRNAs and lncRNAs. A total of 10,975 mRNAs and 2863 lncRNAs were used to construct the coexpression network. As a result, 16 modules were obtained, with the most coexpressed transcripts being 5,085, followed by 2,395 and 1,866. We used the overlapped mRNAs (7,555) of the 3 largest modules (68.83%) to perform GO enrichment and KEGG pathway analysis. The results showed that pathways including “Glycine, serine and threonine metabolism,” “Iron ion homeostasis,” “Mitochondrial respiratory chain complex I,” “Oxidative phosphorylation,” and “Mitochondrial inner membrane” were significantly enriched (Supplementary Table S5). This indicated that lncRNAs in chicken liver play an important role in providing energy for organisms through the mitochondrial respiratory chain in chicken, while also, playing a vital role in antioxidant stress.

Five differentially expressed genes (FGB, APOV1, VTG2, FGG, SULT) were selected to confirm the expression profile by qPCR in 3 hepatic developmental stages. The relative expression of each sample in qPCR was compared with the transformed log2 (FPKM+1) values of RNA-sequencing. With the exception of VTG2, which was only expressed in W5 and W42, the other 4 genes showed the same pattern of expression both in RNA Seq and qPCR (Figure 4). The consistent expression of all the DEGs indicated the reliability of our RNA-Seq data.

Figure 4

Validation of DEGs by qPCR. The x -axis indicates gene names; bars on left represents the log2 (FPKM+1) values of RNA-Seq and marked by solid histogram in the y -axis, bars on the right represents the relative expression of qPCR and marked by oblique histogram; and different developmental stages were distinguished by 3 different colors.

DISCUSSION

The liver, as an important energy homeostasis and lipid metabolism organ, has always been a research hotspot in animals. To investigate the gene expression patterns and molecular genetic mechanisms of chicken liver, previous studies have profiled the transcriptomic data of the chicken liver for many aspects. Desert et al. (2018) have revealed that the liver plays a central role in response to short-term fasting in the chicken by using microarray data, including the lipid and acetyl-CoA metabolisms. Cui et al. (2012) have reported that the liver was dominant in fatty acid synthase by investigating the mRNA expression patterns in liver, breast, and thigh (fat related tissues) for different chicken breeds and developmental stages. Monson's group (2016) has established the liver library of the turkey with AFB1 treatment by RNA-Seq. Bourin et al. (2012) have identified two proteases “cathepsin E-A-like/similar to nothepsin” and “uncharacterized protein LOC419301/similar to porin” that were related to a constituent of the egg yolk; and one antiprotease “ovochymase-2/similar to oviductin” was identified as being related to vitelline membrane in chicken liver. Willson et al. (2018) have compared the hepatic transcriptome of 3 chicken breeds (meat, layer type, and F1 layer × meat type) and found that the FoxO signaling pathway was an active driver of growth between meat and layer chickens. Muret et al. (2017) have profiled the lncRNA repertoire of chicken liver and adipose tissue, and found the DRCH24 gene and the divergent lncRNA were coexpressed strongly, revealing a new lncRNA might be involved in cholesterol synthesis. Another study has investigated the liver transcriptome at embryonic day 16, 20, and 1 D after hatching, and they identified several key genes associated with antioxidant enzyme activity (Yang et al., 2018).

The liver in the adult has well-defined roles in the maintenance of nutrient homeostasis, bile acid synthesis, the metabolism of xenobiotics, and endogenous hormones, and the detoxification of exogenous compounds (Grijalva and Vakili, 2013). Meanwhile, in the developing embryonic period and newborn stage, the liver acts as a hematopoietic organ with a key function in generating blood cells (Guo et al., 2009). However, the systematic lncRNA repertoire across the chicken liver development, especially at time points across embryonic stages and hens after incubation during development and growth, remains to be studied. In the present study, we used RNA-Seq to profile the extensive expression patterns of mRNAs and lncRNAs in the previously unexplored chicken liver development at day 13 of the embryonic stage (E13) and the postnatal stages (at 5 wk and 42 wk of age), with the aim to compare the regulatory differences of 3 developmental stages on RNA level. We obtained 158.00 Gb of clean data in total, and identified 4,127 putative lncRNAs and 13,949 mRNAs. Although the lncRNAs exhibited distinct features (such as fewer exons, shorter transcripts and lower expression levels) from the mRNAs, it was expected that the genomic characteristics of lncRNAs were conserved across vertebrates (Pauli et al., 2012; Iyer et al., 2015; Liu et al., 2017). According to expression level analysis of the lncRNAs and mRNAs, we found that 57 genes were expressed stage-specifically at E13 and no stage specific lncRNAs was expressed for the 3 stages. The lncRNAs and mRNAs were grouped into 3 clusters by age, both in hierarchical clustering, and the Pearson correlation power of the lncRNAs were weaker than the mRNAs. These observations suggested that lncRNAs and mRNAs expressed in stage-specific patterns and the lncRNAs were more stage-specific than mRNAs in chicken liver. Furthermore, the H value illustrated in revealed that both intergenic lncRNAs and exonic/intronic overlapping lncRNAs are more temporally specified than mRNAs. In addition, the expression of alternative splicing mRNAs showed similar expression patterns but weaker power compared with mRNAs.

Multiple comparisons of lncRNAs and mRNAs were performed to identify differently expressed genes and transcripts. In the group E13 vs. W5, we identified 228 lncRNAs and 1,802 mRNAs were differentially expressed, 222 lncRNAs and 1,035 mRNAs in E13 vs. W42, and 89 lncRNAs and 199 mRNAs in W5 vs. W42. Among these 3 comparisons, 7 lncRNAs and 110 mRNAs were differentially expressed in common. Function enrichment results showed that terms involved in cell mitosis and development were assigned to DEGs between E13 vs. W5/W42, while DNA damage repair and immune response were primary enriched in W42.

It was noted that “PPAR signaling pathway,” “Long-chain fatty acid biosynthetic process,” and “Fatty acid metabolism” was found to be enriched in DEGs between E13 and W5. Peroxisome proliferator-activated receptors (PPARs) are nuclear hormone receptors that are activated by fatty acids and their derivatives. PPARs have three subtypes (α, γ, and β/δ), each of which is encoded in a separate gene and binds fatty acids and eicosanoids (Braissant et al., 1996). PPARα plays a role in the clearance of circulating or cellular lipids via the regulation of gene expression involved in lipid metabolism of the liver and skeletal muscle. PPARγ is required for the positive effect of its ligands on modulating lipid metabolism, but the inhibiting effect on inflammation might be receptor independent (Chawla et al., 2001). Pascual et al. (2005) confirmed that PPARγ has vital roles in adipogenesis and glucose homeostasis, it can promote adipocyte differentiation to enhance blood glucose uptake, and repress the transcriptional activation of inflammatory response. PPARβ/δ is involved in lipid oxidation and cell proliferation (Desvergne and Wahli, 1999; Feige et al., 2006). Correspondingly, genes associated with fatty acid biosynthesis and metabolism including ANGPTL4, ELOVL1, NAAA, ACADL, and the ACSL gene family were found to be differentially expressed between E13 and W5. These findings validate previous reports that the fatty acid synthesis mainly occurred in the liver of chickens. Reddy and Rao (2006) indicated that the liver tissue could produce more than 90% of the total fatty acids. The abdominal and thigh adipose tissue only synthesize small quantities of fatty acids, while the fatty acids synthesized in the liver are stored in adipose tissues in the form of lipids (Cui et al., 2012).

Moreover, “Oxidative phosphorylation” was assigned to a significant pathway in the W5 vs. W42 and E13 vs. W42 comparison groups. This is a highly efficient way of releasing energy, which takes place inside the mitochondria and is the terminal process of cellular respiration in eukaryotes (Cross, 2004). Of the enriched DEGs, ATP5G1, ATP5G3, COX7A2, COX10, SLC26A47, and NDUFS6 were important regulatory genes in this pathway. ATP5G1 and ATP5G3 encode a subunit of the mitochondrial ATP synthase; the mitochondrial ATP synthase catalyzes ATP synthesis, utilizing an electrochemical gradient of protons across the inner membrane during oxidative phosphorylation. COX7A2 and COX10 are members of cytochrome c oxidase, which catalyzes the electron transfer from reduced cytochrome c to oxygen. This component is a heteromeric complex consisting of 3 catalytic subunits encoded by mitochondrial genes and multiple structural subunits encoded by nuclear genes. The mitochondrially encoded subunits function in electron transfer, and the nuclear-encoded subunits may function in the regulation and assembly of the complex. This nuclear gene encodes polypeptide 2 of subunit VIIa and polypeptide 2 is present in both muscle and nonmuscle tissues (Arnaudo et al., 1992; Merante et al., 1997; Schüll et al., 2015). Fagerberg et al. (2014) reported that SLC26A47 restricted expression towards the liver. This gene encodes a member from a large family of mitochondrial transporters. The nuclear-encoded carrier protein is embedded in the inner mitochondrial membrane. This member of the family is thought to be an uncoupling protein that uncouples mitochondrial respiration from ATP synthesis by dissipating the transmembrane proton gradient (Palmieri, 2013). NDUFS6 encodes a subunit of the NADH: ubiquinone oxidoreductase, which is the first enzyme complex in the electron transport chain of mitochondria. This complex functions in the transfer of electrons from NADH to the respiratory chain (Lazarou et al., 2007).

Oxidative phosphorylation is a vital part of metabolism; it produces reactive oxygen species, such as superoxide and hydrogen peroxide, which lead to the propagation of free radicals, damaging cells and contributing to disease and, possibly, aging. The enzymes carrying out this metabolic pathway are also the target of many drugs and poisons that inhibit their activities. Unsurprisingly, “Glutathione metabolism” was found in the E13 vs. W5 and E13 vs. W42 comparisons. DEGs were GSTT1L, GPX1, GSTK1, GSTA4, GPX3, GCLC, GSS, and GCLM. Glutathione (GSH) is an important antioxidant in organisms; it can prevent damage to important cellular components caused by reactive oxygen species (ROS) such as free radicals, peroxides, and lipid peroxides (Emri et al., 1997; Wu et al., 2004). Makar et al. (1994) investigated the antioxidant defenses of several antioxidants including GSH in cultured chick astrocytes and neurons, and in the forebrains of 23-day-old chick embryos and measured the activities of enzymes involved in glutathione metabolism. They concluded that astrocytes are resistant to ROS and may play a protective role in the brain. The antioxidant defense of chick embryos were studied by Yang, and they revealed that liver tissues played more important roles in the defense system than the heart, and identified that GSTA2, GSTA4, MGST1, GPX3, and HAO2 were involved in GSH metabolism (Yang et al., 2018). Therefore, differences between E13 vs. W5 and E13 vs. W42 illustrate that antioxidant defense plays a vital role in chicken embryo development.

In conclusion, we profiled the lncRNA and mRNA expression patterns of chicken liver in different developmental stages, and predicted the biological functions of lncRNAs based on the coexpression network. We discovered important pathways involved in fatty acid synthesis, oxidative phosphorylation and glutathione metabolism, and identified corresponding DEGs from comparisons among the 3 stages. These data may further serve as a basis for studying the underlying molecular mechanisms of liver development in chickens as well as other species.