Differential Expression of Genes and DNA Methylation associated with Prenatal Protein Undernutrition by Albumen Removal in an avian model

Previously, long-term effects on body weight and reproductive performance have been demonstrated in the chicken model of prenatal protein undernutrition by albumen removal. Introduction of such persistent alterations in phenotype suggests stable changes in gene expression. Therefore, a genome-wide screening of the hepatic transcriptome by RNA-Seq was performed in adult hens. The albumen-deprived hens were created by partial removal of the albumen from eggs and replacement with saline early during embryonic development. Results were compared to sham-manipulated hens and non-manipulated hens. Grouping of the differentially expressed (DE) genes according to biological functions revealed the involvement of processes such as ‘embryonic and organismal development’ and ‘reproductive system development and function’. Molecular pathways that were altered were ‘amino acid metabolism’, ‘carbohydrate metabolism’ and ‘protein synthesis’. Three key central genes interacting with many DE genes were identified: UBC, NR3C1, and ELAVL1. The DNA methylation of 9 DE genes and 3 key central genes was examined by MeDIP-qPCR. The DNA methylation of a fragment (UBC_3) of the UBC gene was increased in the albumen-deprived hens compared to the non-manipulated hens. In conclusion, these results demonstrated that prenatal protein undernutrition by albumen removal leads to long-term alterations of the hepatic transcriptome in the chicken.

exercise, catch-up growth was observed in the albumen-deprived hens 11 . However when kept in floor pens-a more competitive environment for feed, space and water-the body weight of the albumen-deprived hens remained lower throughout the entire experimental period (55 weeks of age) 7 . Irrespective of the body weight in adulthood, the reproductive performance was markedly diminished by embryonic albumen removal as reflected in the reduced number and weight of the eggs 7,11 . At 10 weeks of age, glucose intolerance was observed in the albumen-deprived hens. This difference however, disappeared in adulthood due to age-related loss in glucose tolerance of the hens 7 .
Induction of an altered phenotype that persists throughout the lifespan implies stable changes in gene expression which would result in altered activities of metabolic pathways 12 . In the low-protein-diet rat model, changes in both hepatic gene expression and DNA methylation have been reported. Six days after weaning, the peroxisome proliferator-activated receptor α (PPARα ) expression is 10.5-fold higher and DNA methylation 20.6% lower, whereas expression of the glucocorticoid receptor (GR) is 200% higher and DNA methylation 22.8% lower 13 . Moreover, these changes persist at least until postnatal day 34 14 . Gong et al. reported an increase in gene expression in the rat of the Insulin-like growth factor 2 (IGF2) in the liver and an increase in DNA methylation of the regulatory region of IGF2 in the liver of male low-protein diet offspring at day 0 15 . Therefore, the objective of the present study was to perform a genome-wide screening for differences in gene expression using RNA-Sequencing (RNA-Seq) in liver samples collected from adult laying hens and differentially expressed genes were grouped according to biological function to discover affected pathways. In addition, it was investigated whether the alterations in gene expression coincided with changes in DNA methylation, in order to examine the possibility of epigenetic mechanisms underlying the observed long-term programming effects.

Results
Physiological results. Detailed physiological results have been published previously 7 . In brief, body weight of the albumen-deprived hens was reduced throughout the entire experimental period (0-55 weeks). In addition, the abdominal fat weight was also reduced in the albumen-deprived hens as compared to the sham-manipulated hens. No differences in absolute or relative liver weight were observed. The reproductive capacity was diminished in the albumen-deprived hens as reflected in the reduced number of eggs and lower egg weight. The plasma triiodothyronine (T 3 ) levels were increased in the albumen-deprived compared to the non-manipulated hens, but not the sham-manipulated hens. An oral glucose tolerance test (OGTT) at 10 weeks of age revealed a decreased glucose tolerance in the albumen-deprived hens. During adulthood, an age-related loss of glucose tolerance was observed in the hens, leading to disappearance of treatment differences in the OGTT.
Genome-wide screening for differentially expressed (DE) genes using RNA-Seq. The logarithms of the fold change and associated P-values of treatment differences are depicted in Supplementary Table S3 online. The heat map of Spearman correlation between all samples using the normalized counts as expression values is shown in Fig. 1. The correlation between all samples was very high (> 0.95) and the biological replicates did not cluster well together. When using P < 0.001 and the log 2 -fold change higher than 1 as cut-off, 156 genes were differentially expressed between the treatments, only 75 of these were previously identified genes (Fig. 2). A heatmap generated based on these differentially expressed genes is depicted in Supplementary Figure 1 and demonstrates a good separation of the three treatments (except for sample sham 3).
Only 3 previously identified genes were differentially expressed between the albumen-deprived hens and both the non-manipulated and the sham-manipulated hens. To proceed with confirmation and validation, an additional 28 previously identified genes where the albumen-deprived hens differed from the sham-manipulated group were included in a list to select genes to validate the RNA-Seq, making a total of 31 genes (Table 1).

Confirmation and validation of DE genes of RNA-Seq via qPCR.
Half of these 31 genes were selected (15 genes) for technical confirmation covering a range of expression levels and biological functions (Fig. 3). Relative expression levels of the albumen-deprived hens versus the non-manipulated and the sham-manipulated hens (n = 3 per group) are displayed as obtained from the RNA-Seq and qPCR results. The fold change estimates by qPCR and by RNA-Seq were strongly correlated (Pearson correlation coefficient was 0.85). The genes of which the expression levels did not fully match are the genes with low expression levels, pointing to decreased sensitivity of the RNA-Seq technique at low expression levels.
To validate the biological significance of the 15 DE genes, sample size was increased to 8 samples per group ( Table 2). As observed from both the RNA-Seq and the qPCR results, 7 genes (46.7%) (TNFSF10, LAPTM4B, TMEM86A, CKS1B, NXPH-2, LRRC3C and BMF) had differentially expression in the liver of the albumen-deprived hens compared to the sham-manipulated hens (P < 0.1) and were thus validated. Eight genes could not be validated. Two genes of these (SEMA6D and H2B-I) had significantly increased expression in the albumen-deprived hens compared to the sham-manipulated hens in the RNA-Seq, but were significantly decreased in the qPCR results. Another 6 genes did not display significant differences in the gene expression measured by qPCR.
Grouping of DE genes according to biological function. The cut-off criteria were loosened (P < 0.005 and Fold change > 1.5) to include more DE genes in our analysis to find key metabolic and biological pathways affected by prenatal protein undernutrition by systems biology analysis using Ingenuity Pathways analysis (IPA) software. Only previously identified genes were included, and all DE genes for which the non-manipulated group differed from the sham-manipulated group were excluded as these do not represent effects of prenatal protein undernutrition. 116 DE genes were obtained, for 13 genes albumen-deprived differed from both non-manipulated and sham-manipulated; for 31 genes albumen-deprived differed from non-manipulated and for 88 genes albumen-deprived differed from sham-manipulated group. Significantly involved biological pathways are listed in Table 3. Two of these networks had key central genes, interacting with many DE genes. The first was a pathway involved in embryonic development, organ development and organ morphology, including 14 DE genes and had ubiquitin C (UBC) as central gene, whereas the second one was involved in cell cycle and carbohydrate metabolism, including 12 DE genes and had glucocorticoid receptor (NR3C1) and Embryonic lethal, abnormal Vision, Drosophila-like 1 (ELAVL1) as central genes (Fig. 4). None of the three central genes (UBC, NR3C1 and ELAVL1) were differentially expressed in the RNA-Seq dataset.
DNA methylation analysis using MeDIP-qPCR. The 9 DE genes of the qPCR (TNFSF10, LAPTM4B, TMEM86A, CKS1B, NXPH-2, LRRC3C, BMF, SEMA6D and H2B-I) were selected for DNA methylation analysis. However, no specific primers for any of the CpG rich fragments of TNFSF10 could be optimised and this gene was therefore excluded from the analysis. In addition, the DNA methylation of the 3 key central regulatory genes (UBC, NR3C1, ELAVL1) identified by the pathway analysis was also examined. For each gene, several CpG rich fragments were examined in the promoter region or around the transcription start site via MeDIP-qPCR ( Table 4). The DNA methylation of most of the examined fragments did not display a treatment effect. Only the DNA methylation of fragment (UBC_3) of the UBC gene was affected by treatment (P = 0.0442). The albumen-deprived hens had significantly more DNA methylation in this fragment than the non-manipulated group, whereas the sham-manipulated group had an intermediate value (Fig. 5).

Discussion
The aim of the present study was to identify differences in gene expression using RNA-Seq and find pathways long-term affected by prenatal protein undernutrition by albumen removal in the liver of the chicken. In addition, it was investigated whether alterations in gene expression coincide with changes in DNA methylation (MeDIP-qPCR).
A first striking observation of the RNA-Seq dataset was the very high correlation between normalized counts as expression levels between all samples (Spearman correlation > 0.95), as demonstrated by the heat map in Fig. 1. This shows that both within as well as between treatments' differences in expression were very small and all samples were very similar. Indeed, posthatch environmental conditions were the same for all groups and therefore differences should only be due to the nutritional programming performed during embryonic development. A possible explanation for the small differences is that the treatments were applied early during embryonic development and measurements were performed in adulthood. The time span between treatment and sampling was nearly 58 weeks. Hence, treatment differences might have been much larger if measurements would have been performed earlier in life. However, at hatch, 2-D DIGE results revealed only 8 differential protein spots between the albumen-deprived hens and either the non-manipulated or the sham-manipulated hens or both 5 . Perhaps, other organs than the liver such as the hypothalamus are the sites of major changes in gene expression causing the  Table 1. List of differentially expressed genes from RNA-Seq. Cut-off criteria were P < 0.001 and log 2 -fold change > 1 between the albumen-deprived hens and the non-manipulated and sham-manipulated hens (3 genes) or between the albumen-deprived hens and the sham-manipulated hens (28 genes). Genes in italic were selected for measurements of gene expression via qPCR.
phenotypic differences. Indeed, the hypothalamus functions as the center of regulation of energy and feed intake and is therefore a good candidate for future research. Using P < 0.001 and log 2 -fold change > 1 as cut-off to determine DE genes in RNA-Seq analysis, 3 genes were DE between the albumen-deprived and both the non-manipulated and sham-manipulated group, 39 between the albumen-deprived and the sham-manipulated and 21 between the albumen-deprived and the non-manipulated group (Fig. 2). To our knowledge, no published studies have examined the effect of low protein diet in mammalian models on the offspring via RNA-Seq. Previously, the effect of low protein diet during gestation on the offspring has been analyzed by microarray in the liver of mouse 16 and rat 17 . In adult rat offspring (day 84), only a small number of genes was affected, 222 genes were upregulated, whereas 89 were downregulated (False discovery rate < 0.05; fold change ≥ 1.5) 17 , although this number was much larger than in the present study. This could point at a distinction between application of strictly nutritional effects and involvement of secondary maternal effects. In addition, in these mammalian models the reduction in protein content in the diet is frequently accompanied by an increase in carbohydrates (e.g. glucose, sucrose or starch) as the diets need to be isocaloric. Therefore, observed effects might just as well be effects of carbohydrate overload instead of protein restriction 18 . Technical arguments to explain the low number of DE genes include lack of annotation of certain genes in the chicken. Indeed, from the present dataset 156 genes were differentially expressed, but only 75 of these were annotated in the databank. Although many of these 'novel genes' are probably mapping artefacts, some of these might represent true new genes. RNA-Seq datasets will become very important in the near future to improve the annotation of the chicken genome and identify more genes.   Expression of 15 DE genes of the RNA-Seq dataset was validated using qPCR. Technically (n = 3), a strong correlation between results from both techniques could be seen. The genes of which the expression levels did not fully match are the genes with low expression levels, pointing to decreased sensitivity of the RNA-Seq technique at low expression levels. Biologically (n = 8), however, only half of the DE genes showed the same results using qPCR. This is in agreement however with the level of consistency observed in other studies 19 . The differences between both techniques may be ascribed to several factors such as interindividual variation as liver samples were collected from different chickens. Many publications of genome-wide expression studies lack validation, which may provide misleading conclusions.
DE genes (P < 0.005 and fold change > 1.5) were grouped to find key metabolic and biological pathways affected by prenatal protein undernutrition using systems biology analysis (Ingenuity Pathways analysis (IPA) software). Seven important physiological system development and functions were affected by the applied  SLC3A1, SLC7A10, TNFSF10, MTTP, LRAT, SLC6A6, DCT, IP6K2,  RACGAP1, SLC20A2, SLC39A14, ABCC9, GRIN2C, SIK1, DIO2,  CD200, NR0B1, ULK1, FOXP2, HPS5   2.91E-05-4.78E-02  20   Small molecule biochemistry  SLC3A1, SLC7A10, ASNS, DIO2, GFPT2, IDO2, SLC6A6, TP53I3,  ABP1, LRAT, MTTP, PDE11A, TNFSF10,    treatment. These pathways include 'embryonic and organismal development' , 'organ morphology and development' , 'tissue development' and 'reproductive system development and function' , as expected by applying a treatment early during embryonic development. These results are in agreement with our previously published data on both the peri-and postnatal development of the albumen-deprived hens [5][6][7]11 . Although no differences were detected in one-day-old chick weight, the proportional carcass weight and the water content of the carcass were increased in the albumen-deprived group 5 . In addition, on embryonic day 20, the plasma thyroxine (T 4 ) concentration was reduced for the albumen-deprived group, indicating a decreased metabolic rate. Body weight and feed intake were reduced during the young to juvenile phase, whereas at adult age the body weight either decreased or increased depending on the posthatch environmental conditions 7,11 . These results are probably the consequence of differences in embryonic development. Still, it is remarkable that this effect was still apparent at adult age. In agreement, maternal protein restriction in the rat offspring affected the 'developmental process' as biological pathway process 17 . Genes involved in the 'reproductive system development and function' pathway were also demonstrated to be affected. Irrespective of posthatch body weight, the reproductive performance of the albumen-deprived hens was seriously diminished as reflected in both a reduced number and weight of the eggs 7,11 . Moreover, these eggs had a different composition (increased proportional yolk and decreased proportional albumen) and were of inferior quality as more second grade eggs were laid by the albumen-deprived group 11 . The reduced reproductive performance is in agreement with the study of Rae et al. linking prenatal undernutrition of ewes to diminished ovulation rate of offspring in adult life 20 .
Additionally, several molecular and cellular functions were affected: 'amino acid metabolism' , 'molecular transport' , 'small molecule biochemistry' , 'cell death and survival' , 'cell-to-cell signaling and interaction' , 'carbohydrate metabolism' and 'protein synthesis' . Lillycrop et al. also found altered molecular functions due to prenatal protein undernutrition in rats: receptor binding, tetrapyrrole binding, and UDP-glycosyltransferase activity, cation and anion transmembrane transporter activity, growth factor activity and ATPase activity 17 . Genes with a wide range of functions were demonstrated to be altered in the present study, which is consistent with our previous studies demonstrating a wide range of both physiological and metabolic alterations. By partial albumen replacement with saline, the embryonic protein availability was decreased and therefore changes in the expression of genes related to the amino acid metabolism could be expected. Moreover, changes in amino acid metabolism were also observed by screening for differential protein abundances in the liver of newly hatched chicks. Affected pathways included valine, methionine, glutamate and cysteine degradation 5 , although the precise genes/proteins affected differed between these two studies. In this respect, changes in protein synthesis were also inferred. A previous study from our group also found indications of an altered protein metabolism in broilers treated by albumen removal before incubation, suggesting a transient increase in muscle proteolysis 10 . Although these results may be interpreted as ordinary and expected, it provides evidence that the applied treatment creates strictly nutritional changes.
Another affected pathway was the 'carbohydrate metabolism' , although few genes were represented here. However, changes in glucose metabolism have been demonstrated before by differential protein abundances in the liver of newly hatched chicks. Affected pathways included the gluconeogenesis, the tricarboxyl acid cycle and pyruvate fermentation to lactate 5 . Furthermore, the albumen-deprived hens exhibited reduced glucose tolerance at 10 weeks of age possibly by a decreased insulin production or increased insulin resistance. However this difference disappeared at adult age due to age-related loss of glucose tolerance in the chicken 7 . Also, in mammalian maternal low protein models, differences in gene expression involved in the carbohydrate metabolism were observed. In rat dams, fed a protein-restricted diet, an increase of the regulatory enzyme of the gluconeogenesis phosphoenolpyruvate carboxykinase PCK1 mRNA and increased activity was detected in liver of the progeny  Table 4. DNA methylation analysis of genes of interest. Analysis was performed for 9 DE genes from RNA-Seq and qPCR and 3 key central genes identified by pathway analysis by MeDIP-qPCR. until 11 months of age, suggesting that programming of the metabolism also extends to the regulation of gene expression 21 .
In offspring from low-protein-diet rat models, alterations in DNA methylation were observed in association with changes in gene expression 13,14 . In the present study, a screening of the DNA methylation of the 9 DE genes and the three key central genes (UBC, NR3C1 and ELAVL1) was performed by MeDIP-qPCR. The data suggest that DNA methylation of fragment UBC_3 of the ubiquitin C (UBC) gene is increased in the albumen-deprived hens compared to the non-manipulated hens. The fact that the gene expression in this region did not differ in the RNA-Seq experiment does not exclude the possibility that these changes in DNA methylation could be related with distal regulation of other genes. The fragment is located around the transcription start site, starting from 31 bp upstream of the 5′ -UTR (Untranslated Region) to 155 bp into this region and 553 bp upstream of the translation start codon (CTG) (Fig. 5). 14 CpG's are contained within this fragment and one or more of these CpGs have increased methylation in the albumen-deprived hens. Ubiquitin C is one of the genes that produces ubiquitin, the precursor of polyubiquitin. Ubiquination is associated with protein degradation, DNA repair, cell cycle regulation, kinase modification, endocytosis, and regulation of other cell signaling pathways. Since the present study only performed a DNA methylation analysis biased towards specific genes, a genome-wide approach will be needed in order to determine the real extent of changes in DNA methylation that could be generated as result of the experimental treatment used. Moreover, other epigenetic modifiers such as histone modification or microRNA could be examined, as demonstrated in a low protein rat model 22,23 .

Conclusions
In conclusion, the present results demonstrate for the first time that prenatal protein undernutrition by albumen removal leads to long-term alterations of the hepatic transcriptome in the chicken. As expected, pathways involved in embryonic development were affected by this treatment. In addition, changes in amino acid metabolism and protein synthesis prove the efficacy of the application of strictly nutritional effects. Limited effects of DNA methylation were observed in the regulation of the currently examined DE genes.

Ethics statement. All experiments were conducted in strict accordance with the European Communities
Council Directive (2010/63/EC) and were approved by the Institutional Ethical Committee of KU Leuven (P002-2012).
Chickens. The set-up of this experiment was previously described by 7 . More information on the technique of albumen removal is available in 11 . In brief, fertilized Isa Brown layer-type eggs (Vepymo, Poppel, Belgium) were randomly divided over the three treatments and incubated together. After one day of incubation, a hole was drilled in the egg, 3 ml of outer thin albumen was removed using a needle and syringe and replaced by about the same volume of sterile saline, followed by sealing of the hole using a drop of paraffin (albumen-deprived group). A sham-manipulated group was mock-treated similar to the albumen-deprived group, except for the actual albumen removal and saline injection. A third group received no treatment (non-manipulated group). A discussion of the technique of albumen removal can be found in 10,11 . After hatch, only the female chicks were reared in floor pens with wood shavings as litter until 55 weeks of age. All floor pens were located in one environmentally controlled room. Room temperature was initially set at 34 °C and was gradually decreased until 20 °C at 5 weeks of age. This temperature was maintained until 55 weeks of age. A 23 h light cycle was initially provided and gradually decreased to 10 h at 6 weeks. At 14 weeks, the light cycle was gradually prolonged to 15 h at 19 weeks to stimulate sexual maturation. The hens received soy-wheat-corn based diets ad libitum formulated according to the developmental requirements based on the Isa Brown Management Guide (Research Diet Services, Wijk bij Duurstede, the Netherlands) 11 . At 55 weeks of age, the laying hens were killed and 8 liver samples per group were randomly collected and snap frozen in liquid nitrogen before storage at − 80 °C. Frozen liver tissue was homogenized by grinding into powder in liquid nitrogen before use.
Genome-wide screening for differential gene expression using RNA-Seq. RNA isolation. Total RNA was extracted from approximately 50 mg of crushed liver collected from hens at 55 weeks of age (n = 8 per group) as previously described by 5 . The RNA concentration and quality (260/280 ratio) was measured using UV-spectroscopy (Implen, Westburg, Leusden, The Netherlands). The integrity of the RNA was electrophoretically verified. In addition, integrity of RNA samples used for RNA-Seq was checked on the BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA).
Sequencing and quality control of the reads. The sequencing was performed in collaboration with the Genomics Core of the UZ Leuven (Belgium). Samples (n = 3 per group, 2 μ g RNA) were prepared by TruSeq library preparation (RS122, Illumina, San Diego, CA, USA) to generate single end unstranded sequencing libraries, according to the manufacturer's recommendations. Sequencing of all samples was carried out in 1 lane of 1 flow cell on the HiSeq2000 (Illumina), using single end chemistry with read lengths of 50 base pairs. Between 17.3 and 26.5 million reads were sequenced for each sample. The quality of the reads was checked using the FASTQC software (Illumina). All the reads passed the Phred score (≥ 28).
Data analysis of the RNA-Seq data. Data analysis was performed in collaboration with the Nucleomics Core of the VIB (Belgium). Reads were mapped to the chicken genome (Galgal4.71) using TopHat (v2.0.8b) 24 . Quality filtering was performed with SAMtools (0.1.19-44428cd) removing all reads from the alignment that are non-primary mappings or have a mapping quality ≤ 20 25 . The number of mapped reads varied between 16 and Scientific RepoRts | 6:20837 | DOI: 10.1038/srep20837 24.8 million reads per sample as calculated by Picard (1.92) (BroadInstitute) resulting in a mapping efficiency of approximately 93%. Per sample, the transcripts were identified from the mapping with Cufflinks v2.1.1 26 . With the Cuffmerge tool, all per-sample transcript lists were merged together with the reference annotation into one file. With Cuffcompare, the identified transcripts were matched with known transcripts and genes, based on overlapping coordinates on the reference genome. About 22.1% novel loci were identified, resulting in a total of 22,082 genes. The number of reads in the alignments that overlap with the gene features were counted using HTSeq-count (0.5.4p3) 27 . The reads that could be attributed to more than one gene (ambiguous, about 1.3 million per sample) or could not be attributed to any gene (no feature, about 1.1 million per sample) were not counted. Exclusion caused by ambiguity excluded about 6.9% and no feature about 6.1% of the mapped reads per sample. The 8,483 genes that had less than 1 counts-per-million were considered 'absent' and therefore removed from the dataset 28 . As such, 13,599 genes were left. Per sample, GC-content was corrected using full quantile normalization on bins of GC-content with the EDASeq package from Bioconductor 29 for within-sample normalization. Between-sample normalization was carried out to correct for library size and RNA composition, as it is known that these are sources of large variation, using full quantile normalization with the EDASeq package from Bioconductor 30 . The statistical modeling assumed the design of the experiment as log (Count) = β 1 × non-manipulated + β 2 × sham-manipulated + β 3 albumen-deprived. For each gene, the coefficients β were estimated with the EdgeR package (3.2.3) 31 of Bioconductor, by fitting a negative binomial generalized linear model (GLM) 28 . The following comparisons were tested: non-manipulated vs. albumen-deprived, sham-manipulated vs. albumen-deprived and sham-manipulated vs. non-manipulated. Differentially expressed genes were selected for confirmation by qPCR based on the criteria that the p-value should be < 0.001 combined with a cut-off fold change of 2 or absolute log 2 -ratio larger than 1 32 .

Identification of relevant biological pathways. A list of differentially expressed genes for pathway
analysis was created by selecting the differentially expressed transcripts with P < 0.005 and fold change > 1.5. This list of differentially expressed genes between the albumen-deprived hens and both non-manipulated and sham-manipulated hens or only the sham-manipulated hens was imported into the Ingenuity Pathway Analysis (IPA; Ingenuity Systems, Redwood City, CA, USA) to identify biological interactions between the gene products. The biological interaction scores were defined by the IPA statistical algorithm based on its knowledge base, and the P value was calculated according to Fisher's exact test.

Confirmation and validation of RNA-Seq results via qPCR. DNase treatment and reverse transcrip-
tion. RNA samples (n = 8 per group) were DNase treated with RQ1 RNase-Free DNase (M6101, Promega, Fitchburg, WI, USA), according to the manufacturer's instructions. The RNA was transcribed into cDNA using the Reverse Transcription system (A3500, Promega, Madison, WI, USA). Denaturation was performed for 3 min at 80 °C followed by 45 min at 42 °C for the reverse transcription. The reaction was terminated by heating the samples for 5 min at 95 °C.
Primer design. The primer design was described previously 5 . In Supplementary Table S1 online, all primers are listed, both for reference genes and the genes of interest. All measurements were performed in the linear range of amplification with correlation coefficient > 0.99. Final primer concentration was 0.3 μ M for all primer pairs, except for BMF, where 0.15 μ M was used.
qRT-PCR. Quantitative real-time PCR (qRT-PCR) measurements were performed as described previously 5 . For the analysis of the qRT-PCR output, the 2 −ΔΔCT method of relative quantification was used 33 . Expression of genes was normalized to the geometric average of the two references genes: β -actin (ACTB) and peptidylprolylisomerase D (PPID).Technical confirmation results were based on the same three samples as used in the RNA-Seq, whereas in the biological validation the sample size was increased to 8 samples per group.

Methylated DNA immunoprecipitation and quantitative real-time PCR (MeDIP-qPCR).
Genomic DNA isolation. One ml of SNET buffer (1% SDS; 0.4M NaCl; 5 mM EDTA; 20 mM Tris.HCl pH 8.0) and 20 μ l of Proteinase K (00120437, Thermo Scientific, Waltham, MA, USA) was added to 100 mg of homogenized liver tissue and incubated overnight in a water bath at 56 °C. The next day 1 ml of a phenol: chloroform: isoamyl alcohol 25:24:1 mixture (A0889, AppliChem, Darmstadt, Germany) was added and thoroughly vortexed. After centrifugation (14,000 rpm, 10 min, 4 °C), the top layer was transferred to a new tube and the process of adding, vortexing and centrifugating was repeated several times with addition of subsequently 1 ml of the phenol: chloroform: isoamyl alcohol 25:24:1 mixture and 1 ml of chloroform (1.02445, Merck Millipore, Readington Township, NJ, USA). Finally, 1 ml of isopropanol (20842, Prolabo, VWR, West Chester, PA, USA) was gently mixed with the top layer, centrifuged (14,000 rpm, 5 min, 4 °C) and the resulting supernatant was discarded. The DNA pellet was washed with 500 μ l 70 ethanol (14,000 rpm, 5 min, 4 °C) and the pellet was dissolved in 150 μ l milliQ water. The DNA concentration and purity was measured on a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA).
Shearing of genomic DNA. 6 μ g of genomic DNA was diluted in 30 μ l milliQ and subsequently fragmented by sonication in icecold water in the Biorupter (Diagenode, Liège, Belgium) for 1 to 4 cycles of 15 times 10 sec ON-10 sec OFF (and briefly spinned down after each cycle) until the fragment size was between 100-600 bp, with the majority of the fragments around 250 bp. Fragment size was verified by electrophoresis on a 2% agarose gel (A9539, Sigma-Aldrich, St. Louis, MO, USA). Methylated DNA immunoprecipitation (MeDIP). This method was previously described by 34 . 5.5 μ g of fragmented genomic DNA was diluted to 400 μ l in TE buffer (10 mM Tris.HCl pH 7.5, 1 mM EDTA) and denatured for 10 minutes at 95 °C, followed by immediate cooling on ice for 5 minutes. 100 μ l of 5X IP buffer (50 mM Na.H 2 PO 4 pH 7, 0.7M NaCl, 0.4 mM Triton X-100) and 2.5 μ l 5-mC monoclonal antibody (C15200006, Diagenode) were added and incubated overnight on a rotating platform at 4 °C. Protein A/G PLUS-Agarose Immunoprecipitation Reagent (sc-2003, Santa Cruz Biotechnology, Inc, Dallas, TX, USA) was pre-washed with 0.1 mg/ml BSA in 1X PBS and resuspended in 40 μ l IP buffer. Beads were added to the DNA-antibody complex and incubated for 2 h on a rotating platform at 4 °C. Beads bound to the DNA-antibody complex were washed 3 times with 1 ml 1X IP buffer. To release the methylated DNA from the beads, the beads were resuspended in 250 μ l digestion buffer (50 mM Tris.HCl pH 8, 10 mM EDTA, 0.5% SDS), 3.5 μ l proteinase K (00120437, Thermo Scientific) and incubated overnight on a rotating platform at 55 °C. The methylated DNA was purified by first mixing with 250 μ l phenol solution (P4557, Sigma-Aldrich) and centrifugation (14,000 rpm, 5 min, room temperature) and repeating this with addition of 250 μ l chloroform:isoamyl alcohol (24:1, C0549, Sigma-Aldrich). Washing was performed by addition of 2 μ l GlycoBlue Coprecipitant (AM9515, Applied Biosystems, Foster City, CA, USA), 20 μ l 5M NaCl and 500 μ l 99.5% icecold ethanol and placement at − 20 °C for at least 1 hour. After centrifugation (14,000 rpm, 15 min, 4 °C), the supernatant was removed and the pellet was washed with 500 μ l of icecold 70% ethanol and centrifuged again. The pellet was resuspended in 30 μ l milliQ. The DNA concentration was measured on the NanoDrop 1000 Spectrophotometer (Thermo Scientific). The yield was about 10-30 ng/μ l.
Whole Genome amplification. 30 ng of each MeDIP sample was amplified using the GenomePlex Complete Whole Genome Amplification kit (WGA2, Sigma-Aldrich) according the manufacturer's instruction, except no fragmentation was performed, as the samples were already sonicated. After amplification, the DNA was cleaned using the GeneJET PCR Purification kit (K0701, Thermo Scientific). The DNA concentration was measured on the NanoDrop 1000 Spectrophotometer (Thermo Scientific) and the samples were diluted to 20 ng/μ l.
Selection of CpG rich regions around transcription start site of interesting genes and primer design. For each gene of interest, several CpG rich (> 4 CpG's) fragments (150-200 bp) were selected in the promoter region (5000 bp upstream of transcription start site) and around the transcription start site and the corresponding primers for these fragments were designed (Supplementary Table S2 online) using Primer Designing Tool (NCBI). The primers were purchased from Life Technologies (Carlsbad, CA, USA) and verified as described previously. Statistical Analysis. qPCR and MeDIP-qPCR data were processed with the statistical software package SAS version 9.3 (SAS Institute Inc., Cary). The gene expression values (2 −ddCt ) and Ct values of the MeDIP-qPCR were analyzed using one-way ANOVA with treatment (non-manipulated, sham-manipulated and albumen-deprived) as independent variable. When a significant effect of treatment was demonstrated, the means were further compared by a posthoc Tukey's test. All data are shown as mean ± SEM.