Introduction

Non-alcoholic fatty liver disease (NAFLD) is a hepatic manifestation of metabolic syndrome with a global prevalence of 24%, including South America and the Middle East, followed by Asia, the USA, and Europe1,2,3. Evidence to support the relationship between nutritional imbalances in the early developmental period and a predisposition for NAFLD in later life is increasing4,5,6,7,8,9. Sandboge et al. examined 1587 aged participants from the Helsinki Birth Cohort Study and showed that birth and childhood body sizes were negatively associated with NAFLD outcomes and also that individuals who had been small in early life and obese as adults were at the highest risk of developing NAFLD6. Faientza et al. reported that NAFLD was detected in 34.8% of children who were born as small for gestational age (SGA), but not in those born as appropriate for gestational age7. Alisi et al. showed that SGA was present in 38.9% of children complicated with NAFLD and in 6.7% of uncomplicated children8. Moreover, Bugianesi et al. found that a low birthweight increased the likelihood of severe steatosis in pediatric NAFLD4. Ours and previous studies demonstrated that maternal global nutrient restrictions primed the deterioration of hepatic steatosis in adulthood in ovine10, rat11, and mouse12 offspring. Collectively, the findings of human cohorts and animal studies support the concept that undernourishment (UN) in utero primes the risk of augmented hepatic fat deposition in later life, particularly with an obesogenic diet. However, the underlying mechanisms remain unclear. Alisi et al. recently demonstrated the importance of pre- and postnatal environmental monitoring for protection against NAFLD9; however, specific therapeutic strategies have not yet been established, particularly after the progression of NAFLD originating from nutritional imbalances during early life.

We established mouse models of UN in utero by maternal caloric retraction and developing the phenotypes of various non-communicable diseases12,13,14,15,16, and subsequently demonstrated that treatments with the hydrophilic secondary bile acid tauroursodeoxycholic acid (TUDCA), an endoplasmic reticulum (ER) stress alleviator, markedly ameliorated developmentally-deteriorated hepatic steatosis12. The TUDCA treatment was only effective for pups with UN in utero. A previous study reported that a treatment with TUDCA improved many acute and chronic diseases17; however, the underlying mechanisms have not yet been elucidated.

In the developmental origins hypothesis, epigenetic modifications, such as DNA methylation and histone modifications, are considered to play critical roles in mediating how the early life environment impacts on later health and susceptibility to non-communicable diseases18,19. Moreover, increasing evidence has revealed the extensive involvement of epigenetic modifications in the pathophysiology of NAFLD20,21.

Therefore, we hypothesize that (1) UN in utero may epigenetically program the expression of some genes, by DNA methylation and/or histone modifications, in the process of the deterioration of hepatic steatosis in offspring in a mouse model of UN in utero under an obesogenic diet, and (2) a treatment with TUDCA remodels these epigenetic modifications, concomitant with the amelioration of advanced hepatic steatosis originating from UN in utero, as an example of epigenetic plasticity in the developmental programming of specific genes.

To prove these hypotheses, we investigated genetic profiles in the liver using a microarray analysis and epigenetic regulation by methyl-binding domain (MBD) protein sequencing and chromatin immunoprecipitation (ChIP) assays.

Results

Deterioration of hepatic steatosis by UN in utero under the obesogenic diet

UN induced a significant deterioration in hepatic steatosis after HFD (cohorts 2 and 3), but not before HFD (cohort 1) (Fig. 1 and Supplementary Fig. S1A–K). The administration of TUDCA significantly ameliorated hepatic steatosis for UN in utero only (Figs 1C, 2A). Figure 2 shows increases in LD sizes along with the deterioration after HFD, and the restoration of their size after the administration of TUDCA. These results were consistent with our previous findings12.

Figure 1
figure 1

Liver and body weight changes by HFD simulating hepatic steatosis. Data are expressed as means and error bars indicate standard deviations (SD) in cohort 1 (A), cohort 2 (B), and cohort 3 (C). Significant differences were observed using the Student’s t-test (A,B) or Steel-Dwass test (C) (p* < 0.05, p** < 0.01, p*** < 0.001).

Figure 2
figure 2

Enlargements in LD are consistent with the deterioration of hepatic steatosis and restoration with TUDCA in cohort 3. Under ×40 magnification; (A) HE stain, (B) Oil Red O stain (inset with ×80 magnification), (C) mean area of LD, (D) mean diameter of LD with error bars indicating SD. Significance is stated as p** < 0.01 calculated by the Steel-Dwass test.

Changes in gene expression profiles with or without the TUDCA treatment

To examine the entire genetic profile, we performed a microarray analysis of all 3 cohorts (at the end of 9, 17, and 22 weeks). A total of 34,472 genes were expressed among all groups. The microarray analysis of cohorts 1 and 2, before and after starting HFD, revealed that 8 and 82 genes were differentially expressed, respectively (Linear fold change ≤−2 or ≥2 fold, ANOVA p < 0.05, UN-Veh vs. NN-Veh) (Supplementary Table S1). Similarly, in cohort 3 at 22 weeks, we examined 53 genes (Linear fold change ≤−2 or ≥2 fold, ANOVA p < 0.05, UN-Veh vs. NN-Veh) shown with a Volcano plot (Fig. 3A). The TUDCA treatment also altered the expression of genes in NN-TU (69 genes; Linear fold change ≤−2 or ≥2 fold, ANOVA p < 0.05) and UN-TU (40 genes; Linear fold change ≤−2 or ≥2 fold, ANOVA p < 0.05) (Fig. 3B). In this cross-sectional assessment, we stratified these genes to select 15 candidate genes in cohort 3 that were significantly up- or down-regulated in UN in utero and restored by TUDCA, as listed in Supplementary Tables S1 and S2. We also performed a longitudinal assessment of microarray data between cohorts 1 and 2. We studied 133 genes, listed in Supplementary Table S3, which showed significant changes in comparisons between UN before HFD (cohort 1) and UN after HFD (cohort 2). We then contrasted them with those unaltered between NN before HFD (cohort 1) and NN after HFD (cohort 2). Therefore, in this present study, we enlisted 9 genes of interest (GOI) in cross-sectional (cohort 3) and longitudinal analyses (between cohorts 1 and 2) (Table 1). We submitted our microarray data to the GEO repository, which is approved under the accession number GSE123733.

Figure 3
figure 3

Genetic profiling of developmentally-programmed hepatic steatosis by a microarray analysis; differentially expressed mRNA in UN in utero and alterations by TUDCA. A Volcano plot of 34,472 genes expressed in all three cohorts. Red and green dots were up- and down-regulated genes, respectively, in cohort 3 (Linear fold change of ≤−2 and ≥2 and ANOVA P < 0.05). Cidea and Cidec genes are encircled on the plot, and were up-regulated in UN-Veh and then down-regulated with the oral gavage of TUDCA.

Table 1 List of genes of interest (GOI).

Cell Death-Inducing DNA Fragmentation Factor-Like Effectors A (Cidea) and C (Cidec) were included in both of the analyses. This particular pattern of expression of these two genes was observed between 9 and 22 weeks (Supplementary Table S4). Relative quantitative RT-PCR showed that the mRNA expression of the Cidea and Cidec genes in cohorts 1, 2, and 3 was consistent with microarray data (Fig. 4A–C). Cidea and Cidec are known to increase LD sizes, thereby contributing to lipid storage. Changes in LD sizes in cohorts 2 and 3 were consistent with those in the gene expression of Cidea and Cidec (Fig. 2B–D and Supplementary Fig. S1L-O). Moreover, the presence of Cidea and Cidec proteins was widely observed in livers with steatosis by immunohistochemistry (Fig. 4D,E).

Figure 4
figure 4

Gene expression and immunohistochemistry of Cidea and Cidec. Quantitative RT-PCR of Cidea and Cidec in cohort 1 (A), cohort 2 (B), and cohort 3 (C). Immunohistochemistry of Cidea (D) and Cidec (E). Significance was observed by the Student’s t-test (A,B) or Steel-Dwass test (C) (p* < 0.05 and p*** < 0.001).

Therefore, we performed an epigenetic assessment of Cidea and Cidec because of their high fold changes and potential ability to regulate hepatic fat deposition.

DNA methylation alterations in Cidea and Cidec genes

Large methylation variations were observed among the four groups and a list of methylated genes was shown in Supplementary Table S5. We examined the entire genome lengths of the Cidea (chr18:67320000–67410000) and Cidec (chr6:113415000–113450000) genes, including up- and downstream regions. No significant differences were observed in methylation peaks among the four groups before and after the TUDCA treatment (Supplementary Fig. S2).

Effects of UN in utero on ChIP assay results around Cidea

In cohort 3, we performed a ChIP assay to examine the acetylation and methylation levels of H3 and H4 around the Cidea and Cidec genes (Fig. 5 and Supplementary Fig. S37). Around the Cidea gene, a modification known to inhibit gene expression22, i.e. the level of dimethylation of H3K27, was significantly suppressed at the 3′ end (Fig. 5A). This result was consistent with the elevated mRNA expression levels of Cidea in the microarray assay as well as qPCR for UN in utero (Fig. 4C).

Figure 5
figure 5

Histone modifications by the ChIP assay. Effects of UN in utero on the dimethylation of H3K27 around Cidea (A) and dimethylation of H3K27 (B) and H3K36 (C) around Cidec. Effects of TUDCA and UN in utero on H3K27 around Cidea (D) and the dimethylation of H3K27 (E), H3K36 (F), and H3K4 (G) around Cidec. Mean values were expressed with SD as error bars and significance (*or #p < 0.05, **or ##p < 0.01, ***or ###p < 0.001) was calculated with the Student’s t-test or Mann-Whitney U test as appropriate. The X-axis shows the region of gene length and the y-axis is % input.

Effects of UN in utero on ChIP assay results around Cidec

UN in utero caused changes in various histone modifications around Cidec, contributing to strong genetic expression. Around the Cidec gene, the levels of gene suppressor modifications, the dimethylation of H3K27 and H3K3623, were lower near the 3′ end (Fig. 5B,C). These results were consistent with the augmented mRNA expression levels of Cidec in the microarray assay as well as qPCR for UN in utero (Fig. 4C).

Effects of the TUDCA treatment on ChIP assay results around Cidea in UN offspring

TUDCA enhanced the suppressor modification H3K27 dimethylation in the entire length, enriched at the 3′ end, to inhibit the mRNA expression of Cidea, as shown by the microarray analysis and qPCR (Fig. 5D). This was consistent with the decrease in the mRNA expression levels of Cidea in the microarray assay as well as qPCR for UN in utero with the TUDCA treatment (Fig. 4C).

Effects of the TUDCA treatment on ChIP assay results around Cidec in UN offspring

The H3K27 and H3K36 dimethylation modifications around Cidec were markedly enhanced by TUDCA in UN in utero, and were also observed near the 3′ end (Fig. 5E,F), restoring changes in mRNA expression. Moreover, lower levels of activator modifications, such as H3K4 dimethylation, were also prominent in TUDCA-treated UN in utero pups (Fig. 5G) at the 3′ end, which is suggestive of remodeling. These results were consistent with the low mRNA expression levels of Cidec in the microarray assay as well as qPCR for UN in utero with the TUDCA treatment (Fig. 4C).

Effects of the TUDCA treatment on ChIP assay results around Cidea and Cidec in NN offspring

The TUDCA treatment exerted suppressive (H3K4 and H3K36 dimethylation) and enhancing (H3K27 dimethylation) effects on the gene expression of Cidea and Cidec in NN offspring (Supplementary Figs S6, S7). We currently cannot provide an explanation for the contribution of histone modifications to gene expression in NN offspring because the TUDCA treatment induced changes that were not consistent for the construction of histones (Fig. 4C).

Discussion

In the present study, we performed a microarray analysis of developmentally-deteriorated hepatic steatosis and its amelioration by TUDCA12. Cross-sectional and longitudinal analyses of microarray data were performed and 9 GOI (Table 1) were commonly selected.

Among the 9 GOI examined, the Cidea and Cidec genes were markedly up-regulated and reported to be closely related to the pathophysiology of metabolic disorders, such as diabetes, obesity, and, most importantly, liver steatosis in humans as well as in animals24,25. The gene enrichment analysis with the functional annotation of 9 GOI showed that Cidea and Cidec had the greatest involvement in the function of apoptotic pathways and also LD and lipid particle pathways (Supplementary Fig. S8). Cidea and Cidec have been shown to increase LD sizes, thereby augmenting fat deposition in hepatocytes24,26, and their strong expression by UN and normalization by the TUDCA treatment (Fig. 4C) explained the deterioration induced by UN and amelioration by TUDCA in hepatic steatosis in our mouse model. Immunohistochemistry revealed the wide distribution of the Cidea and Cidec proteins in LD of various sizes in the mouse fatty liver (Fig. 4D,E). Oil Red O staining showed the significant enlargement of LD with UN in utero, and this was markedly reduced by the TUDCA treatment (Fig. 2B–D). The latest definition of hepatic steatosis involves the excessive accumulation of LD, which was previously hypothesized to be triggered by Cide proteins, including Cidea and Cidec26,27,28. Since 7 other GOIs appeared to be of less biological relevance to the pathophysiology of hepatic steatosis than Cidea and Cidec (Supplementary Table S6), we focused on Cidea and Cidec.

The emerging concept of epigenetics has emphasized the impact of gestational nutrition over the regulation of gene expression in the development of non-communicable diseases18,29. The critical involvement of epigenetic modifications has been noted in the mechanistic background of the developmental origins hypothesis18,19 as well as the pathophysiology of NAFLD20,21,30. Therefore, we investigated whether epigenetic modifications around the Cidea and/or Cidec genes are examples of epigenetic plasticity in the UN-induced deterioration as well as TUDCA-induced amelioration of hepatic steatosis in this animal model.

We then investigated significant differentially methylated sites using overlapping peaks by DNA MBD sequencing. Neither maternal caloric restriction nor the TUDCA treatment had any effect on DNA methylation around entire Cidea and Cidec genes, including up- and downstream regions (Supplementary Fig. S2). Therefore, further studies are warranted to evaluate the role of chromatin modifications.

Environmental factors, such as diet, are known to induce gene expression by histone modifications that include acetylation and methylation at lysine (K) and arginine within the histone tail18,29,31. Therefore, we investigated the mono- and dimethylation of H3K9, H3K27, and H3K36, dimethylation of H3K4, trimethylation of H3K9, H3K27, and H4K20, and acetylation of H3K9 and H4 (Fig. 5 and Supplementary Figs S310). Among 12 modifications, significant changes were observed in the dimethylation of H3K27, H3K36, and H3K4 (Fig. 5).

UN in utero suppressed H3K27 dimethylation in the gene body of Cidea, thereby contributing to the up-regulation of mRNA (Fig. 5A). In contrast, the dimethylation of both H3K27 and H3K36 (Fig. 5B,C) was suppressed in the gene body of Cidec, and positively correlated with up-regulated mRNA expression in microarrays and qPCR. Histone modifications around the gene body region have been reported to regulate gene expression via transcription elongation32,33,34,35. The present results showed a greater likelihood of the suppression of inhibitory effects on transcriptional elongation around the gene body34,35.

TUDCA is a secondary bile acid that acts on the H3 protein under non-enzymatic conditions to induce post-translational modifications in humans31. In the present study, the ChIP assay showed a clear relationship between the TUDCA treatment and its impact on the regulation of histone modifications in offspring. For example, TUDCA augmented the dimethylation of H3K2722 around Cidea, which induced an inhibitory effect on transcriptional elongation through upstream to the gene body (Fig. 5D). Similarly, TUDCA accelerated the inhibitory effects of the dimethylation of H3K27 and H3K3634,35 around Cidec on transcriptional elongation in the gene body (Fig. 5E,F). These modifications resulted in a marked decrease in the mRNA expression of Cidea and Cidec, respectively. Gene loci remodeled by TUDCA were not identical to those by UN (Fig. 5). Moreover, TUDCA induced new histone modifications, such as H3K4 dimethylation around the Cidec gene (Fig. 5G), which inhibited stimulatory effects on transcriptional elongation in the Cidec gene body33, leading to the suppression of Cidec gene expression. Therefore, we speculate that TUDCA may not restore, but remodel histone modifications around Cidea and Cidec. Furthermore, TUDCA may restore gene expression by repressing mRNA transcription elongation concomitant with the amelioration of hepatic steatosis. In contrast, we did not detect significant alterations in H3K9 or H3K27 trimethylation and only slight changes in H4K20 by UN in utero, without any effect of TUDCA (Supplementary Fig. S9A–F), even though these have been regarded as important repressors of mRNA transcription22,36,37,38.

Recent studies demonstrated that H3K27 methylation may regulate transcription repression independent of DNA methylation by environmental influences39, which supports the direct involvement of H3K27 methylation around Cidea and Cidec in their gene expression without significant changes in DNA methylation.

The expression of the Cidea and Cidec genes was dormant from the early period up to 9 weeks in all groups. However, HFD between weeks 17 and 22 gradually induced the strong expression of the Cidea and Cidec genes (Fig. 4) and encoded their expression prominently in pups with UN in utero collateral to the deterioration of hepatic steatosis, as evidenced by lipid profiles (Supplementary Fig. S1).

Interesting results were obtained for histone modifications in cohorts 1 and 2 (Supplementary Fig. S10). In Cidea, although H3K27 dimethylation was not affected by UN in utero at 9 weeks (cohort 1), HFD led to partial inhibition, which persisted long-term through to 22 weeks (Supplementary Fig. S10A,B). H3K36 dimethylation was transiently inhibited after HFD at 17 weeks; however, no significant changes were observed at 9 or 22 weeks (Supplementary Fig. S10C,D). This result suggests modifications involved in the plasticity of histone remodeling. H3K27 dimethylation showed the similar persistency of suppression in Cidec (Supplementary Fig. S10E,F). However, H3K36 dimethylation remained unchanged even after HFD, and long-term environmental exposure led to the suppression of this modification (Supplementary Fig. S10G,H). H3K4 dimethylation in Cidec showed minor suppression at 9 weeks, but no effect between 17 and 22 weeks (Supplementary Fig. S10K,L), suggesting that the TUDCA treatment specifically affects H3K4 dimethylation in Cidec in UN offspring (Fig. 5G). This distinctive pattern of expression indicates that these histone modifications are preprogrammed via an unidentified mechanism by the maternal nutritional state in utero and HFD triggered subsequent effects in a stepwise manner.

In the present study, we did not investigate the direct mechanism by which the TUDCA treatment changed histone modifications. However, the ChIP assay on the Histone 3.3 subunit, an important recruiter of methyl-transferase40, did not show any alterations around the Cidea and Cidec genes after the administration of TUDCA (Supplementary Fig. S9G,H). These results suggest that this subtype deposition does not play major roles in TUDCA-induced alterations in H3 methylation; however, we cannot deny that TUDCA may play a role after the recruitment of methyl-transferase around these genes.

TUDCA is a well-established ER stress alleviator and we previously revealed the significant suppression of the ER stress response in the liver by a treatment with TUDCA12. Since ER stress has been reported to modulate the activities of H3K4 methyltransferases in various organs41,42,43, we speculate that reductions in ER stress by TUDCA may induce histone modifications in our study model through unidentified mechanisms. Cidea and Cidec are present in the ER membrane24 and ER stress has been shown to regulate the formation of LD with metabolic disruption44. We intend to investigate the involvement of ER stress in the plasticity of histone modifications around the Cidea and Cidec genes following a treatment with TUDCA.

The effects of TUDCA on histone modifications around Cidec were more prominent than those around Cidea, suggesting that an unidentified mechanism is involved in the developmental programming of Cide proteins and their phenotypic expression. TUDCA-inducing histone re-modifications were observed in accordance with the normalization of Cidea and Cidec gene expression levels as well as the marked amelioration of hepatic steatosis, even under HFD in UN, but not NN offspring (Supplementary Figs S67).

There were some limitations in the present study. We used the cut-off value of a two-fold change and ANOVA P values in the microarray analysis, in accordance with our previous studies45,46; however, we cannot fully deny the possible contribution of genes with less than two-fold changes. We were unable to assess quantitative changes in Cidea and Cidec protein expression due to the technical limitations of the antibodies available. Furthermore, we also did not examine any knockdown or knock-in in vivo or in vitro models to prove whether Cidea and/or Cidec are the indispensable contributors in this specified animal model. However, a recent finding showed that the partial silencing of the expression of Cide genes improved fatty liver in a mouse model47. In the study of histone modifications around Cidea and Cidec, we performed multiple qPCRs after immunoprecipitation, and due to our technical limitations, we were unable to conduct next-generation sequencing. Therefore, the positions of specific modifications were not mapped.

In conclusion, using the present experimental model, we demonstrated plasticity in developmentally-programmed histone modifications around the Cidea and Cidec genes during the amelioration of hepatic steatosis by a treatment with TUDCA. The present study has provided insights into the therapeutic targets of histone modifications for the future of precision medicine for developmentally-programmed hepatic steatosis.

Materials and Methods

Animal model

The mouse model of UN in utero has been well established in our previous studies12,13,14,15,16. Different time points of 3 independent cohorts were described in Fig. 6, as previously reported12. Pregnant C57Bl/6NCr mice at 7.5 days post-conception (dpc) were purchased from Japan SLC, Inc. (Hamamatsu, Japan) and fed a regular chow diet (formula number D06121301, Research Diets Inc., New Brunswick, NJ) under a 12-h light/dark cycle. On 11.5 dpc, dams were randomly divided into an ad libitum or normal nourished group (group NN dams; n = 20) and a caloric restriction or undernourished group (group UN dams; n = 20). Group NN dams were fed powdered regular chow ad libitum and group UN dams were fed a 40% calorie reduction i.e. 60% of the average daily intake of group NN dams, from 11.5 dpc to the day before the delivery of pups (18.5 dpc). We used male offspring only and adjusted the number to 8 pups per litter and cross-fostered (on 1.5 days of age) them with NN dams (up to 23.5 days of age)12. Pups were then fed the regular chow diet for 1 week (up to 9 weeks) followed by a high-fat diet (HFD) containing 60% lipids (formula number D12492, Research Diets Inc.) up to 22 weeks in order to mimic an obesogenic diet, as previously described12. We initiated different experimental procedures with 3 cohorts as follows: cohort 1 (at 9 weeks; without HFD), cohort 2 (at 17 weeks; with HFD), and cohort 3 (at 22 weeks; HFD with or without the TUDCA treatment)12. In cohort 3, TUDCA (Merck Japan Ltd., Tokyo, Japan) was orally administered by gastric lavage at 0.5 g/kg body weight per day while pups were on HFD, according to our previous study12 (Fig. 6). All animal procedures were approved by the Institutional Animal Research Committee of Hamamatsu University School of Medicine (H20–014) and conducted in accordance with the ARRIVE guidelines and the standards of humane animal care by the criteria outlined in the “Guide for the Care and Use of Laboratory Animals”.

Figure 6
figure 6

Schematic presentation of the undernourishment in utero mouse model: male offspring under an obesogenic diet. Cohort 1: After weaning, NN pups from ad libitum NN dams and UN pups from calorie-restricted UN dams were fed a regular chow diet up to 9 weeks. Cohort 2: High-fat diet (HFD) between 9 and 17 weeks. Cohort 3: Pups were randomly subdivided into vehicle groups (NN-Veh and UN-Veh) and TUDCA-treated groups (NN-TU and UN-TU) still under HFD up to 22 weeks.

Blood and tissue sampling

The time points of sampling are shown in Fig. 6 as previously described12. Eight randomly selected pups per group (cohorts 1, 2, and 3) were used in subsequent analyses, except for the microarray analysis of cohort 1 (3 randomly selected pups per group) and cohorts 2 and 3 (4 randomly selected pups per group) as well as DNA MBD sequencing (4 randomly selected pups per group in cohort 3). Total lipids were extracted from liver tissue and measured as previously described12.

Microarray analysis

Aliquots (100 ng) of total RNA obtained from 3 animals per group on week 9 (cohort 1) and 4 animals per group on weeks 17 (cohort 2) and 22 (cohort 3) were individually converted to cRNA and labeled with a Gene ChIP® Poly-A RNA Control Kit, WT Amplification Kit, and Gene ChIP® WT Terminal Labeling Kit (Affymetrix, Santa Clara, CA) according to the manufacturer’s instructions. Hybridization, washing, and staining were performed using Affymetrix® MoGene2.1 ST Array Strips and a GeneAtlas® Hybridization Wash and Stain Kit for WT Assay Strips (Affymetrix), according to the manufacturer’s protocols. After washing, MoGene2.1 Array Strips were analyzed using a GeneAtlas Imaging Station (Affymetrix). Data analyses were performed using Expression Console (Affymetrix) and Transcriptome Analysis Console (Affymetrix). The cut-off point: ≤−2 or ≥2 of a linear fold change and ANOVA P values were used as described in our previous studies45,46.

Relative quantitative real-time RT-PCR analysis

The gene expression levels of Cidea and Cidec were measured by relative quantitation in real-time RT-PCR using the Thunderbird SYBR qPCR Mix (TOYOBO, Osaka, Japan). The primers used were purchased from Sigma-Aldrich Japan and were as follows: Cidea; forward primer: CTGTAGCTGTGCCCTGGTTA, reverse primer: CGGGACAGTTCCTGGTCTAT and Cidec; forward primer: ATTCTGAGTCACCCAGGCC, reverse primer: AAATGAGAACAAGAGAGGCAGC.

Histopathological analysis

Immunohistochemistry was performed as previously described16 using the primary antibodies of Cidea (1:400; Novus Biologicals, U.S.A.) and Cidec (1:400; Bioss Antibodies Inc., Massachusetts, U.S.A.). Frozen sections were used for Oil Red O staining. We then measured the mean size (µm²) of lipid droplets (LD) and the mean diameter of each LD using WinRoof ver7.4 software at a scale of a 50-µm magnification. Measurements of 20 randomly selected LD in 4 areas of interest per slide near the central vein were performed in 3.88-mm² fields.

DNA MBD sequencing

MBD sequencing was performed on cohort 3 only (NN-Veh, NN-TU, UN-Veh, and UN-TU; n = 4 each)48. In brief, DNA extracts were sheared to an average length of 200 bp using Covaris S2. Fragment distribution was assessed by the Agilent 2100 Bioanalyzer. The MethylCAp kit (Diagenode, Belgium) was used to enrich methylated fragments from sheared DNA. A library was prepared using the NEBNext® Ultra™ DNA Library Prep Kit for Illumina® (New England Biolabs, Frankfurt am Main, Germany) combined with the NEBNext® Multiplex Oligos for Illumina®-Index Primers Set 1 (New England Biolabs, Frankfurt am Main, Germany). Fragments of approximately 300–450 bp were excised and purified. Paired-end sequencing was performed on Illumina HiSeq4000 PE50 (2×50 bp) and reads were mapped using Bowtie 2 (v2.1.0) software in the sensitive mode. Only concordantly and paired mapped reads located within 400 bp of each other on the mouse reference genome build GRCm38 (mm10) were retained using Bowtie 0.12.7 and peaks were called using MACS14 (v1.4.2). Sites with potential differential methylation were analyzed using DiffBind R-package v2.0.9.

ChIP assay and quantitative real-time PCR

A ChIP assay was performed on cohort 3 (n = 8 for the 4 groups). Liver tissue samples were homogenized using fixative buffer (1% formaldehyde, 4.5 mM Hepes, pH 8.0, 9 mM NaCl, 0.09 mM EDTA, and 0.04 mM EGTA in phosphate-buffered saline49) and incubated at 37 °C for 30 min. The fixation reaction was terminated by the addition of glycine to a final concentration of 1.5 M. After being washed in fluorescence-activated cell sorting solution (2% bovine serum and 0.05% NaN3 in PBS), samples were sonicated in SDS lysis buffer (16 mM Tris/HCl, pH 8.0, 10 mM EDTA, 1% SDS, and complete mini) to yield DNA fragments ranging in size between 200 and 500 bp, as confirmed by electrophoresis using a 2% agarose gel. Protein concentrations were adjusted to 2 mg/mL after sonication. ChIP assays were performed using 400 µg of protein in SDS-lysis buffer50, and 2 µg each of specific antibodies (Supplementary Table S7). Extracted DNA was then subjected to quantitative PCR using primers corresponding to the indicated sites in the enhancer/promoter and transcribed regions (−2000 bp, −1000 bp, −500 bp, +1 bp, +1000 bp, +2000 bp, +5000 bp, and over +5000 bp) of target genes. Primers for the ChIP assay are listed in Supplementary Table S8. All ChIP signals were calculated by the 2-DDCt method and were normalized to the corresponding input signals as follows: ChIP signals (% input) ¼ 2 (Ct of the input–Ct of the IP sample) × 100, where IP is immunoprecipitation. Notably, the IP of two groups was performed on different days. Non-specific antibody binding to protein-DNA fragments, as indicated by the percentage input using IgG, was negligible (IgG signals <0.1).

Statistical analysis

Data are expressed as means ± standard deviations (SDs). The significance of differences between two mean values was assessed using the Student’s t-test or Mann-Whitney U test as appropriate. The significance of differences among more than four mean values was assessed with the Steel-Dwass test or a one-way ANOVA (unpaired). A p value of less than 0.05 was regarded as significant.