extfieldsofstudy: Medicine Biology

text:
A combined reference panel from the 1000 Genomes and UK10K projects improved rare variant imputation in European and Chinese samples

Imputation using the 1000 Genomes haplotype reference panel has been widely adapted to estimate genotypes in genome wide association studies. To evaluate imputation quality with a relatively larger reference panel and a reference panel composed of different ethnic populations, we conducted imputations in the Framingham Heart Study and the North Chinese Study using a combined reference panel from the 1000 Genomes (N = 1,092) and UK10K (N = 3,781) projects. For rare variants with 0.01% < MAF ≤ 0.5%, imputation in the Framingham Heart Study with the combined reference panel increased well-imputed genotypes (with imputation quality score ≥0.4) from 62.9% to 76.1% when compared to imputation with the 1000 Genomes. For the North Chinese samples, imputation of rare variants with 0.01% < MAF ≤ 0.5% with the combined reference panel increased well-imputed genotypes by from 49.8% to 61.8%. The predominant European ancestry of the UK10K and the combined reference panels may explain why there was less of an increase in imputation success in the North Chinese samples. Our results underscore the importance and potential of larger reference panels to impute rare variants, while recognizing that increasing ethnic specific variants in reference panels may result in better imputation for genotypes in some ethnic groups.

In the current study, we used a relatively new reference panel with 4,873 whole genome sequences (1,092 from 1000 G and 3,781 from UK10K Project). We evaluated the imputation accuracy for both common and rare SNPs in a large-scale Caucasian Cohort study, the Framingham Heart Study (FHS), by comparing the imputed genotypes and actual genotypes from SNP arrays. In addition, previous studies also demonstrated that using ethnicity-specific reference panels to impute genotypes of the targeted samples with the same ancestral genetic background might improve imputation accuracy especially for those ethnic-specific variants 14 . The improvement of imputation accuracy is due to the fact that LD patterns of the ethnic-specific variants may not be captured by other ethnic groups with different ancestral genetic background. However, for the SNPs that are not ethnic-specific, using a reference panel from other ethnic groups or mixed ethnic groups may still provide accurate imputation. To evaluate this hypothesis, we performed whole genome imputation to generate un-genotyped variants in an East Asian cohort, North Chinese Study (NCS), using our new reference panel (1000 G + UK10K, with a majority of sequences from samples of European-descent) and estimated imputation quality of SNPs grouped into different ancestral-distance between East Asians and Europeans estimated by population differentiation score (fixation index, "F ST ").

Results
Genotype imputation of FHS data. We performed whole genome imputation in 8,700 FHS participants with 536,217 genome-wide genotypes from Affymetrix 550 K SNP-chips. The imputations with 1000 G and the combined 1000 G + UK10K generated 29,605,432 and 39,849,120 variants including the input genotypes (Table 1). There were 33% more variants generated with the 1000 G + UK10K panel than with the 1000 G panel. The imputations with 1000 G and 1000 G + UK10k added 15,245,172 and 21,449,101 variants to 105,796 low-frequency (MAF 1-5%) and rare (MAF 0.01-1%) actual genotypes, and added 8,068,486 and 11,639,300 variants to 1,578 actual genotypes with MAF less than 0.01%. Among imputed variants, 65.5% and 67.6% were variants with 0 < MAF ≤ 1% when imputed from the 1000 G and the 1000 G + UK10K panels respectively. The most abundant imputed variants were the variants with 0.01% < MAF ≤ 0.5%. The complete numbers of actual genotypes and imputed variants stratified by MAFs are listed in Table 1.

Quality of FHS imputation.
We evaluated the quality of the FHS imputation by median squared correlations (r 2 ) between actual allelic dosages from the OMNI5 chip and allelic dosages from imputations with 1000 G and 1000 G + UK10K. Figure 1 displays the r 2 of 2,752,953 variants within categories of MAF estimated using the 1000 G + UK10K imputation results. The median r 2 of imputation with 1000 G + UK10K was higher than the median r 2 of 1000 G in most of MAF ranges except the very rare alleles (0 < MAF ≤ 0.01%), but the median r 2 in the very rare MAF range was very low, e.g. 0.01 and 0 in 1000 G and 1000 G + UK10K imputations. The median r 2 for SNPs with MAF > 5% were 0.95 and 0.96 in imputations of 1000 G and 1000 G + UK10K, respectively. Imputation using the 1000 G + UK10K had the greatest r 2 improvement of 0.22 for variants with MAF 0.5-1%, while variants with MAF 0.01-0.5% had 0.13 increase in the median r 2 . Figure S1 shows a table with detailed breakdown of FHS imputation quality. We also evaluated imputation quality when excluding SNPs with MAF ≤ 2% from the measured genotypes. The median r 2 of imputation with 1000 G reference panel were 0.16, 0.27, and 0.4 within MAF ranges 0.01-0.5%, 0.5-1%, and 1-2%, respectively. The median r 2 of imputation with 1000 G + UK10K reference panel were 0.29, 0.48, and 0.57 within the same MAF ranges. The detailed results are shown in Supplementary Table S3. In addition, we also performed imputations by using DuoHMM to take into account family relations in FHS samples. With or without family information, the median r 2 of 1000 G + UK10K imputation both improved by 0.19 for variants of MAF 1-2%. Family information slightly improved 1000 G + UK10K imputation from 0.57 to 0.64. The results are shown in Supplementary Table S4.
Increased proportion of well-imputed variants in FHS imputation with 1000 G + UK10K reference panel. We compared proportions of well-imputed variants across categories of imputation quality scores (INFO score, estimated by IMPUTE2) in all MAF ranges. Since most GWAS analyses commonly exclude low quality imputed variants with INFO scores lower than 0.4, we used INFO scores of 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 to evaluate proportions of well-imputed variants. As shown in Table 2, overall, imputation with 1000 G + UK10K had a larger proportion (increased from 0.4% to 27.5%) of well-imputed variants across all INFO and MAF ranges. Imputed variants with INFO score ≥ 0.4 and MAF > 0.5% all had well-imputed variants over 90%. Among imputed variants with INFO score ≥ 0.4 and MAF 0.01-0.5%, imputation with 1000 G + UK10K produced 76.1% well-imputed variants, which represented a greater increase (13.2%) compared to imputation with 1000 G in the same INFO range (62.9%). Imputed variants with INFO score ≥ 0.4 and 0 < MAF ≤ 0.01% had a 9.1% increase in well-imputed variants. Although imputed variants with INFO score ≥ 0.4 and MAF > 5% had only a 0.4% increase in well-imputed variants, their well-imputed variants were both over 99%.
Across all INFO and MAF ranges, the largest increase in well-imputed variants was 27.5% in the category of INFO score ≥ 0.6 and MAF 0.5-1%. Among imputed variants with INFO score ≥ 0.7, there was a 25.4% increase in well-imputed variants with MAFs 0.5-1%. For the imputed variants with INFO score ≥ 0.7 and MAF 0.5-1%, 64.2% were well-imputed with 1000 G + UK10K panel, and only 38.8% were well-imputed with 1000 G. If we used a more stringent INFO score cutoff of 0.8, variants imputed using the 1000 G + UK10K had 16.9% and 16.8% increases in MAF ranges of 0.5-1% and 1-5%. In the most stringent INFO score of 0.9, variants with MAF 1-5% had 11.9% increase compared to imputation with 1000 G. In the most rare MAF range of 0 < MAF ≤ 0.01%, if we used INFO ≥ 0.7, 0.8, and 0.9, imputation with 1000 G + UK10K had 3.6% to 11.6% well-imputed variants, but imputation with 1000 G did not produce any well-imputed variants.
Distribution of the well-imputed variants with functional annotation. We investigated the proportion of 1000 G + UK10K well-imputed (INFO score ≥ 0.4) variants according to their functional roles. The well-imputed variants were classified using ANNOVAR 16 with human genome version hg19. The classes of variants included those in exons, introns, non-coding RNA and intergenic regions. Table 3 shows the proportions of well-imputed variants of different types stratified by MAF. The proportion was the number of imputed variants with INFO ≥ 0.4 divided by total number of imputed variants. In general, the distributions of well-imputed variants were similar regardless of the assigned functional roles. The variants with MAF greater than 0.5% had a relatively large percentage of well-imputed variants, and the well-imputed proportions were at least 92.7% even for variants in non-coding and intergenic regions. The variants with MAF less than 0.5% had lower proportions of well-imputed variants. Almost all variants annotated as exonic splicing were well-imputed even in the very rare MAF range. Except exonic splicing, the well-imputed variants were evenly distributed in all MAF ranges.
Genotype imputation of NCS data. We performed whole genome imputation in 3,042 NCS participants with 527,984 genome-wide genotypes from Affymetrix Axiom CHB 1 (657 K) array. The imputations with 1000 G and 1000 G + UK10K generated 29,198,055 and 40,611,832 variants including the input genotypes ( Table 4). The imputation with 1000 G + UK10K had 39% more variants comparing to the imputation with 1000 G only. The imputations with 1000 G and 1000 G + UK10k added 14,040,190 and 16,151,464 variants to 16,173 low-frequency (MAF 1-5%) and rare (MAF 0.01-1%) actual genotypes, and added 9,753,722 and 18,276,702 variants to none actual genotype with MAF less than 0.01%. Among imputed variants, 68% and 66.9% were variants with  Table 2. Proportion of well-imputed variants in FHS imputation results. The imputations were performed with 1000 G and 1000 G + UK10K reference panels.
0 < MAF ≤ 1% when imputed from the 1000 G and the 1000 G + UK10K panels respectively. The most abundant imputed variants with 1000 G and 1000 G + UK10K were variants with MAF 0.01-0.5% and MAF 0-0.01%, respectively. The complete numbers of actual genotypes and imputed variants stratified by MAFs are listed in Table 4.

Quality of NCS imputation.
Due to lack of a second source of genotype data to evaluate imputation quality in Chinese sample, we used squared correlation (R 2 ) between actual allelic dosages and imputed allelic dosages by treating the actual genotypes as unknown genotypes. We stratified the R 2 of 514,707 variants with MAFs estimated in 1000 G + UK10K imputation results in Fig. 2. We did not observe increased median R 2 when performing imputation with 1000 G + UK10K. The median R 2 were all above 0.92 among the low-frequent and common variants.
We then evaluated imputation quality of rare variants (MAF < 1%) by median squared correlation (r 2 , different from R 2 that provided by IMPUTE2) between actual allelic dosages and imputed dosages, 335 SNPs with MAF < 1% (which were excluded before imputation). We observed imputation with 1000 G + UK10K decreased imputation quality by 0.19 in terms of median r 2 . The median r 2 of imputation with 1000 G reference panel is 0.51 while the median r 2 of imputation with 1000 G + UK10K reference panel is 0.32.

Increased proportion of well-imputed variants in NCS imputation with 1000 G + UK10K reference panel.
We compared the proportions of well-imputed variants across categories of imputation quality in a similar fashion to the FHS imputation above. We used INFO scores of 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 to evaluate proportions of well-imputed variants. As shown in Table 5, overall, imputation with 1000 G + UK10K had a larger proportion (increased from 0.2% to 14%) of well-imputed variants across all INFO and MAF ranges. Imputed variants with INFO score ≥ 0.4 and MAF > 0.5% all had well-imputed variants over 87%. Among imputed variants with INFO score ≥ 0.4 and MAF 0.01-0.5%, imputation with 1000 G + UK10K produced 61.8% well-imputed variants, which was the largest increase (12%) compared to imputation with 1000 G in the INFO range (49.8%). Imputed variants with INFO score ≥ 0.4 and 0 < MAF ≤ 0.01% had a 0.8% increase in well-imputed variants. Although imputed variants with INFO score ≥ 0.4 and MAF > 5% had only a 0.4% increase in well-imputed variants, their well-imputed variants are both over 99.5%.
Across all INFO and MAF ranges, the largest increase in well-imputed variants was 14% in the category of INFO score ≥ 0.5 and MAF 0.01-0.5%. Among imputed variants with INFO score ≥ 0.6, there was a 12.7% increase in well-imputed variants with MAFs 0.5-1%. For the imputed variants with INFO score ≥ 0.6 and MAF 0.5-1%, 60.3% were well-imputed with 1000 G + UK10K panel, and only 47.6% were well-imputed with 1000 G. If we used a more stringent INFO score cutoff of 0.8, variants imputed using the 1000 G + UK10K had 5.8% increases in MAF ranges of 1-5%. In the most stringent INFO score of 0.9, variants with MAF 1-5% had 4% increase compared to imputation with 1000 G. In the most rare MAF range of 0 < MAF ≤ 0.01%, both imputations with 1000 G and 1000 G + UK10K had less than 1% well-imputed variants even when we used the most loosen INFO of 0.4.

Imputation quality and Population Difference Estimated by fixation index (F ST ). Compared to
imputation results in the FHS, we observed smaller increases in well-imputed rare variants in NCS imputation results. The proportion of well-imputed variants with MAF 0.1-0.5% increased by 21% in the FHS imputation in contrast to only a 9% increase in NCS when imputations were done using the 1000 G + UK10K panel. We first checked the population structure of all individuals by principal component analysis (PCA). In the PCA result shown in Fig. S2, FHS data overlapped well with the 1000 G European and UK10K data, and the NCS data overlapped with the 1000 G Asian data. F ST , which describes the population differentiation due to genetic structure, was then calculated for 14,339,925 variants existing both in FHS and NCS imputation results to investigate the relationship between imputation quality and F ST . Figure 3 shows the INFO score of imputed variants in both the FHS and NCS cohorts. The pattern of the plots suggested a decreasing quality of the NCS imputation results (dots shifted toward left) whereas the FHS imputation results maintained good quality (dots remained above the line of identity). For example, in the plot with MAF 0.1-0.5% and F ST 0.09-0.1, most imputed variants in FHS had INFO scores higher than 0.75, and imputed variants in NCS had INFO scores widely distributed (from 0.38 to higher). Computing cost. Genome-wide imputation of FHS genotypes with 1000 G reference panel took 2,422 computing hours in total (with 773 parallel jobs), and imputation with 1000 G + UK10K took 4,688 computing hours (with 908 parallel jobs).

Discussion
Our ability to comprehensively ascertain both common and rare genetic variants may bring us a step closer to understanding the genetic architecture of a disease. In addition, rare genetic variants may explain missing heritability that common variants fail to explain. Thus, identifying rare variants associated with phenotypes and diseases become a more important task when investigating genetic risk factors. Whole genome sequencing is the best way to find rare variants in the genome. However, the cost of whole genome sequencing for a significant amount of individuals at this time is still expensive, and there remains considerable interest by large-scale enterprises, such as GEFOS 17 and UK Biobank 18 , in analyzing imputed genotypes. Our work suggests that increasing sample size of the reference panel increases imputation quality, especially for rare variants. Huang et al. recently demonstrated the same finding when comparing imputation results between 1000 G and 1000 G + UK10K 19 . In addition, our findings also suggest that even if the reference panel is not ethnic-specific for the targeted samples, the imputation quality is still improved. Comparing the improvement of the imputation quality between Chinese samples and Caucasian samples by using 1000 G + UK10K reference panel, we found that the improvement in imputation quality of Chinese samples was not as good as in the Caucasian samples. This was not surprising since the 1000 G + UK10K reference panel has an increased proportion of similar (Caucasian) ethnic samples (adding 3,000 Caucasian samples from the UK10K project). We suggest that to increase imputation quality, the best strategy is to use reference panels with the specific ethnic composition of the study population being imputed as this may improve the imputation of relatively rare SNPs.
In this study, we investigated genotype imputation improvement in terms of quality and quantity especially on rare variants with the 1000 G + UK10K reference panel in Caucasian and Chinese cohorts. The imputations with 1000 G + UK10K reference added 9,783,688 and 11,413,777 additional variants in total to the imputation using the 1000 G reference panel alone in FHS and NCS, respectively. Our findings suggested that a larger reference panel improved imputation performance, especially for less-common and rare variants, even for Chinese samples. On the other hand, in the FHS, the imputation quality for common variants (MAF > 5%) only improved by 0.01 according to the r 2 estimated between actual allelic dosages and imputed allelic dosages. One potential explanation for this finding is that using the 1000 G as a reference panel already provides high quality imputation (r 2 = 0.95) for common variants. In NCS imputation, for common variants, we did not see much overall improvement in quality in the imputation with the 1000 G + UK10K reference panel. Our observation confirms a recent study that demonstrated the value of study-specific reference panels 20 . However, even when using INFO score ≥ 0.4 as cutoff, the 1000 G + UK10K panel still increased the fraction of well-imputed variants in NCS samples. Number of well-imputed SNPs (with INFO score ≥ 0.4, using 1000 G + UK10K reference) with MAF 0.01-0.5% increased from 62.9% to 76.1% and 49.8% to 61.8% in FHS and NCS, respectively. We still observed increased proportions of well-imputed variants when we used more stringent INFO score cutoff.   We used variants' F ST scores representing the population differentiation between Caucasian and Chinese to explain why the imputation of NCS samples did not have comparable qualities to the imputation of FHS samples based on 1000 G + UK10K reference panel. Our findings suggest that some rare variants are commonly existing in both ethnicities and these variants may retain haplotype structure across ethnicities which allow us to impute some rare variants with high quality given the underlying genetic background is different between reference panel and targeted samples.
There are potential limitations in this study. First, we compared imputation quality between two different-ethnicity samples, Chinese and Caucasian. However, the difference between these two samples include sample size, genotype platform and density of SNPs, which may all affect imputation quality, in addition to the differences in genetic background. Our results may not be applicable to other ethnic groups, such as African populations due to their short haplotype blocks compared to Caucasian and East Asian populations. Second, major imputation results of family-based FHS were not performed using family information. However, results of a small set FHS imputation with family information did not change our conclusion on imputation improvement with 1000 G + UK10K reference panel. Third, this study only used genotyping data generated from Affymetrix chips to impute genotypes. We cannot estimate the influences that may possibly come from different density of variants or including rare variants in their genotyping platforms. We thus tested a small set of FHS imputation excluding variants with MAF ≤ 2% in the input measured genotypes. The imputation quality without rare variants was consistent with the findings when actual genotypes with MAF ≤ 2% were included. For 335 rare variants of Chinese samples (with MAF < 1%), we observed decreased imputation quality when 1000 G + UK10K were used. More data on Chinese-specific rare variants are needed to confirm this source for decrease in quality. We didn't evaluate whether even larger reference panels such as the human haplotypes provided by the Haplotype Reference  Table 5. Proportion of well-imputed variants in NCS imputation results. The imputations were performed with 1000 G and 1000 G + UK10K reference panels. Consortium (HRC, http://www.haplotype-reference-consortium.org/), will improve a lot of the imputation quality of the relatively rare variants. The HRC was not yet available when we performed our analyses. HRC panel consists of 64,976 haplotypes with 39,235,157 SNPs. Majority of the HRC references are Caucasian samples, thus, we expect such larger reference panel may further increase imputation quality for the rare variants in Caucasian samples, but may not necessarily improve imputation quality of non-Caucasian samples.
Although the imputation quality for rare variants improved in our study, still, more than half of the rare variants were not imputed with considerable high quality (having INFO score < 0.4); thus, these variants cannot be used in association analyses. It will be necessary to estimate how large the sample size of the reference panel should be, in order to have more than 80% of imputed rare variants with INFO score ≥ 0.4. We were also not able to evaluate the effect of sample size in the targeted samples to the imputation quality.
In summary, we evaluated imputation quality by imputing un-genotyped SNPs in a Caucasian sample and a North Chinese sample using a combined reference panel (N = 3781) by combing whole genome sequences from 1000 G and UK10K projects. We found that our proposed reference panel outperformed the 1000 G reference panel, especially for relatively rare variants in the Caucasian sample, suggesting that increasing sample size in a reference panel, especially ethnic specific sample size (the same ethnicity of the targeted sample) would improve imputation quality. Nevertheless, the imputation quality was still improved in common variants for non ethnic-specific samples (for example, imputation on North Chinese samples in current study) although the improvement is limited.

Methods
Reference panel: 1000 Genome project (1000 G). The publicly available database of the Phase 1 version 3 of the 1000 G contains 1,092 individuals including the following populations: 181 Admixed American, 246 African, 286 East Asian, and 379 European 7 . The data were downloaded from the IMPUTE website (ALL_1000 G_phase1integrated_v3_impute_macGT1.tgz, build date: 2012-08-25). In the 1000 G data set, monomorphics, singleton and doubleton single nucleotide variants (SNVs with minor allele count less than two) were removed. A total of 30,061,896 SNVs were used for imputations.
Reference panel: UK10K project (UK10K). The database from the UK10K project had a total of 3,781 individuals who underwent whole genome sequencing from the Avon Longitudinal Study of Parents and Children (ALSPAC) and the King's College London Department of Twin Research and Genetic Epidemiology Twins Registry (TWINSUK). The UK10K data (IDs: EGAD00001000195 and EGAD00001000194) were downloaded from the European Genome-Phenome Archive. The two data sets were first merge by PLINK 21 and then haplotype phased by SHAPEIT2 22 . We excluded monomorphics, singleton and doubleton single nucleotide variants (SNVs with minor allele count less than two) from the UK10K data set. A total of 27,793,639 SNVs were used for imputations.  Table S1 shows the three reference panels.

Study Sample and Genotyping of Framingham Heart Study (FHS).
The FHS contains family members and unrelated participants. The participants are Caucasian of European descent. Two sets of genotyping data, Affy550K and OMNI5 were generated by the FHS SHARe (SNP Health Association Resource) project. The first set, Affy550K, with 549,827 SNPs (Affymetrix 500 K mapping array plus Affymetrix 50 K supplemental array) came from genotyping in over 9,274 FHS subjects from over 900 families 24 . By estimation, we expected at least 80% genomic coverage (pair-wise genotype correlation r 2 > 0.8) of the HapMap Phase I + II common SNPs (minor allele frequency, MAF ≥ 0.05) for the Caucasian population 25 . We re-confirmed the kinship relationship among the Framingham pedigrees using all available SNPs. PEDCHECK 26 and RELPAIR 27 programs were used to check for errors in Mendelian inheritance and inconsistencies within pedigrees. We excluded 793 individuals with an average SNP call rate of < 0.97. We also excluded individual SNPs with a call rate of < 0.95 (34,868 SNPs; 6.3%); Hardy-Weinberg equilibrium (HWE) test p-value < 1 × 10 −6 (8,531 SNPs; 1.6%); MAF < 0.01 (66,829 SNPs; 12.2%); or unknown genomic annotation (6, Haplotype phasing and genotype imputation tools. Input genotype data sets for imputation were first pre-phased with SHAPEIT2 22 , and then missing genotypes were imputed with IMPUTE2 2 by 5-Mb region at a time against the 1000 G or the combined reference panel. To evaluate the effect of pedigree information to imputation, we post-process phased haplotypes with duoHMM to incorporate pedigree information.

Estimation of overall genotype imputation quality by squared Pearson correlation (r 2 and R 2 ).
To estimate overall genotype imputation quality, we calculated squared Pearson correlation (r 2 ) between actual (discrete) allelic dosage and the imputed (continuous) allelic dosage over a defined set of samples. For FHS, we excluded overlapping genotypes both existing in Affy550K and OMNI5 data, and the remaining genotypes overlapping between 1000 G and 1000 G + UK10K panels were used to examine imputation quality by r 2 . For NCS, due to lack of a second set of genotyping data, we estimated imputation quality by squared correlation (R 2 ) between input actual allelic dosage and the allelic dosage imputed for the input genotypes. The imputed allelic dosages were generated with the input genotype masked. The R 2 measurements of NCS were provided by IMPUTE2. Principle component analysis (PCA). We performed PCA to find directions of the largest variability of the genotype data. We used the first two principal components to reflect the geographical distance of all individuals from reference panels and our study samples. The PCA was performed using chromosome 20.

Population differentiation measured by fixation index (F ST ).
We used a F ST to describe the population differentiation due to genetic structure at the level of the SNP. F ST of zero indicates that the SNP has no divergence between populations. To analyze the imputation results of FHS and NCS with the combined reference panel, we used pairwise F ST scores between Utah Residents with Northern and Western European ancestry (CEU) and Han Chinese in Beijing (CHB) calculated by the method of Weir and Cockerham 29 . The F ST scores of CEU and CHB were calculated using the 1000 G project data 30 .


extfieldsofstudy: Biology Medicine

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