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

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.


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

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). 1C). The weight and the relative weight of abdominal fat and liver increased significantly in the fatty liver group of chickens ( Figure 1D,E).  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.

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,

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 For the 516 upregulated genes, 15 pathways were identified. Eight of which were lipidand 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).

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.

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 coexpressed 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,

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.
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 cisacting 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. 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.  (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.

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.

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.

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

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.

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

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

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 (R 2 ) between qRT-PCR and RNA-seq was acquired from Pearson analysis.

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.

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.

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.

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

Conflicts of Interest:
The authors declare no conflict of interest.