Next Article in Journal
Reperfusion Arrhythmias Increase after Superior Cervical Ganglionectomy Due to Conduction Disorders and Changes in Repolarization
Previous Article in Journal
Secretome of Senescent Adipose-Derived Mesenchymal Stem Cells Negatively Regulates Angiogenesis
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 21
Issue 5
10.3390/ijms21051800
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide Detection of Key Genes and Epigenetic Markers for Chicken Fatty Liver

by
Xiaodong Tan
1,2,†,
Ranran Liu
1,2,†,
Siyuan Xing
1,2,3,
Yonghong Zhang
1,
Qinghe Li
1,2,
Maiqing Zheng
1,2,
Guiping Zhao
1,2,* and
Jie Wen
1,2,*
1
Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing 100193, China
2
State Key Laboratory of Animal Nutrition, Beijing 100193, China
3
Animal Breeding and Genomics, Wageningen University & Research, 6708 PB Wageningen, The Netherlands
*
Authors to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Int. J. Mol. Sci. 2020, 21(5), 1800; https://doi.org/10.3390/ijms21051800
Submission received: 13 January 2020 / Revised: 13 February 2020 / Accepted: 19 February 2020 / Published: 5 March 2020
(This article belongs to the Section Molecular Genetics and Genomics)
Browse Figures
Versions Notes

Abstract

:
Chickens are one of the most important sources of meat worldwide, and the occurrence of fatty liver syndrome (FLS) is closely related to production efficiency. However, the potential mechanism of FLS remains poorly understood. An integrated analysis of data from whole-genome bisulfite sequencing and long noncoding RNA (lncRNA) sequencing was conducted. A total of 1177 differentially expressed genes (DEGs) and 1442 differentially methylated genes (DMGs) were found. There were 72% of 83 lipid- and glucose-related genes upregulated; 81% of 150 immune-related genes were downregulated in fatty livers. Part of those genes was within differentially methylated regions (DMRs). Besides, sixty-seven lncRNAs were identified differentially expressed and divided into 13 clusters based on their expression pattern. Some lipid- and glucose-related lncRNAs (e.g., LNC_006756, LNC_012355, and LNC_005024) and immune-related lncRNAs (e.g., LNC_010111, LNC_010862, and LNC_001272) were found through a co-expression network and functional annotation. From the expression and epigenetic profiles, 23 target genes (e.g., HAO1, ABCD3, and BLMH) were found to be hub genes that were regulated by both methylation and lncRNAs. We have provided comprehensive epigenetic and transcriptomic profiles on FLS in chicken, and the identification of key genes and epigenetic markers will expand our understanding of the molecular mechanism of chicken FLS.

1. Introduction

Chickens are one of the most important sources of meat worldwide and comprise 32.63% of all meat consumption, with more than 66.6 billion meat-type chickens produced in the world in 2017 [1]. With the experience of decades of breeding, the growth rate of chicken has been greatly improved. Nevertheless, excessive fat accretion is a crucial problem during the production, which could result in low feed conversion ratio, high cost of chicken production, as well as the excessive pollution to the environment.
For chickens, the liver is the core organ for lipid synthesis [2,3]. Lipid homeostasis is closely dependent on some hepatic metabolic pathways, including lipid absorption, lipid synthesis, β-oxidation, and lipoprotein transport, but the disorder of these pathways could lead to the fatty liver [4,5]. Fatty liver syndrome (FLS) is characterized by increased lipid accumulation in the liver, which is different from fatty liver hemorrhagic syndrome (FLHS) [6]. Both the mortality rate and egg production are quite different between the two syndromes [7,8]. But if chickens with FLS do not receive timely treatment, the birds will develop FLHS.
In the standardized farming of chickens, FLS is an inevitable problem, and it was reported that this incidence reached 4% or even close to 20% [9,10]. The study of FLS is difficult because of its obscure features and irregular occurrence. Usually, nutrients supplement is widely applied for the prevention of FLS, such as choline [11], betaine [12], polyunsaturated fatty acid [13], and conjugated linoleic acid [14]. However, this approach is costly and not precise enough. Therefore, the research on the molecular mechanism of FLS may be a better solution.
The environment can interplay with genetics in such a manner that traits can be modified [15,16]. This form of trait modification is influenced by epigenetic factors that mediate the effect of the environment on genetics. An example of the epigenetic modification is the biochemical modification of DNA by cytosine methylation of CpG dinucleotides. Variants in hepatic DNA methylation have been linked to lipid metabolism and fatty liver [17,18]. Hepatic methylation of PPARγ is higher in non-alcoholic fatty liver disease (NAFLD) subjects and correlated with plasma fasting insulin levels [19]. Site-specific changes in the methylation level of FAS have been found in male mice fed a high-fat diet [20]. Such results suggest that aberrant methylation of genomic DNA is an epigenetic modification that may relate to abnormal transcription in animals with fatty liver.
The importance of long noncoding RNA (lncRNA) has become a focus of hepatic lipid metabolism and fatty liver [21,22]. A large number of lncRNAs have been identified in the livers of NAFLD patients [23]. FLRL8, FLRL3, and FLRL7 have been demonstrated to show an underlying effect on the Peroxisome Proliferators-activated Receptors (PPAR) signaling pathway by interaction with candidate genes related to lipogenesis and lipid transport, such as FABP5, LPL, and FADS2 [22]. There are few investigations in animals with the one exception that Li et al. has identified, a metabolism-related lncRNA (lncLTR), which is regulated by estrogen and associated with plasma triglyceride in hens [24].
Both DNA methylation and lncRNA are closed related to fatty liver, the analysis of which is a valid approach to explore the molecular mechanism. For chickens under the standardized environment, considering that the FLS is inevitable and nutrition supplement is not precise enough, the research based on the molecular mechanism is urgently needed. Therefore, we started this investigation to detect the key genes and epigenetic markers associated with chicken FLS. Notably, all chickens were reared under the same environment. Whole-genome bisulfite sequencing and lncRNA/mRNA sequencing were performed, and the epigenetic and expression profiles were used to provide the new research targets, which could further the understanding of the molecular mechanism of chicken FLS.

2. Results

2.1. The Slaughter Performance and Serum Biochemical Indices of Chickens with Fatty Liver

To assess the occurrence of fatty liver, we examined the extent of fatty degeneration by H&E and Oil Red O staining. Typical liver features are shown in Figure 1A with distinct pathological changes in chickens with fatty liver. Histological analysis indicated that more than 1/3 of hepatocytes had steatosis and massive accumulations of fat droplets within hepatocytes, as well as enlarged hepatocytes and abnormal liver leaflets (Figure 1A(a–c)). No signs of these abnormalities were observed in normal chickens (Figure 1A(d–f)).
The slaughter performance and serum biochemical indices are shown in Figure 1B–E. We tested the serum lipid-related index and found the content of triglyceride (TG), total cholesterol (TC), and low-density lipoprotein (LDL) to be significantly higher in the fatty liver group, while the serum high-density lipoprotein (HDL) content did not differ between the two groups (Figure 1B). For slaughter performance, body weight (BW) and eviscerated weight (EW) were significantly higher in the fatty liver group. The dressed weight (DW) had a similar trend between groups (p = 0.066, Figure 1C). The weight and the relative weight of abdominal fat and liver increased significantly in the fatty liver group of chickens (Figure 1D,E).

2.2. Transcriptome Profiling Analysis of Liver

Hepatic gene expression profiles were assessed by RNA-seq for both the fatty liver and control groups (Figure 2A). We identified 1177 differentially expressed genes (DEGs), including 516 upregulated and 661 downregulated genes that differed between the fatty liver and control groups (Figure 2B, Table S1). Then we focused on the KEGG enrichment analysis with upregulated and downregulated genes, respectively.
For the 516 upregulated genes, 15 pathways were identified. Eight of which were lipid- and glucose-related pathways, including the Tricarboxylic Acid Cycle (TCA) cycle, biosynthesis of unsaturated fatty acids, fatty acid metabolism, the PPAR signaling pathway, and fatty acid elongation (Figure 2C). For the 661 downregulated genes, 18 pathways were predicted (Figure 2D). Most of these pathways were related to immune function, including cell adhesion molecules (CAMs), the toll-like receptor signaling pathway, and phagosome.
With special interests in DEGs involved in lipid- and glucose-related pathways and immune-related pathways. Sixty of 83 DEGs (72%) were upregulated and annotated to lipid and glucose metabolism. One hundred and twenty-one out of 150 DEGs (81%) were downregulated and annotated to immune function (Figure 2E).
RNA-seq results were validated by quantitative real-time PCR (qRT-PCR) analysis of 10 randomly selected DEGs (PGM1, ADI1, RGS1, PARD6G, DLAT, FADS2, CYP8B1, HAO1, PLIN1, SCD). Correlation analysis of expression results between qRT-PCR and RNA-seq showed a high relevance (R2 = 0.9419, Figure 2F), confirming the reliability of the RNA-seq data.

2.3. Integration Analysis of Methylome and Transcriptome

Overt differences in whole-genome DNA methylation were found between chickens with (n = 3) and without (n = 4) fatty liver. In comparison to control chickens, hepatic methylation levels of fatty liver chickens were lower for up- and downstream of the gene body and various functional regions (Figure 3A,B).
We identified a total of 3041 differentially methylated regions (DMRs) and 1442 differentially methylated genes (DMGs) by comparison of fatty liver and control groups (Table S2). Compared to the control group, the number of hypo-methylated (hypo) DMRs was greater in the fatty liver group (Figure 3C), while the methylation levels of DMRs overlapping with various gene regions were lower in the fatty liver group (p < 0.01, Figure 3D).
DMGs in the promoter were identified, including 80 hyper-methylated (hyper) genes and 139 hypo-methylated genes. KEGG enrichment analysis of those DMGs showed that seven pathways were significantly enriched (Figure 3E). Of the seven pathways, two (carbon metabolism and TCA cycle) were enriched by both DEGs and DMGs (Table 1).
DMGs in gene body regions were found, including 496 hypermethylated genes and 889 hypomethylated genes. Thirteen pathways were significantly enriched with those DMGs (Figure 3F). Three of the 13 pathways (calcium signaling pathway, p53 signaling pathway, and AGE-RAGE signaling pathway in diabetic complications) were also enriched by DEGs (Table 1).
A total of 127 genes were identified as both DEGs and DMGs (Table S3). By correlation analysis, it was found to be significantly associated between expression and methylation of 35 genes (p < 0.05), while the correlation of 17 genes between expression and methylation did not reach the statistical significance (p < 0.1) (Figure 3G,H). For DMRs overlapped with lipid- and glucose-related DEGs and immune-related DEGs, 12 were found in 11 lipid- and glucose-related DEGs (e.g., HAO1, PDK3, and ABCD3), and 16 were found in 14 immune-related DEGs (e.g., BLMH, MARCH1, and CD80) (Figure 3I), for example, immune-related gene MARCH1 had hyper DMRs, which were significantly associated with gene expression. The lipid-related gene ARNTL had a similar regulatory relationship, and ABCD3 may also have been affected by DNA methylation.

2.4. Integration Analysis of the LncRNA and the mRNA Profiles

The striking difference in global lncRNA expression profiles was found between the fatty liver and control groups (Figure 4A). Sixty-seven DE lncRNAs were identified (Figure 4B, Table S4), and those lncRNAs were divided into thirteen clusters by Pearson analysis. Potential cis and trans target genes were predicted and 56 differentially expressed cis target genes and 544 trans target genes were detected.
Correlations were sought between thirteen clusters of lncRNAs and DEGs. A set of 1290 co-expressed lncRNA/mRNA pairs and 544 trans target genes were identified (Table S5), from which a co-expressed network was constructed (Figure S1). Within the network, we found LNC_006756, LNC_001439, LNC_010098, LNC_010111, and LNC_001531 to have the highest degree, suggesting they may play a central role in the process of chicken FLS.
The putative biological function of lncRNAs was predicted by Gene Ontology (GO) enrichment analysis with co-expression DEGs by each lncRNA cluster (Figure 4C). We found that DEGs related to cluster I were associated with lipid metabolism. And LNC_006756, LNC_012355, and LNC_005024 were linked to many lipid metabolism-related genes, including FABP1, PLIN1, MMP1, EPHX2, SLC27A4, and HAO1 (r > 0.95) (Figure 4D). Furthermore, DEGs related to cluster IV and X were associated with immune function. LNC_010111, LNC_010862, and LNC_001272 were found to be linked to immune-related genes, such as DOCK2, MARCH1, PTPRC, and CD4 (Figure 4E). These results indicated that abnormal expression of lncRNAs had a major effect on the gene regulation underlying fatty liver development.
Chromosomal co-expressed genes within 300 kbs upstream and downstream of differentially expressed lncRNAs were assessed as potential cis genes. Fifty-nine lncRNAs were found to be cis-acting on 545 target genes, only 56 of which were differentially expressed between the two groups (Table S6). LNC_008671, LNC_009039, and LNC_010494 were found to have over 10 differentially expressed cis target genes.

2.5. Integration Analysis of the Methylation, the LncRNA, and the mRNA Profiles

By integration of DMR and lncRNA profiles, a list of DMRs was identified that overlapped with the genomic positions of 564 lncRNAs (Table S7). The methylation differences of most lncRNAs were between −0.5 and 0.5. Of those, the expression of five lncRNAs was significantly different between the two groups (Figure S2).
DMGs and targets of lncRNA regulated in the cis and trans way were overlapped to explore the candidate genes, which were related to both DNA methylation and lncRNAs. A total of 23 target genes were detected, including lipid- and glucose-related genes and immune-related genes (e.g., HAO1, ABCD3, and BLMH) (Table 2). Most target genes were associated with more than one lncRNA, with methylation differences mainly distributed in the gene body region.

3. Discussion

A high-fat diet (HFD) is a common and valid means by which to induce a fatty liver model. For the chicken line used in this study, only the Jingxing–Huang (JXH) chickens of the first generation (F0) were induced with HFD and offsprings were fed a normal diet. From the F1 to F3 generations, the incidence of FLS in the fatty liver group was twice as high as that in the control group [10]. No obvious genetic differentiation was observed between the two groups. As in mammals, offspring from parents with metabolic syndrome present a higher risk of metabolic abnormality [25,26]. Herein, the key genes and epigenetic markers for fatty liver may contribute to diagnosing FLS, and perhaps to explain the transgenerational effect of chicken FLS.
DEGs between chickens with and without fatty liver were assessed. Upregulated DEGs were mainly enriched in glucose and lipid metabolism pathways. The process of hepatic lipid synthesis was increased in the fatty liver group, with fat deposition rapidly increased [27]. We also found that downregulated DEGs were mainly enriched in immune pathways, such as the toll-like receptor signaling pathway. Previous studies showed female TLR4−/− mice to have increased obesity but to be partially protected against HFD induced insulin resistance, possibly owing to reduced expression of inflammatory genes in the liver. The downregulated expression of TLRs in the liver may also increase obesity and play a crucial role in the occurrence of fatty liver. Gluckman et al. and Bruce et al. have suggested that DNA methylation modification is associated with fatty liver [17,18].
Exploring epigenetic modification with relation to gene expression has provided a new model by which to associate a genomic function to biological phenotypes and metabolism processes [17,18,28]. Previous studies have shown methylation changes to candidate genes to be related to NAFLD [29]. The integrated analysis herein permitted the detection of genome-wide changes in DNA methylation and in gene expression in chicken fatty liver.
DNA methylation in the promoter is negatively correlated with transcription, while DNA methylation of the gene body can affect transcription either positively or negatively. Jung et al. suggested that overexpression of PDK3 promoted elevated levels of glucose aerobic oxidation, which has an important effect on liver disease [30,31]. In this study, this gene had higher expression and intron decreasing methylation in chickens with fatty liver. Reduced expression of ARNTL has been observed in those with obesity [32]. The gene is involved in the process of fat deposition with high levels of circulating fatty acids in the liver [33], which indicates ARNTL has a possible regulatory role in the process of hepatic lipid accumulation. Two bidirectional DMRs within ARNTL were significantly associated with reduced expression levels, which is consistent with previous reports. Elevated TCA cycle flux is often observed with fatty liver [34,35], but an epigenetic association with genes of the TCA cycle has not been reported. TCA cycle is the ultimate pathway of nutrient metabolism; therefore, our results provide a possible mechanism by which nutrients relieve chicken fatty liver syndrome. In this study, many other key genes (e.g., IGF2BP1, ADI1, and HADHA) were also detected, but the function of these genes has not been reported to relate to fatty liver. Further investigation of these results is warranted.
In addition to gene methylation modification, lncRNAs also have an effect on physiological processes by regulation of gene expression and protein function [36]. Over 1000 lncRNAs have been reported to be associated with fatty liver [23]. Our results identified 67 differentially expressed lncRNAs, only two of which were annotated. A study of a NAFLD mouse model induced by an HFD showed over 290 lncRNAs to be differentially expressed [22]. PER2, a regulator of circadian rhythm, was positively regulated by lncRNA FLRL6 [22], which indicated that a regulatory relationship may be important in NAFLD progression because PER2 has an effect on hepatic lipid metabolism through PPARγ [37]. In this study, PER2 was positively regulated by LNC_010240. We, therefore, inferred that LNC_010240 had a similar function to lncRNA FLRL6 in liver lipid metabolism. Although a large number of lncRNAs were detected, their regulatory mechanisms remain largely unknown [38]. Guo et al. classified lncRNAs to four classes and found the first class to respond to HFD and the third class to respond to liver or metabolic disease [39]. We also divided lncRNAs into thirteen clusters based on the expression pattern that contributed to the function. We found the lncRNA of cluster I to be mainly annotated to lipid metabolism pathways, especially LNC_006756, LNC_012355, and LNC_005024, whose target genes are involved in liver lipogenesis and transport. The lncRNAs of cluster IV and X were annotated to immune-related pathways, including LNC_010111, LNC_010862, and LNC_001272. These results provide candidate lncRNAs that are associated with chicken fatty liver. Besides that, it has been reported that lncRNA has a regulation effect on target genes by the cis way [40]. We found that LNC_008671, LNC_009039, and LNC_010494 had over 10 differentially expressed cis target genes, which suggested that these lncRNAs had a pivotal role in the process of fatty liver generation by regulation of gene expression.
Increasing evidence suggests potential crosstalk between DNA methylation and lncRNA interaction networks [41,42]. We integrated the methylation and lncRNA profile to detect common target genes. For target genes of lncRNAs, we found 23 genes (e.g., ABCD3, BLMH, and HAO1) that may be positively or negatively regulated by DNA methylation and may underlie the regulation of chicken fatty liver. These results provided powerful evidence that epigenetic changes play a critical role in the development of chicken fatty liver by regulation of gene expression. ABCD3 is a transporter of very-long-chain fatty acids that cause lipid accumulation in peroxisomes [43]. Overexpression of ABCD3 results in β-oxidation of palmitic acid [44]. The regulation of ABCD3 is an important target of lipid transport and metabolism. Currently, no studies have explored epigenetic modifications of ABCD3. We found methylation of ABCD3 to be associated with transcriptional repression and that the expression of ABCD3 may be regulated by LNC_012355. These results provided a candidate lncRNA and methylation marker for a new regulatory mechanism of abnormal lipid metabolism. HAO1, a liver-specific peroxisomal enzyme with high fatty acid oxidase activity, is targeted to peroxisomes. Abnormal alterations of this gene are connected to hepatic steatosis [45,46,47]. We detected a DMR in the gene body of HAO1 and a candidate LNC_006765. But a clear association of epigenetic modification and transcription of HAO1 requires further investigation. Both ABCD3 and HAO1 are associated with peroxisome, which indicates the peroxisome maybe a central location in the process of cellular lipid homeostasis. BLMH is an aminohydrolase and is related to some liver diseases. It has been reported that BLMH can interact with lipoprotein and play an important role in fatty liver via homocysteine metabolism [48]. Besides, BLMH is also correlated with hepatocellular carcinoma [49]. However, no studies focused on the epigenetic modification of this gene. In this study, we discovered a candidate lncRNA and increased methylation level on the gene body of this gene, which provided a new idea to explore the regulatory mechanism of this gene. In addition, many other genes were also found to be regulated by epigenetic factors and related to liver disease, such as MARCH1 [50,51], LIMD2 [52], SLC39A8 [52], etc. Meanwhile, some genes have not been proved to be included in the process of fatty liver or other liver disease, such as ASPA, CCDC18, etc. Which provided new targets for the epigenetic study of fatty liver.
In conclusion, we have provided a comprehensive epigenetic and transcriptomic profile for chicken FLS. Besides, the targets (e.g., HAO1, ABCD3, and BLMH) and epigenetic markers identified by the integration analysis could contribute to the understanding of the molecular mechanism of FLS, and perhaps indicate a new research focus of chicken FLS.

4. Materials and Methods

4.1. Ethics Statement

All the animal experiments in this study were conducted in accordance with the guidelines for experimental animals established by the Ministry of Science and Technology (Beijing, China). Ethical approval on animal survival was given by the Science Research Department (responsible for animal welfare) of the Institute of Animal Sciences (IAS, Beijing, China), the Chinese Academy of Agricultural Sciences (CAAS, Beijing, China) with the following reference number: IASCAAS-AE-03.

4.2. Animal Model and Environment

One dwarf chicken line of JXH chicken, highly susceptible to fatty liver, was obtained from the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences. These chickens were evaluated, as reported by Zhang [10]. One fatty liver group and one control group were generated from this line. Briefly, the initial JXH chickens (F0 generation) were randomly assigned to two groups and were fed HFD (fatty liver group) and basal diet (control diet), respectively. The offsprings (F1-F3) were produced by male chickens with (fatty liver group) or without (control group) fatty liver. All the offsprings (F1-F3) in two groups were fed a basal diet. For this study, male chickens from the F3 generation of the fatty liver group (n = 70) and the control group (n = 52) (Figure S3) were used. All chickens were fed the same basal diet formulated according to NRC (1994) and NY/T (33-2004), and raised in three-story step cages (one chicken per cage) under the recommended environmental conditions. Feed and water were provided ad libitum during the study.

4.3. Sample Collection

Blood samples were collected before the tissue sample collection. The serum was isolated after incubation for 1 h at 37 °C and stored at −80 °C. In the 42nd week after hatching, all the male chickens were euthanized by carbon dioxide anesthesia and exsanguination by severing the carotid artery after 12-h fasting. The traits (BW, DW, EW, AFW, and LW) were recorded from part of chickens (two out of five chickens at random). The liver was removed and collected. A portion of the liver was snap-frozen and stored at −80 °C for future RNA and methylation analysis, another portion of the liver was fixed in 4% paraformaldehyde for histological analysis.

4.4. Serum Biochemical Analysis and Liver Histology

The concentration of TG, TC, HDL, and LDL was tested with colorimetric reagents (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). The paraformaldehyde-fixed liver was sectioned and stained with H&E or oil red O (Beijing Xuebang Science and Technology Co., Ltd., Beijing, China).

4.5. Evaluation of fatty liver

Fatty liver was assessed by pathological examination of liver H&E and oil red O stained sections [10]. Fatty liver was identified when the extent of fatty degeneration was greater than 33%, as well as excessive deposition of lipid in hepatocyte. The liver phenotype was also assessed.

4.6. Sequencing and Identification of Differentially Expressed LncRNAs and mRNAs

Total RNA was isolated from livers of male chickens with (n = 3, four samples were used to sequencing but one failed) and without (n = 4) fatty liver and assessed by agarose gel electrophoresis (Figure S4). Complementary DNA (cDNA) libraries were prepared using a NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (Ipswich, MA, USA). After removing rRNA, 150 bp paired-end reads were generated with the Illumina HiSeq X Ten platform (Novogene, Beijing, China). High-quality clean reads were obtained by the removal of low-quality reads from raw data with an in-hour script. An index reference genome (Gallus 5.0) was built and clean reads were aligned to the reference genome using HISAT2 v2.0.4 (https://ccb.jhu.edu/software/hisat2/index.shtml) [53]. Mapped reads were assembled by StringTie (http://ccb.jhu.edu/software/stringtie/) [54]. The expression normalization was performed by DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) [55], and DEGs were defined as |FC| > 1.5 with p < 0.05. We identified lncRNA as previously described [56] using coding potential analysis with CNCI, CPC, and PFAM [57,58,59]. Differentially expressed lncRNAs were defined with DESeq2 (adjust p-value < 0.05, |FC| > 1.5).
All DEGs were annotated from the Ensembl database (http://asia.ensembl.org/index.html) [60], and a statistic for lipid- and glucose-related DEGs was conducted based on their biological process, as well as immune-related DEGs.

4.7. Quantitative Real-Time PCR

Total RNA from 26 frozen liver tissues (12 from the fatty liver group and 14 from the control group) was obtained using an RNA Isolation Kit (Tiangen, Beijing, China). The mRNA was converted to cDNA with a FastQuant RT Kit (Tiangen) following the manufacturer’s instructions. Primers of selected genes were designed based on chicken coding region sequences from the Ensembl database, which are shown in Table 3. ACTB and RPS6 were selected as reference genes for normalization. The qRT-PCR was conducted in triplicate with the SYBR Premix Ex TaqTM reagent Kit (TAKARA, Kusatsu, Japan) with the QuantStuio 7 Flex Real-Time PCR System (Waltham, MA, USA), with the following program: 95 °C for 3 min, 40 cycles of 95 °C for 3 s, annealing temperature for 34 s. Results were analyzed by the 2−ΔΔCt method [61]. The correlation coefficient (R2) between qRT-PCR and RNA-seq was acquired from Pearson analysis.

4.8. Construction and Analysis of lncRNA-mRNA Network

Pearson correlation analysis of lncRNAs was performed, from which thirteen clusters were classified by hierarchical clustering methods [39], with the Pearson correlation coefficient (PCC) > 0.7 between any two lncRNAs in one cluster. Then the cis and trans target genes of differentially expressed lncRNAs were predicted. For trans target genes, we calculated the PCC and significant p-value for the expression levels of each lncRNA-mRNA pair. A lncRNA-mRNA network was constructed with the trans target genes (p < 0.01, |PCC| > 0.95) using Cytoscape software [62]. Gene Ontology annotation was applied with KOBAS (http://kobas.cbi.pku.edu.cn/) [63] to elucidate the function of each cluster of differentially expressed lncRNAs. For cis target genes, we identified chromosomal co-expressed genes within 300 kbps upstream and downstream of differentially expressed lncRNAs.

4.9. Whole-Genome Bisulfite Sequencing and DMGs Identification

Genomic DNA was obtained from the same liver samples (n = 3 in the fatty liver group and n = 4 in the control group) as used for mRNA sequencing. It was assessed by agarose gel electrophoresis (Figure S5). Library preparation was performed as previously described [64] and sequenced with the Illumina HiSeq X Ten platform (Novogene, Beijing, China). Subsequently, 150 bp paired-end reads were generated. The quality of raw data was assessed by FastQC. Clean data were filtered with specific conditions, as previously reported [65]. Before analysis, the transformed chicken reference genome was bisulfite-converted (C to T and G to A). Clean reads were fully bisulfite-converted (C to T and G to A) and mapped to the converted genome version using Bismark [66].
Methylation level (ML) for each C site was analyzed by the sliding-window approach with the sum of methylated and unmethylated reads counted, where: ML (C) = reads (mC) / (reads (mC) + reads (umC)). Herein the focus was on the methylation level for each CpG site. DMRs were identified using DSS software [67,68,69]. Spatial correlation and biological replicates were both considered in the detection process. DMGs were defined as the gene body and promoter region that overlapped with DMRs.
Pearson correlation analysis was used to assess the overlapped genes of DMGs and DEGs. A p-value of 0.05 was set as the threshold for significant correlation. Genes with a p-value of 0.1 were also taken into consideration.

4.10. KEGG Pathways Analysis

Pathway enrichment analysis was performed with ClusterProfiler [70] to explore the function of DEGs and DMGs, respectively. A p-value of 0.05 was set as the threshold for significant enrichment.

4.11. Statistical Analysis

SPSS 25.0 (SPSS, Chicago, IL, USA) was used for statistical analysis. Data are shown as mean ± standard error. Comparisons were performed by Student’s t-test or Wilcoxon signed-rank test. A p-value < 0.05 (*) and p-value < 0.01 (**) implied a statistically significant difference and highly significant difference, respectively. Graphics were drawn using GraphPad Prism 7 (GraphPad Software, San Diego, CA, USA).

Supplementary Materials

The following are available online at https://www.mdpi.com/1422-0067/21/5/1800/s1.

Author Contributions

X.T. and R.L. conceived the experiments, collected samples, and analyzed the sequencing data, wrote and revised the manuscript. S.X. contributed to collect samples, analyzing the sequencing data and revising the manuscript. Y.Z. participated in experiment management and sample collection. Q.L. and M.Z. participated in animal reproduction work and data explanation. G.Z. and J.W. contributed to conceiving and supervising experiments, as well as making the manuscript revision. All authors reviewed the results and approved the final version of the manuscript.

Funding

The research was funded by the National Key Research and Development Program of China (No. 2018YFD0500401); the Project of State Key Laboratory of Animal Nutrition (2004DA125184G1607); the Modern Agro-industry Technology Research System (CARS-41), and the Agricultural Science and Technology Innovation Program (ASTIP-IAS04; ASTIP-IAS-TS-15). The funding agencies were not involved in the experimental design, analysis and interpretation of the data or writing of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

DEGDifferentially Expressed Gene
DMGDifferentially Methylated Gene
DMRDifferentially Methylated Region
FCFold Change
FLSFatty Liver Syndrome
HFDHigh Fat Diet
LncRNALong Noncoding RNA
NAFLDNon-Alcoholic Fatty Liver Disease
PCCPearson Correlation Coefficient
HDLHigh-Density Lipoprotein
LDLLow-Density lipoprotein
TCTotal Cholesterol
TGTriglyceride

References

  1. Faostat. Available online: http://www.fao.org/faostat/en/#home (accessed on 20 February 2020).
  2. Hermier, D. Lipoprotein metabolism and fattening in poultry. J. Nutr. 1997, 127, 805s–808s. [Google Scholar] [CrossRef] [PubMed]
  3. Saadoun, A.; Leclercq, B. In vivo lipogenesis of genetically lean and fat chickens: Effects of nutritional state and dietary fat. J. Nutr. 1987, 117, 428–435. [Google Scholar] [CrossRef] [PubMed]
  4. Burdge, G.C.; Lillycrop, K.A. Fatty acids and epigenetics. Curr. Opin. Clin. Nutr. Metab. Care 2014, 17, 156–161. [Google Scholar] [CrossRef] [PubMed]
  5. Ferreira, D.M.; Simao, A.L.; Rodrigues, C.M.; Castro, R.E. Revisiting the metabolic syndrome and paving the way for microRNAs in non-alcoholic fatty liver disease. FEBS J. 2014, 281, 2503–2524. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Wolford, J.H.; Polin, D. Lipid accumulation and hemorrhage in livers of laying chickens. A study on fatty liver-hemorrhagic syndrome (FLHS). Poult. Sci. 1972, 51, 1707–1713. [Google Scholar] [CrossRef]
  7. Deacon, L. The fatty liver syndrome—History and early observations. In Proceedings of the 23rd Annual Texas Nutrition Conference, College Park, MD, USA, October 1968; pp. 124–126. [Google Scholar]
  8. Reedy, L.M. Some clinical observations on the fatty liver syndrome (FLS) in laying hens. Proc. Tex. Nutr. Conf. 1968, 23, 1046. [Google Scholar]
  9. Grimes, T.M. Causes of disease in two commercial flocks of laying hens. Aust. Vet. J. 1975, 51, 337–343. [Google Scholar] [CrossRef]
  10. Zhang, Y.; Liu, Z.; Liu, R.; Wang, J.; Zheng, M.; Li, Q.; Cui, H.; Zhao, G.; Wen, J. Alteration of Hepatic Gene Expression along with the Inherited Phenotype of Acquired Fatty Liver in Chicken. Genes 2018, 9, 199. [Google Scholar] [CrossRef] [Green Version]
  11. Sherriff, J.L.; O’Sullivan, T.A.; Properzi, C.; Oddo, J.L.; Adams, L.A. Choline, Its Potential Role in Nonalcoholic Fatty Liver Disease, and the Case for Human and Bacterial Genes. Adv. Nutr. 2016, 7, 5–13. [Google Scholar] [CrossRef] [Green Version]
  12. Veskovic, M.; Mladenovic, D.; Milenkovic, M.; Tosic, J.; Borozan, S.; Gopcevic, K.; Labudovic-Borovic, M.; Dragutinovic, V.; Vucevic, D.; Jorgacevic, B.; et al. Betaine modulates oxidative stress, inflammation, apoptosis, autophagy, and Akt/mTOR signaling in methionine-choline deficiency-induced fatty liver disease. Eur. J. Pharmacol. 2019, 848, 39–48. [Google Scholar] [CrossRef]
  13. Dongiovanni, P.; Valenti, L. A nutrigenomic approach to non-alcoholic fatty liver disease. Int. J. Mol. Sci. 2017, 18, 1534. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Kim, J.H.; Kim, Y.; Kim, Y.J.; Park, Y. Conjugated Linoleic Acid: Potential Health Benefits as a Functional Food Ingredient. Annu. Rev. Food Sci. Technol. 2016, 7, 221–244. [Google Scholar] [CrossRef] [PubMed]
  15. McCarty, R. Cross-fostering: Elucidating the effects of genexenvironment interactions on phenotypic development. Neurosci. Biobehav. Rev. 2017, 73, 219–254. [Google Scholar] [CrossRef] [PubMed]
  16. Huang, T.; Hu, F.B. Gene-environment interactions and obesity: Recent developments and future directions. BMC Med. Genom. 2015, 8. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Gluckman, P.D. Epigenetics and metabolism in 2011: Epigenetics, the life-course and metabolic disease. Nat. Rev. Endocrinol. 2011, 8, 74. [Google Scholar] [CrossRef] [PubMed]
  18. Jiang, M.; Zhang, Y.; Liu, M.; Lan, M.S.; Fei, J.; Fan, W.; Gao, X.; Lu, D. Hypermethylation of hepatic glucokinase and L-type pyruvate kinase promoters in high-fat diet–induced obese rats. Endocrinology 2011, 152, 1284–1289. [Google Scholar] [CrossRef] [Green Version]
  19. Sookoian, S.; Rosselli, M.S.; Gemma, C.; Burgueno, A.L.; Fernandez Gianotti, T.; Castano, G.O.; Pirola, C.J. Epigenetic regulation of insulin resistance in nonalcoholic fatty liver disease: Impact of liver methylation of the peroxisome proliferator-activated receptor gamma coactivator 1alpha promoter. Hepatology 2010, 52, 1992–2000. [Google Scholar] [CrossRef]
  20. Cordero, P.; Gomez-Uriz, A.M.; Campion, J.; Milagro, F.I.; Martinez, J.A. Dietary supplementation with methyl donors reduces fatty liver and modifies the fatty acid synthase DNA methylation profile in rats fed an obesogenic diet. Genes Nutr. 2013, 8, 105–113. [Google Scholar] [CrossRef] [Green Version]
  21. Sulaiman, S.A.; Muhsin, N.I.A.; Jamal, R. Regulatory Non-coding RNAs Network in Non-alcoholic Fatty Liver Disease. Front. Physiol. 2019, 10, 279. [Google Scholar] [CrossRef]
  22. Chen, Y.; Huang, H.; Xu, C.; Yu, C.; Li, Y. Long Non-Coding RNA Profiling in a Non-Alcoholic Fatty Liver Disease Rodent Model: New Insight into Pathogenesis. Int. J. Mol. Sci. 2017, 18, 21. [Google Scholar] [CrossRef] [Green Version]
  23. Sun, C.; Liu, X.; Yi, Z.; Xiao, X.; Yang, M.; Hu, G.; Liu, H.; Liao, L.; Huang, F. Genome-wide analysis of long noncoding RNA expression profiles in patients with non-alcoholic fatty liver disease. IUBMB Life 2015, 67, 847–852. [Google Scholar] [CrossRef] [PubMed]
  24. Li, H.; Gu, Z.; Yang, L.; Tian, Y.; Kang, X.; Liu, X. Transcriptome Profile Analysis Reveals an Estrogen Induced LncRNA Associated with Lipid Metabolism and Carcass Traits in Chickens (Gallus Gallus). Cell. Physiol. Biochem. Int. J. Exp. Cell. Physiol. Biochem. Pharmacol. 2018, 50, 1638–1658. [Google Scholar] [CrossRef] [PubMed]
  25. Wei, Y.; Yang, C.R.; Wei, Y.P.; Zhao, Z.A.; Hou, Y.; Schatten, H.; Sun, Q.Y. Paternally induced transgenerational inheritance of susceptibility to diabetes in mammals. Proc. Natl. Acad. Sci. USA 2014, 111, 1873–1878. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Dhana, K.; Haines, J.; Liu, G.; Zhang, C.; Wang, X.; Field, A.E.; Chavarro, J.E.; Sun, Q. Association between maternal adherence to healthy lifestyle practices and risk of obesity in offspring: Results from two prospective cohort studies of mother-child pairs in the United States. BMJ 2018, 362, k2486. [Google Scholar] [CrossRef] [Green Version]
  27. Liu, Z.; Li, Q.; Liu, R.; Zhao, G.; Zhang, Y.; Zheng, M.; Cui, H.; Li, P.; Cui, X.; Liu, J.; et al. Expression and methylation of microsomal triglyceride transfer protein and acetyl-CoA carboxylase are associated with fatty liver syndrome in chicken. Poult. Sci. 2016, 95, 1387–1395. [Google Scholar] [CrossRef]
  28. Pan, X.; Gong, D.; Gao, F.; Sangild, P.T. Diet-dependent changes in the intestinal DNA methylome after introduction of enteral feeding in preterm pigs. Epigenomics 2018, 10, 395–408. [Google Scholar] [CrossRef]
  29. Hardy, T.; Zeybel, M.; Day, C.P.; Dipper, C.; Masson, S.; McPherson, S.; Henderson, E.; Tiniakos, D.; White, S.; French, J.; et al. Plasma DNA methylation: A potential biomarker for stratification of liver fibrosis in non-alcoholic fatty liver disease. Gut 2017, 66, 1321–1328. [Google Scholar] [CrossRef]
  30. Jung, G.S.; Jeon, J.H.; Choi, Y.K.; Jang, S.Y.; Park, S.Y.; Kim, S.W.; Byun, J.K.; Kim, M.K.; Lee, S.; Shin, E.C.; et al. Pyruvate dehydrogenase kinase regulates hepatitis C virus replication. Sci. Rep. 2016, 6, 30846. [Google Scholar] [CrossRef] [Green Version]
  31. Degenhardt, T.; Saramaki, A.; Malinen, M.; Rieck, M.; Vaisanen, S.; Huotari, A.; Herzig, K.H.; Muller, R.; Carlberg, C. Three members of the human pyruvate dehydrogenase kinase gene family are direct targets of the peroxisome proliferator-activated receptor beta/delta. J. Mol. Biol. 2007, 372, 341–355. [Google Scholar] [CrossRef]
  32. Paschos, G.K.; Ibrahim, S.; Song, W.L.; Kunieda, T.; Grant, G.; Reyes, T.M.; Bradfield, C.A.; Vaughan, C.H.; Eiden, M.; Masoodi, M.; et al. Obesity in mice with adipocyte-specific deletion of clock component Arntl. Nat. Med. 2012, 18, 1768–1777. [Google Scholar] [CrossRef] [Green Version]
  33. Shimba, S.; Ogawa, T.; Hitosugi, S.; Ichihashi, Y.; Nakadaira, Y.; Kobayashi, M.; Tezuka, M.; Kosuge, Y.; Ishige, K.; Ito, Y.; et al. Deficient of a clock gene, brain and muscle Arnt-like protein-1 (BMAL1), induces dyslipidemia and ectopic fat formation. PLoS ONE 2011, 6, e25231. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Satapati, S.; Sunny, N.E.; Kucejova, B.; Fu, X.; He, T.T.; Mendez-Lucas, A.; Shelton, J.M.; Perales, J.C.; Browning, J.D.; Burgess, S.C. Elevated TCA cycle function in the pathology of diet-induced hepatic insulin resistance and fatty liver. J. Lipid Res. 2012, 53, 1080–1092. [Google Scholar] [CrossRef] [Green Version]
  35. White, H.M. The Role of TCA Cycle Anaplerosis in Ketosis and Fatty Liver in Periparturient Dairy Cows. Animals 2015, 5, 793–802. [Google Scholar] [CrossRef] [PubMed]
  36. Zhang, K.; Shi, Z.M.; Chang, Y.N.; Hu, Z.M.; Qi, H.X.; Hong, W. The ways of action of long non-coding RNAs in cytoplasm and nucleus. Gene 2014, 547, 1–9. [Google Scholar] [CrossRef] [PubMed]
  37. Grimaldi, B.; Bellet, M.M.; Katada, S.; Astarita, G.; Hirayama, J.; Amin, R.H.; Granneman, J.G.; Piomelli, D.; Leff, T.; Sassone-Corsi, P. PER2 controls lipid metabolism by direct regulation of PPARgamma. Cell Metab. 2010, 12, 509–520. [Google Scholar] [CrossRef] [Green Version]
  38. Chen, L.L.; Carmichael, G.G. Long noncoding RNAs in mammalian cells: What, where, and why? Wiley Interdiscip. Rev. RNA 2010, 1, 2–21. [Google Scholar] [CrossRef]
  39. Guo, J.; Zhou, Y.; Cheng, Y.; Fang, W.; Hu, G.; Wei, J.; Lin, Y.; Man, Y.; Guo, L.; Sun, M.; et al. Metformin-Induced Changes of the Coding Transcriptome and Non-Coding RNAs in the Livers of Non-Alcoholic Fatty Liver Disease Mice. Cell. Physiol. Biochem. Int. J. Exp. Cell. Physiol. Biochem. Pharmacol. 2018, 45, 1487–1505. [Google Scholar] [CrossRef]
  40. Yan, P.; Luo, S.; Lu, J.Y.; Shen, X. Cis- and trans-acting lncRNAs in pluripotency and reprogramming. Curr. Opin. Genet. Dev. 2017, 46, 170–178. [Google Scholar] [CrossRef]
  41. Wang, L.; Zhao, Y.; Bao, X.; Zhu, X.; Kwok, Y.K.; Sun, K.; Chen, X.; Huang, Y.; Jauch, R.; Esteban, M.A.; et al. LncRNA Dum interacts with Dnmts to regulate Dppa2 expression during myogenic differentiation and muscle regeneration. Cell Res. 2015, 25, 335–350. [Google Scholar] [CrossRef] [Green Version]
  42. O’Leary, V.B.; Ovsepian, S.V.; Carrascosa, L.G.; Buske, F.A.; Radulovic, V.; Niyazi, M.; Moertl, S.; Trau, M.; Atkinson, M.J.; Anastasov, N. PARTICLE, a Triplex-Forming Long ncRNA, Regulates Locus-Specific Methylation in Response to Low-Dose Irradiation. Cell Rep. 2015, 11, 474–485. [Google Scholar] [CrossRef] [Green Version]
  43. Morita, M.; Imanaka, T. Peroxisomal ABC transporters: Structure, function and role in disease. Biochim. Biophys. Acta 2012, 1822, 1387–1396. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Imanaka, T.; Aihara, K.; Takano, T.; Yamashita, A.; Sato, R.; Suzuki, Y.; Yokota, S.; Osumi, T. Characterization of the 70-kDa peroxisomal membrane protein, an ATP binding cassette transporter. J. Biol. Chem. 1999, 274, 11968–11976. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Schwam, H.; Michelson, S.; Randall, W.C.; Sondey, J.M.; Hirschmann, R. Purification and characterization of human liver glycolate oxidase. Molecular weight, subunit, and kinetic properties. Biochemistry 1979, 18, 2828–2833. [Google Scholar] [CrossRef] [PubMed]
  46. Recalcati, S.; Tacchini, L.; Alberghini, A.; Conte, D.; Cairo, G. Oxidative stress-mediated down-regulation of rat hydroxyacid oxidase 1, a liver-specific peroxisomal enzyme. Hepatology 2003, 38, 1159–1166. [Google Scholar] [CrossRef]
  47. De Craemer, D.; Pauwels, M.; Van den Branden, C. Alterations of peroxisomes in steatosis of the human liver: A quantitative study. Hepatology 1995, 22, 744–752. [Google Scholar] [CrossRef]
  48. Suszynska-Zajczyk, J.; Wroblewski, J.; Utyro, O.; Luczak, M.; Marczak, L.; Jakubowski, H. Bleomycin hydrolase and hyperhomocysteinemia modulate the expression of mouse proteins involved in liver homeostasis. Amino Acids 2014, 46, 1471–1480. [Google Scholar] [CrossRef]
  49. Okamura, Y.; Nomoto, S.; Hayashi, M.; Hishida, M.; Nishikawa, Y.; Yamada, S.; Fujii, T.; Sugimoto, H.; Takeda, S.; Kodera, Y.; et al. Identification of the bleomycin hydrolase gene as a methylated tumor suppressor gene in hepatocellular carcinoma using a novel triple-combination array method. Cancer Lett. 2011, 312, 150–157. [Google Scholar] [CrossRef]
  50. Xie, L.; Li, M.; Liu, D.; Wang, X.; Wang, P.; Dai, H.; Yang, W.; Liu, W.; Hu, X.; Zhao, M. Secalonic Acid-F, a Novel Mycotoxin, Represses the Progression of Hepatocellular Carcinoma via MARCH1 Regulation of the PI3K/AKT/beta-catenin Signaling Pathway. Molecules 2019, 24, 393. [Google Scholar] [CrossRef] [Green Version]
  51. Bhagwandin, C.; Ashbeck, E.L.; Whalen, M.; Bandola-Simon, J.; Roche, P.A.; Szajman, A.; Truong, S.M.; Wertheim, B.C.; Klimentidis, Y.C.; Ishido, S.; et al. The E3 ubiquitin ligase MARCH1 regulates glucose-tolerance and lipid storage in a sex-specific manner. PLoS ONE 2018, 13, e0204898. [Google Scholar] [CrossRef]
  52. Liu, L.; Geng, X.; Cai, Y.; Copple, B.; Yoshinaga, M.; Shen, J.; Nebert, D.W.; Wang, H.; Liu, Z. Hepatic ZIP8 deficiency is associated with disrupted selenium homeostasis, liver pathology, and tumor formation. Am. J. Physiol. Gastrointest. Liver Physiol. 2018, 315, G569–G579. [Google Scholar] [CrossRef]
  53. Kim, D.; Paggi, J.M.; Park, C.; Bennett, C.; Salzberg, S.L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 2019, 37, 907–915. [Google Scholar] [CrossRef] [PubMed]
  54. Pertea, M.; Kim, D.; Pertea, G.M.; Leek, J.T.; Salzberg, S.L. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 2016, 11, 1650. [Google Scholar] [CrossRef] [PubMed]
  55. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Zhou, J.; Xiong, Q.; Chen, H.; Yang, C.; Fan, Y. Identification of the Spinal Expression Profile of Non-coding RNAs Involved in Neuropathic Pain Following Spared Nerve Injury by Sequence Analysis. Front. Mol. Neurosci. 2017, 10, 91. [Google Scholar] [CrossRef]
  57. Sun, L.; Luo, H.; Bu, D.; Zhao, G.; Yu, K.; Zhang, C.; Liu, Y.; Chen, R.; Zhao, Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013, 41, e166. [Google Scholar] [CrossRef]
  58. Punta, M.; Coggill, P.C.; Eberhardt, R.Y.; Mistry, J.; Tate, J.; Boursnell, C.; Pang, N.; Forslund, K.; Ceric, G.; Clements, J. The Pfam protein families database. Nucleic Acids Res. 2011, 40, D290–D301. [Google Scholar] [CrossRef] [PubMed]
  59. Kong, L.; Zhang, Y.; Ye, Z.-Q.; Liu, X.-Q.; Zhao, S.-Q.; Wei, L.; Gao, G. CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007, 35, W345–W349. [Google Scholar] [CrossRef]
  60. Zerbino, D.R.; Achuthan, P.; Akanni, W.; Amode, M.R.; Barrell, D.; Bhai, J.; Billis, K.; Cummins, C.; Gall, A.; Giron, C.G.; et al. Ensembl 2018. Nucleic Acids Res. 2018, 46, D754–D761. [Google Scholar] [CrossRef]
  61. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  62. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
  63. Xie, C.; Mao, X.; Huang, J.; Ding, Y.; Wu, J.; Dong, S.; Kong, L.; Gao, G.; Li, C.Y.; Wei, L. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011, 39, W316–W322. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Zhang, S.; Qin, C.; Cao, G.; Guo, L.; Feng, C.; Zhang, W. Genome-wide analysis of DNA methylation profiles in a senescence-accelerated mouse prone 8 brain using whole-genome bisulfite sequencing. Bioinformatics 2017, 33, 1591–1595. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Li, C.; Li, Y.; Zhou, G.; Gao, Y.; Ma, S.; Chen, Y.; Song, J.; Wang, X. Whole-genome bisulfite sequencing of goat skins identifies signatures associated with hair cycling. BMC Genom. 2018, 19, 638. [Google Scholar] [CrossRef] [PubMed]
  66. Felix, K.; Andrews, S.R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 2011, 27, 1571–1572. [Google Scholar]
  67. Hao, F.; Conneely, K.N.; Hao, W. A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data. Nucleic Acids Res. 2014, 42, e69. [Google Scholar]
  68. Hao, W.; Tianlei, X.; Hao, F.; Li, C.; Ben, L.; Bing, Y.; Zhaohui, Q.; Peng, J.; Conneely, K.N. Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates. Nucleic Acids Res. 2015, 43, e141. [Google Scholar]
  69. Park, Y.; Wu, H. Differential methylation analysis for BS-seq data under general experimental design. Bioinformatics 2016, 32, 1–10. [Google Scholar] [CrossRef] [Green Version]
  70. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. clusterProfiler: An R package for comparing biological themes among gene clusters. Omics A J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef]
Ijms 21 01800 g001 550
Figure 1. Liver histopathology and phenotypic identification. (A) The plots of liver histopathology. The plots of a–c represent the phenotype, HE staining, and Oil red O staining of fatty liver, respectively, while the plots of d–f represent the same items of normal liver (250 µm). (B) Contents of lipids in serum. Includes high-density lipoprotein (HDL), low-density lipoprotein (LDL), total cholesterol (TC) and triglyceride (TG), n > 10 per group. (C) Bodyweight (BW), dressed weight (DW) and eviscerated weight (EW) of chickens in two groups, n > 20 per group. (D) Abdominal fat weight (AFW) and liver weight (LW) of chickens in two groups, n > 20 per group. E Relative weight of abdominal fat (AFW %) and liver (LW %) of chickens in two groups, n > 20 per group. * represents p < 0.05, ** represents p < 0.01.
Figure 1. Liver histopathology and phenotypic identification. (A) The plots of liver histopathology. The plots of a–c represent the phenotype, HE staining, and Oil red O staining of fatty liver, respectively, while the plots of d–f represent the same items of normal liver (250 µm). (B) Contents of lipids in serum. Includes high-density lipoprotein (HDL), low-density lipoprotein (LDL), total cholesterol (TC) and triglyceride (TG), n > 10 per group. (C) Bodyweight (BW), dressed weight (DW) and eviscerated weight (EW) of chickens in two groups, n > 20 per group. (D) Abdominal fat weight (AFW) and liver weight (LW) of chickens in two groups, n > 20 per group. E Relative weight of abdominal fat (AFW %) and liver (LW %) of chickens in two groups, n > 20 per group. * represents p < 0.05, ** represents p < 0.01.
Ijms 21 01800 g001
Ijms 21 01800 g002 550
Figure 2. Transcriptomic changes and enrichment analysis in normal liver and fatty liver. (A) Heatmap of differentially expressed genes (DEGs) via hierarchical cluster analysis. Different rows correspond to different genes, red and blue strips represent up- and downregulation, respectively. (B) Volcano plot revealing DEGs with various p-values and fold changes (FC). Vertical line, FC = 1.5; horizontal line, p-value = 0.05. Red points, upregulated genes; blue points, downregulated genes. (C) Pathways of enrichment analysis with upregulated DEGs. Rich_factor, the ratio of imputed genes to background genes. (D) Pathways of enrichment analysis with downregulated DEGs. (E) The expression profile of lipid- and glucose-related DEGs and immune-related DEGs. (F) Verification results of qPCR for 10 randomly selected DEGs.
Figure 2. Transcriptomic changes and enrichment analysis in normal liver and fatty liver. (A) Heatmap of differentially expressed genes (DEGs) via hierarchical cluster analysis. Different rows correspond to different genes, red and blue strips represent up- and downregulation, respectively. (B) Volcano plot revealing DEGs with various p-values and fold changes (FC). Vertical line, FC = 1.5; horizontal line, p-value = 0.05. Red points, upregulated genes; blue points, downregulated genes. (C) Pathways of enrichment analysis with upregulated DEGs. Rich_factor, the ratio of imputed genes to background genes. (D) Pathways of enrichment analysis with downregulated DEGs. (E) The expression profile of lipid- and glucose-related DEGs and immune-related DEGs. (F) Verification results of qPCR for 10 randomly selected DEGs.
Ijms 21 01800 g002
Ijms 21 01800 g003 550
Figure 3. Methylation profile and enrichment analysis in normal liver and fatty liver. (A) Distribution of methylation in the gene body, upstream and downstream. Gene body, from transcription start site (TSS) to transcription end site (TES); upstream2k, two thousand base pairs of the upstream region from TSS; downstream2k, two thousand base pairs of the downstream region from TES. (B) Distribution of methylation in various regions, including promoter, 5′UTR, 3′UTR, intron, exon, and repeat region. (C) Statistics for differentially methylated region (DMR) number in the genome-wide range, gene body and promoter region. Hyper means the methylation level of DMR in the fatty liver group is higher than that in the control group, while hypo means the methylation level of DMR in the fatty liver group is lower than that in the control group. (D) Difference level of DMR overlapping with genes. ** represents p-value < 0.01. (E) Pathways of enrichment analysis with DMGs in the promoter. Rich_factor, the ratio of imputed genes to background genes. (F) Pathways of enrichment analysis with DMGs in the gene body. (G) Methylation and expression level of common genes with DMR in the promoter region. (H) Methylation and expression level of common genes with DMR in the gene body region. Type I: differential methylation (DM) > 0, log2FC > 0.585; type II: DM > 0, log2FC < −0.585; type III: DM < 0, log2FC < −0.585; type IV: DM < 0, log2FC > 0.585; type V: DMGs with no significant difference of expression level. (I) Methylation and expression level of lipid- and glucose-related genes and immune-related genes.
Figure 3. Methylation profile and enrichment analysis in normal liver and fatty liver. (A) Distribution of methylation in the gene body, upstream and downstream. Gene body, from transcription start site (TSS) to transcription end site (TES); upstream2k, two thousand base pairs of the upstream region from TSS; downstream2k, two thousand base pairs of the downstream region from TES. (B) Distribution of methylation in various regions, including promoter, 5′UTR, 3′UTR, intron, exon, and repeat region. (C) Statistics for differentially methylated region (DMR) number in the genome-wide range, gene body and promoter region. Hyper means the methylation level of DMR in the fatty liver group is higher than that in the control group, while hypo means the methylation level of DMR in the fatty liver group is lower than that in the control group. (D) Difference level of DMR overlapping with genes. ** represents p-value < 0.01. (E) Pathways of enrichment analysis with DMGs in the promoter. Rich_factor, the ratio of imputed genes to background genes. (F) Pathways of enrichment analysis with DMGs in the gene body. (G) Methylation and expression level of common genes with DMR in the promoter region. (H) Methylation and expression level of common genes with DMR in the gene body region. Type I: differential methylation (DM) > 0, log2FC > 0.585; type II: DM > 0, log2FC < −0.585; type III: DM < 0, log2FC < −0.585; type IV: DM < 0, log2FC > 0.585; type V: DMGs with no significant difference of expression level. (I) Methylation and expression level of lipid- and glucose-related genes and immune-related genes.
Ijms 21 01800 g003
Ijms 21 01800 g004 550
Figure 4. Expression profile of lncRNAs and interaction between lncRNAs and DEGs. (A) Heatmap of differentially expressed lncRNAs via hierarchical cluster analysis. Red and blue strips represent up- and downregulation, respectively. (B) Volcano plot revealing differentially expressed lncRNAs with various adjusted p-values and FC. Vertical line, FC = 1.5; horizontal line, adjusted p-value = 0.05. Red points, upregulated lncRNAs; blue points, downregulated lncRNAs. (C) Gene Ontology (GO) annotation with target DEGs of lncRNAs from each cluster. Boxes with various colors represent 13 clusters of lncRNAs. (D) Co-expression network between DEGs and lipid-related long noncoding RNAs (lncRNAs). The diamonds with various colors represent different clusters of lncRNAs, the circles with light green represent co-expression DEGs; grey line means the positive correlation between lncRNAs and DEGs, while blue line means the negative correlation between lncRNAs and DEGs. (E) Co-expression network between DEGs and immune-related lncRNAs.
Figure 4. Expression profile of lncRNAs and interaction between lncRNAs and DEGs. (A) Heatmap of differentially expressed lncRNAs via hierarchical cluster analysis. Red and blue strips represent up- and downregulation, respectively. (B) Volcano plot revealing differentially expressed lncRNAs with various adjusted p-values and FC. Vertical line, FC = 1.5; horizontal line, adjusted p-value = 0.05. Red points, upregulated lncRNAs; blue points, downregulated lncRNAs. (C) Gene Ontology (GO) annotation with target DEGs of lncRNAs from each cluster. Boxes with various colors represent 13 clusters of lncRNAs. (D) Co-expression network between DEGs and lipid-related long noncoding RNAs (lncRNAs). The diamonds with various colors represent different clusters of lncRNAs, the circles with light green represent co-expression DEGs; grey line means the positive correlation between lncRNAs and DEGs, while blue line means the negative correlation between lncRNAs and DEGs. (E) Co-expression network between DEGs and immune-related lncRNAs.
Ijms 21 01800 g004
Table
Table 1. Pathways enriched by both differentially expressed genes (DEGs) and differentially methylated genes (DMGs).
Table 1. Pathways enriched by both differentially expressed genes (DEGs) and differentially methylated genes (DMGs).
IDPathwayTendency of DEGp-Value 1DMRp-Value 2
gga01200Carbon metabolismup8.53 × 10−4promoter3.69 × 10−2
gga00020Citrate cycle (TCA cycle)up1.70 × 10−2promoter4.82 ×10−2
gga04115p53 signaling pathwayup2.88 ×10−2gene body3.76 × 10−2
gga04020Calcium signaling pathwaydown8.94 × 10−3gene body1.35 × 10−2
gga04933AGE-RAGE signaling pathway in diabetic complicationsdown4.41 × 10−2gene body3.09 × 10−2
1p-value of enrichment analysis with DEGs; 2p-value of enrichment analysis with DMGs.
Table
Table 2. Target genes regulated by both long noncoding RNA (lncRNA) and DNA methylation.
Table 2. Target genes regulated by both long noncoding RNA (lncRNA) and DNA methylation.
lncRNARegulationGeneLog2FCDMRMethylation Difference
LNC_008609,
LNC_008671		transLIMD2−0.83gene body−0.39
LNC_012679transBLMH10.67gene body0.35
LNC_008303transASPA1.41promoter−0.16
LNC_012355trans, cisABCD321.05gene body−0.49
LNC_010111trans, cisCCDC18−0.74gene body0.35
LNC_006756transHAO120.99gene body0.33
LNC_010111, LNC_010862transFLVCR2−1.16gene body−0.4
LNC_010073, LNC_010240transFAM13A0.89gene body−0.18
LNC_002556trans, cisENSGALG00000010639−0.73gene body0.52
LNC_009039trans, cisENSGALG00000011528−0.91gene body0.6
LNC_000820transSLC39A80.73gene body−0.2
LNC_010111transMYO161.44gene body−0.32
LNC_000333transCOTL1−0.84gene body−0.36
LNC_007320, LNC_007320transCELF2−0.72gene body−0.24, −0.31
LNC_005357, LNC_007350, LNC_010111, LNC_010862transRAC21−0.78gene body−0.17
LNC_001439, LNC_001531, LNC_007015, LNC_010098transJAM20.59gene body0.38
LNC_001714, LNC_001742, LNC_006829, LNC_012722transWDPCP1.92gene body−0.24
LNC_001439, LNC_001531, LNC_005357, LNC_007015, LNC_010098, LNC_010111transENSGALG000000339190.84gene body0.11
LNC_001272, LNC_002705, LNC_003079, LNC_007151, LNC_010862, LNC_011070transDOCK21−0.62gene body0.19
LNC_002705, LNC_003079, LNC_008608, LNC_012083, LNC_012722transDIP2C0.92gene body−0.22
LNC_001272, LNC_002705, LNC_003079, LNC_007151, LNC_008608, LNC_010862, LNC_011070transGALNT17−1.32gene body−0.18
ENSGALT00000085791, LNC_001272, LNC_007151, LNC_007350, LNC_010862, LNC_011070transMARCH11−0.78gene body−0.27, −0.33
LNC_001272, LNC_001439, LNC_002556, LNC_005357, LNC_007151, LNC_007350, LNC_008609, LNC_010111, LNC_010494transMEGF113.15gene body−0.17
1 Immune-related genes; 2 Lipid- and glucose-related genes.
Table
Table 3. Specific primers for qPCR.
Table 3. Specific primers for qPCR.
Gene IDGenePrimer SequenceProduct Size (bp)
ENSGALG00000002549RGS1F:5′-AGGATTTACGAGGAGTTTGT-3′105
R:5′-TGTGTGAGTTGGGTCTTG-3′
ENSGALG00000033511CYP8B1F:5′-GGATAAGTGAACAAGACCAGTA-3′132
R:5′-GATACAAGAGGAGCCAGAAG-3′
ENSGALG00000007904DLATF:5′-TTGCTCTCCCTGCTCTGT-3′127
R:5′-CCTATTGTGGCTTTATCTGTCT-3′
ENSGALG00000013594PARD6GF:5′-GCCAACAGCCATAACCTT-3′184
R:5′-CCTCTTCGTCACTCTCCA-3′
ENSGALG00000005739SCDF:5′-GGCTGACAAAGTGGTGATG-3′137
R:5′-GGATGGCTGGAATGAAGA-3′
ENSGALG00000007178FADS2F:5′-CTGAGGAAGACAGCAGAGGACAT-3′153
R:5′-GCAGGCAAGGATTAGAGTTGTG-3′
ENSGALG00000025796ADI1F:5′-ACATGGACGAGTCCCAGGAG-3′113
R:5′-AGCATCCAATCTGCGGTAGG-3′
ENSGALG00000011016PGM1F:5′-ACGGTGAAAACCAAGGCGT-3′103
R:5′-TGAAGTTCTCGGCGTAGTGG-3′
ENSGALG00000023395PLIN1F:5′-GCAATCCAGGGCTTACAG-3′171
R:5′-ATCCAGACGACCAGTTCC-3′
ENSGALG00000008845HAO1F:5′-CGGTTTGTGTTGCTGATTT-3′116
R:5′-TGCTGCTACATTATCTGCTA-3′
ENSGALG00000015082RPS6F:5′-GAGCGCAACGTGAGAACATT-3′92
R:5′-CGGACAACATAGCCCTTCCA-3′
ENSGALG00000009621ACTBF:5′-GAGAAATTGTGCGTGACATCA-3′152
R:5′-CCTGAACCTCTCATTGCCA-3′
F: forward, R: reverse.

Share and Cite

MDPI and ACS Style

Tan, X.; Liu, R.; Xing, S.; Zhang, Y.; Li, Q.; Zheng, M.; Zhao, G.; Wen, J. Genome-Wide Detection of Key Genes and Epigenetic Markers for Chicken Fatty Liver. Int. J. Mol. Sci. 2020, 21, 1800. https://doi.org/10.3390/ijms21051800

AMA Style

Tan X, Liu R, Xing S, Zhang Y, Li Q, Zheng M, Zhao G, Wen J. Genome-Wide Detection of Key Genes and Epigenetic Markers for Chicken Fatty Liver. International Journal of Molecular Sciences. 2020; 21(5):1800. https://doi.org/10.3390/ijms21051800

Chicago/Turabian Style

Tan, Xiaodong, Ranran Liu, Siyuan Xing, Yonghong Zhang, Qinghe Li, Maiqing Zheng, Guiping Zhao, and Jie Wen. 2020. "Genome-Wide Detection of Key Genes and Epigenetic Markers for Chicken Fatty Liver" International Journal of Molecular Sciences 21, no. 5: 1800. https://doi.org/10.3390/ijms21051800

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop