Long Non-coding RNA and mRNA Profile of Liver Tissue During Four Developmental Stages in the Chicken

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 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(Hussain et al., , 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.

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-Zero TM 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.

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 package 1 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. 1 http://www.r-project.org/

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 tool 2 . 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 RNAsequencing 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.

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 highquality 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).
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 intronoverlapping 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 stagespecific 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, exonoverlapping, and intron-overlapping, showed increased temporal specificity ( Figure 2E).

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.

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.
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 coexpression 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.

DISCUSSION
It is well established that the liver is the major organ of lipid biosynthesis in the chicken 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., 2015Li et al., , 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.

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.