Embryonic transcriptome and proteome analyses on hepatic lipid metabolism in chickens divergently selected for abdominal fat content

Background In avian species, liver is the main site of de novo lipogenesis, and hepatic lipid metabolism relates closely to adipose fat deposition. Using our fat and lean chicken lines of striking differences in abdominal fat content, post-hatch lipid metabolism in both liver and adipose tissues has been studied extensively. However, whether molecular discrepancy for hepatic lipid metabolism exists in chicken embryos remains obscure. Results We performed transcriptome and proteome profiling on chicken livers at five embryonic stages (E7, E12, E14, E17 and E21) between the fat and lean chicken lines. At each stage, 521, 141, 882, 979 and 169 differentially expressed genes were found by the digital gene expression, respectively, which were significantly enriched in the metabolic, PPAR signaling and fatty acid metabolism pathways. Quantitative proteomics analysis found 20 differentially expressed proteins related to lipid metabolism, PPAR signaling, fat digestion and absorption, and oxidative phosphorylation pathways. Combined analysis showed that genes and proteins related to lipid transport (intestinal fatty acid-binding protein, nucleoside diphosphate kinase, and apolipoprotein A-I), lipid clearance (heat shock protein beta-1) and energy metabolism (NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 10 and succinate dehydrogenase flavoprotein subunit) were significantly differentially expressed between the two lines. Conclusions For hepatic lipid metabolism at embryonic stages, molecular differences related to lipid transport, lipid clearance and energy metabolism exist between the fat and lean chicken lines, which might contribute to the striking differences of abdominal fat deposition at post-hatch stages. Electronic supplementary material The online version of this article (10.1186/s12864-018-4776-9) contains supplementary material, which is available to authorized users.


Background
Broiler is the most efficient meat-producing farm animal, and contributes to alleviating the challenge of food security imposed upon the human society [1,2]. For over half a century, commercial broiler has been selected intensively for growth rate and feed efficiency [1]. However, intensive selection on fast growth rate brings along adverse outcomes, such as obesity and related metabolic syndromes [3]. Excessive fat deposition is undesirable, since it degrades meat quality, decreases feed efficiency, and increases production and health cost [4]. Currently, to reduce fat deposition is still a main objective of commercial broiler selection and breeding program [1,5].
Unlike in mammals, chicken lipogenesis is very limited in the adipose tissue [6], and more than 70% of de novo fatty acid synthesis takes place in the liver instead [7]. Fatty acids synthesized in the liver are incorporated into triacylglycerols, and secreted as very low-density lipoprotein (VLDL). After hydrolysis by the lipoprotein lipase (LPL), fatty acids released from VLDL penetrate adipocytes, where they are re-esterified and stored as triglycerides [8]. Thus, accumulation of triacylglycerols in adipocytes is closely related to lipid metabolism in the liver [9].
The Northeast Agricultural University broiler lines (NEAUHLF) originate from a commercial Arbor Acres grandsire line, and are under divergent selection for abdominal fat content based on abdominal fat percentage (AFP) and plasma VLDL concentration since 1996 [10]. NEAUHLF is a unique animal model to study the molecular mechanism of adipose tissue growth and development. In previous studies, using adipose and liver tissues from the fat and lean chickens at various post-hatch stages (1, 4 and 7 weeks of age), we have found a number of differentially expressed genes (DEGs) and proteins (DEPs) related to lipid metabolism by the microarray and proteomics methods, such as peroxisome proliferatoractivated receptor gamma (PPARγ), liver basic fatty acids binding protein (LBFABP), LPL, adipocyte fatty acid-binding protein (AFABP), apolipoprotein A-I (ApoA-I), and long-chain acyl-coenzyme A dehydrogenase (ACADL) [10][11][12][13].
Embryonic stage now occupies nearly one third of the time to market size for broilers, and is vital to the post-hatch performance of broilers. In order to see if molecular differences of hepatic lipid metabolism existed between our two chicken lines at embryonic stages, we performed digital gene expression profiling and quantitative proteomics on the livers of chicken embryos taken from five embryonic stages, embryonic day 7 (E7), E12, E14, E17, and day 1 after hatch (E21). We identified DEGs and DEPs associated with hepatic lipid metabolism at embryonic stages between the two chicken lines, which could help explain the striking differences of post-hatch abdominal fat content.

Transcriptome analysis on embryonic livers
Between our broiler lines divergently selected for AFP, significant differences exist since generation 4 ( Fig. 1). From the fat and lean chicken lines at generation 14, which had 4.5-fold difference of AFP at 7 weeks of age [14], we collected hepatic tissues from embryos at five important developmental stages (E7, E12, E14, E17 and E21). Total RNAs isolated from hepatic tissues were submitted for sequencing by the digital gene expression profiling technology. For each of the 10 sequenced libraries, an average of 8500 genes (54.8% of all annotated protein-coding genes) was found, and a large number of novel transcripts were also detected (from 36,569 to 45,021) (Additional file 1).
To identify genes differentially expressed (false discovery rate ≤ 0.001 and fold changes ≥2) and potentially involved in hepatic lipid metabolism, we compared the gene expression profiles between the fat and lean chicken lines. We found 521, 141, 882, 979 and 169 DEGs (Additional files 2, 3, 4, 5 and 6) for the five embryonic stages, respectively (Fig. 2a). Furthermore, after analyzing the number of DEGs among the five different stages, we observed that E17 shared the largest number of DEGs with other embryonic stages (Fig. 2b). However, there were only 6 DEGs common to all five developmental stages, which were involved in cell metabolism, and cell apoptosis pathways, respectively ( Fig. 2b and Table 1).
Gene ontology (GO) and KEGG analysis were performed to analyze biological function of DEGs. DEGs at E7 were enriched in the GO term, the ubiquitin-protein ligase activity. No GO term was found at E12. At E14, the highly enriched GO terms were related to biological processes, such as translation, metabolic process, catabolic process, and RNA binding. At E17, GO terms, such as carboxylic acid metabolic process, oxoacid metabolic process, cellular ketone metabolic process, oxidoreductase activity, and catalytic activity, were enriched. The GO term metallopeptidase activity was enriched for E21 (Fig. 3).
Similar to GO analysis, no significant signaling pathway was found for E12 in KEGG analysis (Fig. 4a). While for the remaining 4 embryonic stages (E7, E14, E17, and E21), a common signaling pathway, the metabolic pathway, was found. At E14, pathways mainly related to aminoacyl-tRNA biosynthesis, ribosome, spliceosome and lysosome were enriched. A large number of signaling pathways were enriched for DEGs at E17, including those directly related to lipid metabolism, such as fatty acid metabolism and biosynthesis, PPAR signaling pathway, and glycolysis ( Fig. 4b).
At E21, the amino acid metabolism (tryptophan metabolism and lysine degradation) pathway was enriched.
A total of 26 differentially expressed protein spots were found at the five embryonic stages between the fat and lean lines, and after removing the 4 protein spots that repeated multiple times, 22 differentially expressed protein spots were identified by MALDI-TOF-MS. Except that 2 protein spots (4, 6) could not be identified, the remaining 20 protein spots matched to 17 proteins in the chicken protein database, and to 3 other proteins in protein databases of other species (Table 2).
These 17 proteins matching the chicken protein database included 3 proteins which had been identified  Table 2). For the 14 different DEPs, after GO analysis, 9 of them could be classified into 3 categories with regard to the molecular function: catalytic activity (67%), binding (22%), and transporter activity (11%). However, for the remaining 5, we could not find any hits related to the molcular function. These DEPs were mainly involved in the PPAR signaling pathway, citrate cycle (TCA cycle), fat digestion and absorption, oxidative phosphorylation, aminoacyl-tRNA biosynthesis and MAPK signaling pathways.

Integrated analysis on transcriptome and proteome data
For the 14 DEPs identified by the comparative proteomics, 2 of them, alpha-fetoprotein and OVAL, had no data in the digital gene expression experiment. The transcriptional abundances of the remaining 12 DEPs were validated by qRT-PCR, at the same time points when significant differences of protein abundances were found. However, only were three genes, lysyl-tRNA synthetase (KARS), tyrosyl-tRNA synthetase (YARS) and intestinal fatty acid-binding protein (FABP2), validated to be significantly differentially expressed at E7 and E17, respectively (Fig. 8). When comparing the digital gene expression and qRT-PCR results, KARS and YARS had opposite expression trends at E7, and FABP2 had a similar expression trend at E17. In addition, similar expression patterns to their digital gene expression results were obtained for 5 genes, succinate dehydrogenase [ubiquinone] flavoprotein subunit (SDHA) at E7, thioredoxin domain-containing protein 5 (TXNDC5) at E7, sulfotransferase (SULT) at E7 and E12, inosine triphosphatase (ITPA) at E7, and ApoA-I at E17. By contrast, opposite expression levels were found for 4 other genes, NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 10 (NDUFB10) at E7, heat shock protein beta-1 (HSPB1) at E7 and E12, coronin-1C (CORO1C) at E14, and NDK at E7 and E21. Furthermore, a general correlation analyses showed that digital gene expression and qRT-PCR results were consistent with each other (correlation coefficients r = 0.86 and 0.84 for the lean and fat lines, respectively).
To validate the proteomics results, the western blot experiment was performed for ApoA-I and FABP2. We selected these two proteins, based on the facts that ApoA-I antibody is prepared in the lab and ready-to-use, FABP2 antibody is commercially available and the antigenic epitope of the commercial FABP2 antibody is relatively conserved (72% similarity to human). Western blot results confirmed that ApoA-I and FABP2 had significantly differential abundance at E17 (Fig. 9a, b).
Furthermore, joint analysis on mRNA and protein expression data based on qRT-PCR results and proteomics results was performed, to see if the mRNA and protein abundance levels of 12 DEPs were consistent with each other. The 12 DEPs were classified into two main groups. Group I contained 3 proteins significantly differentially expressed at both transcriptional and protein expression levels, which could be divided into two subgroups. Subgroup 1 included 2 proteins, in which YARS and FABP2 were significantly higher in the fat birds at E7 and E17. In contrast, in subgroup 2, the protein level of KARS was opposite to its transcriptional level at E7. In Group II, 9 DEPs were not significantly differentially expressed at the transcriptional levels, which could also be classified into 2 subgroups. Subgroup 1 has 6 DEPs with well-matched tendency of transcriptional levels, including TXNDC5, SULT, ITPA, NDUFB10, NDK (protein spot 58) at E7, HSPB1 at E12. Subgroup 2 has 6 DEPs with opposite transcriptional trends, including SDHA, HSPB1 at E7, SULT at E12,  (Table 3).

Discussion
Nowadays, to reduce fat deposition is still an important goal for commercial chicken breeding program [5]. In avian species, the liver is the main site of de novo fatty acid synthesis [7], while the adipose tissue serves mainly as a storage tissue [6,15]. Accumulation of triacylglycerols in adipose tissue is directly related with hepatic lipogenesis [9]. Using NEAUHLF, a suitable animal model to study the molecular mechanism of adipose tissue growth and development, we previously found genes and molecular pathways important for hepatic lipid metabolism (PPARγ, LBFABP, ApoA-I, AFABP and glycolmetabolism) in the liver and adipose tissues at 1, 4 and 7 weeks of age by microarray and proteomics analyses [10][11][12][13]. However, whether there are differences on hepatic lipid metabolism between the two lines in embryonic stages remains unknown. In the current study, we compared transcriptome and proteome profiling on hepatic tissues sampled from 5 different embryonic stages between the two lines. Integrated mRNA and protein expression data showed that 8 DEPs (2 and 6 in subgroups 1 of Group I and II, respectively), YARS, TXNDC5, SULT, ITPA, NDUFB10, NDK (protein spot 58) at E7, HSPB1 at E12, and FABP2 at E17, had similar trend of   Proteins related to lipid transport, FABP2, NDK and ApoA-I, were found to be differentially expressed between the two lines in the present study. FABP2, also known as intestinal FABP (I-FABP), had higher protein and mRNA abundances in the fat broilers at stages of rapid growth and development (E12, E14, E17 and E21) found by both transcriptome and proteome analyses,  which were also validated to be significantly differentially expressed between the two lines at E17 by western blot and qRT-PCR. FABP2 is involved in lipid metabolism, especially in the uptake, intracellular metabolism and transport of long chain fatty acids [17,18]. FABP2 may also influence mitochondrial fatty acid oxidation and free cholesterol transport by regulating gene expression and interaction with nuclear receptors [19]. The higher expression level of FABP2 in the fat birds may imply that the fat birds have stronger capacity of lipid transport.
In the present study, mass spectrometry results showed NDK had two isoforms (protein spots 40 and 58). The abundance of protein spot 40 in the fat line was significantly higher than that in the lean line at E7 and E21, and the transcriptional level of NDK was higher in the lean line compared to the fat line at E7 and E21. However, for protein spot 58 at E7, its mRNA and protein abundance was down-regulated in the fat birds. Human studies also indicate that NDK have two isoforms. Though the two isoforms are closely related in amino acid sequences (88% identity), they display significant differences in cellular functions [20]. NDK catalyzes phosphoryl transfer from a nucleoside triphosphate to a nucleoside diphosphate, and functions in the metabolic pathway [21]. NDK regulates synaptic vesicle internalization, where the dynamin GTPase is required to function [22], and is indispensable in energy metabolism in development [23]. Recently, it was found that NDK is a lipid-dependent mitochondrial switch in both phosphortransfer and inter-membrane cardiolipin transfer, which relates to apoptotic signaling and other putative functions, potentially important in lipid metabolism [24,25].
In addition, for ApoA-I that is also involved in lipid transport, we found that the lean birds had a significantly higher ApoA-I protein abundance at E17. However, no transcriptional difference for ApoA-I in the embryonic liver was found between the fat and lean broilers at E17, suggesting the ApoA-I may be regulated at the translational level. We previously examined ApoA-I and its association with fat deposition using genetics, gene expression and proteomics methods in adipose tissues [12,13]. ApoA-I is a major component of high-density lipoprotein (HDL) in the plasma [26], and can promote cholesterol efflux from peripheral tissues to the liver to keep body cholesterol in balance [27]. ApoA-I can not only function as a key lipoprotein to transport cholesterol, but can also inhibit fatty acid synthesis in mice [28]. Thus, we speculate that ApoA-I can influence embryonic liver lipid metabolism, and contribute to the striking differences of abdominal fat deposition between the fat and lean chicken lines.
HSPB1 could be related to lipid clearance, and protein levels were significantly higher in the lean birds at E7 and E12 as revealed by quantitative proteomics. At E12, Fig. 9 Western blots for ApoA-I and FABP2. Four hepatic tissue samples at E17 were assayed for the lean and fat chicken lines, respectively. For ApoA-I, the lean line were significantly higher than the fat line (a, b). In contrast, for FABP2, the fat line had strikingly higher protein level than the lean line (a, b) both transcriptional and protein expression levels of HSPB1 agreed with each other. But at E7, the protein level of HSPB1 was opposite to its transcriptional level, which may be due to post-translational modification [29]. HSPB1 (also known as HSP27) serves as an ATPindependent chaperone [30]. In mammals, it has been reported that obese subjects had higher anti-HSP27 antibody levels [31], and induction of HSP27 may blunt the adverse effect of fat overexposure on insulin function [32]. In diabetic hearts, the phosphorylation of HSP27 enhanced lipoprotein lipase activity and promoted hydrolysis of triglyceride-rich lipoproteins to fatty acids [33]. Phosphorylated HSP27 promotes autophagy and hepatic lipid clearance via autophagy-lysosome pathway in human hepatic cells [34]. We previously found that HSP27 protein was down-regulated in the abdominal adipose tissue of fat birds [12,13], and here we found that HSP27 was down-regulated in the liver tissue of fat We identified also 2 other proteins, SULT and TXNDC5, potentially important for lipid metabolism in chicken embryos. The protein and transcriptional abundance of SULT was all higher in the lean line at E7, whereas the protein level of SULT was opposite to its transcriptional level at E12. In humans, SULT shows nuclear translocation and can be post-translationally modified [35]. The different expression patterns of SULT in protein and transcription levels could also be due to post-translational modification at E12 [35]. Sulfotransferase is a transferase enzyme known to catalyze the transfer of a sulfo group from a donor molecule to an acceptor group of numerous substrates [36], and reactive groups for a sulfonation via sulfotransferases may be part of a protein, lipid, or steroid [37]. It is capable of responding to inflammatory cues and controlling lipid metabolism by PPARγ in humans [38]. The different expression of SULT between the two lines at E7 and E12 showed that SULT may participate in hepatic lipid metabolism via PPARγ in the chicken. Both transcriptional and protein expression levels of TXNDC5 was higher in the fat chickens at E7. TXNDC5 has a protein disulphide isomerase-like domain, and belongs to the thioredoxin family, which is thought to catalyze disulphide formation to aid protein folding or to regulate protein function against endoplasmic reticulum stress induced by oxidative insults. TXNDC5 protein and mRNA levels were significantly associated with hepatic fat content in ApoEknockout mice [39]. TXNDC5 modulated adiponectin signalling by interacting with adiponectin receptor 1 (AdipoR1) [40] and contributed to increased risk of hepatocellular carcinoma development [41]. Moreover, we found also that proteins in the oxidative stress pathway were differentially expressed. Oxidative stress can alter the expression of ApoB, and VLDL secretion [42]. Therefore, TXNDC5 could couple with the control of ApoB levels by the oxidative stress pathway, to exert its effect on subsequent hepatic lipid metabolism.
Succinate dehydrogenase (SDH) and NDUFB10 could participate in energy metabolism. The protein abundance of SDHA was higher in the fat line, whereas the protein level was opposite to its transcriptional level at E7. The mRNA and protein levels of NDUFB10 were higher in the lean lines at E7. SDH is known as succinate-ubiquinone oxidoreductase, a complex of the mitochondrial respiratory chain. The complex is composed of four nuclear-encoded subunits (SDHA, SDHB, SDHC and SDHD). It has a role in citric acid cycle and mitochondrial energy generation in mammals [43]. NDUFB10 is a subunit of the NADH dehydrogenase (ubiquinone) complex, located in the mitochondrial inner membrane [44]. The different expression of SDHA and NDUFB10 indicated that the fat birds and lean birds had differences in regard to energy metabolism.
We found two Aminoacyl-tRNA synthetases (KARS and YARS), which are important for amino acid synthesis, were differentially expressed between the fat and lean birds in liver tissue at E7. Both the protein and mRNA expression levels of YARS were significantly higher in the fat birds at E7. The protein level of KARS was significantly higher in the fat birds at E7 and the mRNA level was significantly higher in the lean birds at E7. They catalyze the aminoacylation of tRNA by their cognate amino acid [45]. A number of studies reported that the functions of ARSs were also associated with RNA splicing, immune responses, angiogenesis and cell fate determination besides protein synthesis [46]. But the functions of the two proteins in chicken hepatic lipid metabolism are not very clear.
Two proteins (ITPA, CORO1C) were found to be differentially expressed between the two lines at E17 and E14, respectively. ITPA had higher protein and mRNA abundances in the lean broilers at E7. The protein abundance of CORO1C was significantly higher in the lean broilers at E14, whereas transcriptional abundance was higher in the fat broilers. However, their functions on lipid metabolism remain to be investigated.
In poultry, fat deposition depends on the availability of plasma triglycerides, which are transported as components of lipoproteins [9]. Fattening therefore is in connection with three aspects of lipid metabolism: (1) lipid synthesis; (2) lipid transport; and (3) lipid utilization. It is reported previously that the liver of the avian embryo has the capacity for lipoprotein synthesis, secretion and β-oxidation [47], though lipogenesis within the embryonic liver tissue is low [48]. Moreover, the chick embryo liver has a very high capacity for β-oxidation, and fatty acid oxidation provides most of the energy that is required for embryo development [49]. As mentioned above, we found that genes/proteins related to lipid transport and energy metabolism were differentially expressed in the embryonic liver between the fat and lean lines. As a result, differences of transport and utilization of lipids as well, will appear. This could be one of the underlying reasons for the significant difference of fat deposition between our two chicken lines, starting at 7 days posthatch [50].

Conclusions
Molecular differences related to lipid transport, lipid clearance and energy metabolism exist for hepatic lipid metabolism at embryonic stages between the fat and lean chicken lines, which might contribute to the striking differences of abdominal fat deposition at post-hatch stages.

Chicken embryos
The fertilized chicken eggs were chosen from the 14th generation of Northeast Agricultural University broiler lines divergently selected for high and low abdominal fat content (NEAUHLF), 200 each for the fat and lean lines, respectively. The fertilized chicken eggs were all hatched in the same conditions at Northeast Agricultural University hatchery.

Collection of liver tissues
Liver tissues were collected from chicken embryos at the five embryonic stages, E7, E12, E14, E17 and E21. The start of the incubation period was referred to as "E1" (1day-old embryos) and "E21", for newly hatched birds. Liver tissues were collected from chicken embryos under aseptic conditions, frozen immediately in liquid nitrogen, and stored at − 80°C.

RNA sample preparation and digital gene expression
Except 30 liver samples from embryos at E7 were pooled together for RNA extraction, for the remaining four embryonic stages, 15 samples were used, respectively. Total RNAs from the 10 samples were extracted using Trizol (Invitrogen), and RNA quality and concentration were evaluated, to ensure RNA integrity number (RIN) ≥ 9.
Sequencing libraries for digital gene expression analyses were prepared according to the following procedures. Messenger RNAs were purified from 6 μg total RNA, and reversely transcribed into cDNA, which was then digested using the restriction enzyme NlaIII (recognition sites: 5′-…CATG/GTAC…-3′). Adapter 1 (Illumina) was then ligated to the 5′-end, and the restriction enzyme MmeI was then added, which can recognize the nucleotide sequences composed of adapter 1 and CATG, and then cut at the downstream 17 bp site, to produce tags labelled with adapter 1. Then adapter 2 (Illumina) was ligated with these tags, and 15 cycles of PCR were performed. PCR products were run on 6% TBE PAGE and purified, which were then submitted to Illumina HiSeq™ 2000 for digital gene expression analysis (BGI, Shengzhen).

Bioinformatics analyses of sequencing data
Procedures for analyzing the digital gene expression profiling data were briefly described as follows. Short reads generated were assessed first for their sequence quality and reads of poor quality were discarded. Libraries containing all possible indices were built. Assessment of sequencing quality was performed, such as distribution of distinct tags (Additional file 8), and proportion of clean tags (all > 96%) (Additional file 9). Reference tag library was created by analyzing mRNAs containing the restriction enzyme NlaIII recognition sites (CATG) and 17 bases of the reference gene sequences, which was then used for the alignment of generated short reads, with only one base mismatch allowed. Unambiguous tags mapped to one gene were annotated and novel transcripts were discovered by comparing to known transcripts in the database. Number of genes and its relationship with the sequencing volume was evaluated (2 M required to identify maximum number of genes, and all 10 libraries > 3 M) (Additional file 10). Annotated genes and their expression levels were determined and used for downstream differential expression and gene pathway analyses. All genes expressed at a level centred around 10 TPM (transcripts per million clean tags). Furthermore, antisense transcripts were also found, but without significant differences in numbers between lines (Additional file 11). Differentially expressed genes were identified according to the methods described previously [51], and false discovery rate was used for multiple testing correction. We used a criteria of FDR ≤ 0.001 and fold changes ≥2. GO and pathway enrichment analysis was performed by using the GO database (http://www.geneontology. org/) [52] and the KEGG pathway (http://www.genome.jp/kegg/) [53].

Protein sample preparation and 2-D DIGE
There were three replicates for samples at each stage, and each replicate contained a mixture of liver tissues from more than five different animals, respectively. Total proteins were extracted from liver tissues using Trizol (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol with minor modifications as follows. The liver tissue samples were dissolved in lysis buffer containing 7 M urea, 2 M thiourea, 4% CHAPS, 10 mM Tris (pH 8.5) and 1 × protease inhibitor cocktails (Roche Diagnostics GmbH, Mannheim, Germany). Then, to remove insoluble materials, samples were centrifuged at 25,000 g for 30 min. Total proteins were purified by 2D Clean-Up Kit (GE Healthcare, Chalfont St Giles, UK), and protein concentration was determined by 2D Quant Kit (Amersham Biosciences Corp., Piscataway, NJ). Protein stock solutions were kept at a final concentration of 5 mg/mL. We compared and analyzed 30 samples on a total of 15 gels (two samples on one gel). Details of the experimental design were briefly described as follows. Protein samples labeled with Cy3 or Cy5 fluorescent dyes were loaded on the same gel, together with an internal standard labeled with Cy2. The labeling reaction was carried out using 400 pmol dyes for 50 μg protein, according to the manufacturer's instructions (GE Healthcare). Each sample was labeled three times to minimize the influence of dye and systematic errors. Then, samples labeled with three different dyes were mixed, and an equal volume of 2 × sample buffer (7 M urea, 2 M thiourea, 4% CHAPS, 50 mM DTT, 2% pharmalytes 3-10 NL) were added. Final sample volume was brought to 350 μL, with additional sample dissolved in a rehydration buffer (7 M urea, 2 M thiourea, 2% CHAPS, 2% pharmalytes 3-10 NL, 20 mM DTT). Firstdimension electrophoresis was conducted with the IPG-phor3 isoelectric focusing system (GE Healthcare) using IPG strips (18 cm, pH 3-10 NL), with a total focusing time of 8 kVh at 20°C. Prior to SDS-PAGE, each strip was equilibrated with 15 mL equilibration buffer A (6 M urea, 50 mM Tris-HCl, pH 8.8, 30% glycerol, 2% SDS, 0. 002% Bromophenol blue, 10 mM DTT) on a rocking table for 15 min, followed by a treatment of another 15 min in 15 mL equilibration buffer B (6 M urea, 50 mM Tris-HCl, pH 8.8, 30% glycerol, 2% SDS, 0.002% bromophenol blue, 25 mg/mL iodoacetamide). The strips were then loaded onto the 12.5% acrylamide gels, and gels were run under a constant power at 12°C first with 2 W/strip for 60 min, and then 15 W/strip, until the bromophenol blue reached the bottoms of the gels.

Scanning and image analysis
Gels were scanned by the Typhoon 9400 scanner (GE Healthcare). Cy2, Cy3, and Cy5 images for each gel were taken at 488/520, 532/580 and 633/670 nm excitation/ emission wavelengths, respectively, adjusting the pixel resolution to 100 mm. All gels were scanned at 50 nm resolution, and the intensity was adjusted to ensure that the maximum volume of each image was within 50,000-80,000. Images were cropped to remove areas extraneous to the gel image, using Image Quant V 5.2 (Amersham Biosciences, UK). Image analysis was performed with DeCyder 6.5 (GE Healthcare). The DeCyder BVA module was used to performing comparative cross-gel statistical analysis of all spots, permitting the detection of differentially expressed spots between experimental groups (t-test, P < 0.05). Protein spots with a fold change of at least 1.5 were analyzed.

Spot picking and in-gel digestion
Gels were fixed and stained with Coomassie brilliant blue (CBB). Proteins of interests, as defined by the 2D-DIGE/DeCyder analysis, were excised from the CBBstained gels. Gel pieces were added into 100 mmol/L NH 4 HCO 3 solution buffer in 30% acetonitrile, to decolor for 15 min. The precipitate was collected and washed in 100% acetonitrile and put aside at room temperature for 5 min. After vacuum drying, the gel pieces were incubated with modified trypsin (sequencing grade) at a final concentration of 50 ng/uL at 4°C for 60 min, and then treated in 50 mmol/L NH 4 HCO 3 for 16 h at 37°C. Digested peptide mixtures were extracted twice with 0. 1% TFA in 60% acetonitrile. Then, the extracted solutions were blended, lyophilized and kept at − 20°C for further identification by MS.

MALDI-TOF-TOF MS analysis and database search
For MALDI-TOF-TOF MS analysis, 1 uL of sample was mixed with 1 uL of matrix and loaded onto the MALDI-TOF slides, and the spot number and sample name were recorded. MALDI-TOF-TOF MS analysis was performed on a 4800 Plus MALDI-TOF-TOF™ analyzer (Applied Biosystems, Foster City, CA, USA). The obtained spectra of proteins were submitted for online database searching against NCBInr and Swiss-Prot databases, using MAS-COT program (http://www.matrixscience.com). The following parameters were adopted when searching protein databases: 0.1 Da mass tolerance for peptides and 0.3 Da mass tolerance of TOF-TOF fragments, allowed 50 ppm mass outlier error, reduced Min S/N to 20. Only significant hits were accepted, as defined by the MASCOT probability analysis (P < 0.05). Positive hits with either protein score confidence interval (CI) % or Ion CI% greater than 95 were considered significant. GO analysis was performed by using the GO database (http://www. geneontology.org/) [52].

qRT-PCR
Liver samples were prepared as aforementioned. For E7, at least 15 liver samples were combined into one biological replicate. For E12 and E14, 5 liver samples were mixed together. For E17 and E21, only were 3 liver samples grouped together. Three biological replicates were assayed for the validation of 12 DEGs, while three (E7 and E12), four (E14 and E17), and five (E21) biological replicates were assayed for 12 DEPs, respectively. Primers were designed to span introns and the sequences were listed as in Additional file 12.
Total RNA of liver tissues (50-100 mg) was isolated using Trizol (Invitrogen, Carlsbad, California) according to the manufacturer's recommendations. RNA concentration was determined by spectrophotometry, and RNA quality was checked on 1.5% agarose gel. cDNAs were reverse transcribed from 1 μg of total RNA at 25°C for 5 min, 42°C for 60 min, and 70°C for 15 min, using 0. 5 μL an oligo(dT) primer (Takara, Daliang, China), and 1.0 μL ImProm-II reverse transcriptase (Promega, Madison, WI) in a final volume of 20 μL. TATA box-binding protein (TBP) was used as an internal reference to normalize the expression data.
qRT-PCR reactions were performed using FastStart Universal SYBR Green Master kit (Roche) on QuantStu-dio™ Real-time PCR System following cycling conditions: 1 cycle at 50°C for 2 min and 95°C for 10 min, 40 cycles at 95°C for 15 s and 60°C for 1 min, and a melting curve analysis (95°C for 15 s, 60°C for 1 min and 95°C for 30 s). Relative expression levels were calculated using