1 Introduction

Nonalcoholic fatty liver disease (NAFLD) is a metabolic disease characterized by excessive liver fat accumulation (fat content exceeds 5% of liver weight) and inflammation (1–3). NAFLD, associated with underlying metabolic syndrome, including type 2 diabetes, high blood pressure, elevated triglycerides and cholesterol, was recently renamed “Metabolic dysfunction-Associated Fatty Liver Disease” (MAFLD) (4). The worldwide prevalence of NAFLD is increasing annually, with recent (2020) rates of 25% among individuals aged 2-70+ (5–8). NAFLD alone can have severe health consequences, and when left untreated, NAFLD can progress into more severe end-stage liver disease such as hepatocellular carcinoma (HCC). HCC is one of the leading causes of cancer-related death worldwide (6, 7). Therefore, it is critical to elucidate the early causes and mechanisms underlying NAFLD.

The most marked characteristic of NAFLD is liver fat accumulation caused by metabolic dyshomeostasis (8), whereby energy production is disrupted in favor of energy storage. While metabolic dyshomeostasis can occur at any point in time, in utero development has long been recognized as a sensitized window of exposure (9) when offspring are particularly vulnerable to nutrient fluctuations (10, 11) that can alter this balance. For example, depletion of vitamin D levels during in utero development in Wistar rat models caused severe liver fat accumulation in adulthood (12). This effect was attributed to underlying mitochondrial dysfunction and abnormal liver lipid metabolism (12). Sprague-Dawley rat models also showed that maternal vitamin D deficiency during pregnancy is associated with offspring insulin resistance, a key regulator of energy storage in the liver (13). In this model, the insulin resistance was attributed to deficiency-induced upregulation of inflammatory cytokines in the liver and serum (13).

Chronic inflammation is another key characteristic of NAFLD (3), which when untreated left can progress to fibrosis and HCC (14, 15). In Sprague-Dawley rat models of developmental vitamin D deficiency (DVD), DVD induced elevated liver and pancreatic oxidative stress levels in adult offspring despite vitamin D repletion (16). Following a similar exposure model, Masako et al. found that DVD in C57BL/6J mouse models caused permanent changes in the proportions of inflammatory cells in the adult liver and expression of genes regulating lipid metabolism (17). One hypothesis is that vitamin D deficiency during in utero development creates an inflammatory state in the fetus. Supporting this hypothesis, vitamin D has been shown to regulate placental inflammation during human pregnancy (18) and in C57BL/6J mice, developmental vitamin D deficiency promotes infiltration and activation of liver macrophages (17). However, the role of DVD in NAFLD-related liver inflammation has not been directly investigated.

Severe maternal vitamin D deficiency in humans during pregnancy impairs fetal growth and development, causing low birth weight and small for gestational age (SGA) (19, 20). Both low birth weight and SGA are recognized risk factors for NAFLD, especially under conditions of rapid postnatal weight gain (14, 15). This is in part due to underlying disruptions in metabolic-endocrine maintenance (21, 22), making the offspring more susceptible to liver disease outcomes.

Previously, we found that DVD exposure in Collaborative Cross (CC) mice caused increased adult body weight and fat mass despite full recovery of vitamin D sufficiency in adulthood (23). Furthermore, a comparison of offspring from reciprocal crosses between strains CC011/Unc (CC011) and CC001/Unc (CC001) showed that this effect differed based on the parental origin of the genome (POG) (23). F1 males from CC001-dams x CC011-sires were susceptible to DVD-induced increased adult adiposity, while F1 males from CC011-dams x CC001-sires were resistant to DVD-induced increased adult adiposity (23).

Here, we leveraged this unique POG model in an ancillary study to investigate a novel role for vitamin D in the developmental origins of liver metabolic dysfunction. Using liver samples collected from our previously published DVD model (23), we defined DVD-induced changes in liver transcriptional pathways that regulate energy metabolism, cholesterol biosynthesis, inflammation, growth & development, and liver detoxification; examined evidence for a role in perturbing liver metabolic processes; and investigated the role of epigenetic mechanisms in the persistence of these effects. In addition, the use of our Collaborative Cross POG model allowed us to define interindividual differences in susceptibility to DVD-induced adult liver disease.

2 Materials and methods2.1 Animal husbandry, dietary treatment, and breeding

All animals were handled in accordance with the Guide for the Care and Use of Laboratory Animals under an approved animal use protocol at the University of North Carolina (UNC) at Chapel Hill. As previously published (23), Collaborative Cross inbred mouse strains, CC001/Unc (CC001) and CC011/Unc (CC011), were purchased from the UNC Systems Genetics Core Facility (Chapel Hill, NC) (24). Inbred CC001 and CC011 virgin dams, aged 4-6 weeks, were placed on one of two isocaloric purified diets for 5 weeks prior to the start of breeding, which is the timeline shown to induce deficiency in mice (25): VDS (vitamin D sufficient diet, 1000 IU/kg vitamin D3, AIN-93G, 110700, Dyets Inc., PA) or VDD (vitamin D depleted diet, 0 IU/kg vitamin D3, AIN-93G, #119266, Dyets Inc., PA). Sires stayed on the VDS diet except when mating with VDD dams (14 days). Dams remained on diet during and after reciprocal mating to CC001 and CC011 sires to generate male F1 offspring that were primarily genetically identical other than different parental origin of the genome (POG) and different mitochondrial and Y chromosomes. (CC001-dam xCC011-sire)F1 and (CC011-dam x CC001-sire)F1 males were maternally exposed to either VDS or VDD diets from conception to weaning. At PND21, mice were weaned onto standard rodent chow diet (2400 IU/kg vitamin D3; Teklad diet #8604, Harlan Laboratories, Germany) and remained on this vitamin D sufficient diet for 5-6 weeks, which is a timeline shown to be sufficient for vitamin D repletion in mice (25). At 8-9 weeks of age, all F1 adult male offspring were euthanized by CO2, and livers were collected and flash-frozen in liquid nitrogen and stored at −80°C. Throughout the study, all mice were housed and maintained at a vivarium temperature of 21-23°C with a 12-h light cycle and ad libitum access to sterilized water and rodent chow (23).

2.2 Sample selection (total RNA-sequencing, metabolomics, bisulfite-sequencing)

For RNA-seq, six adult male F1 offspring liver samples were selected from at least three different litters for each diet and POG group (n=6/diet/POG, 24 samples total). A subset of these samples (n=3/diet/POG from 3 litters, 12 samples total) were used for the metabolomics and bisulfite-seq experiments. Whole liver samples were pulverized and mixed while frozen and then aliquoted while frozen for each experiment.

2.3 Total RNA-Seq

Total RNA was isolated from adult male F1 offspring livers using Trizol reagent following the manufacturer’s protocol (#155960926, Life Technologies, NY). All remaining sample prep and sequencing were performed by the UNC High-Throughput Sequencing Facility (HTSF) at UNC Chapel Hill. RNA was prepped for RNA-Seq using the RNA Clean & concentratorTM-5 kit according to the manufacturer’s protocol (#R1013, Zymo Research, CA). All RNA analytes were assayed for RNA integrity, concentration, and fragment size. Samples for total RNA-Seq were quantified using RNA Qubit (#Q32855, Invitrogen, MA). All samples had a RIN>7.0 (average RIN = 8.9). 500 ng of input RNA was used for library preparation. Fragmentation was performed for 6 minutes at 85°C following the KAPA Stranded RNA-Seq Kit with RiboErase protocol (#8098131702, Illumina Inc., CA). Indexed libraries were prepared and run on NovaSeq XP (#20043131, Illumina Inc., CA) to generate an average of 221 million paired-end reads (100 bp) per sample library. The read length was 2 x 100bp.

Raw sequencing reads were filtered for a quality score > 20 in at least 90% of bases using fastq_quality_filter (version 0.0.14) (26). Sequence adapters were then trimmed using Cutadapt (27) (version 1.12). We then adopted a 2-pass alignment approach utilizing CC001 and CC011 pseudo-references (retrieved from http://www.csbio.unc.edu/CCstatus/index.py?run=Pseudo) to obtain and quantify strain-specific alignments. All samples were first aligned to a hybrid reference genome in which genomic sequences for CC001 and CC011 (based on GRCm37/mm9 coordinates) were concatenated (28), and each chromosome name (chrM) was prepended with the appropriate CC strain. Alignments were performed using STAR (29) (version 2.7.3a), supplying GENCODE vM25 gene annotations as a guide. The splice junction output files were then concatenated, chrM was removed, and then one sample was re-aligned to the initial hybrid reference to insert these newly detected splice junctions and create the final 2-pass hybrid reference genome. All samples were then re-processed onto this 2-pass reference to generate a final set of alignments. Transcript abundance for each strain-gene was then estimated using salmon (30) (version 1.6), and differential expression was detected using DESeq2 (31) (version 1.26.0), where strain and isoform were aggregated by gene symbol.

Differentially expressed genes were identified in two ways using DESeq2 (31) (version 1.26.0): (1) Main diet effects independent of POG - samples from both POGs were combined and data was modeled for diet effects adjusted for POG (~ POG + diet) (Supplementary Table 1A); and (2) POG-specific diet effects (Supplementary Tables 1B, C) - samples were stratified by POG and modeled for diet effects on each POG separately.

Differentially expressed genes (DEGs) with a nominal p<0.05 were imported into PANTHER (32). DEGs that mapped to the PANTHER mus musculus reference genome (version 2021_03) were assigned to pathways and then tested for pathway overrepresentation using a Fisher’s exact test (PANTHER (32, 33), version 17.0, released 2022-10-13) with FDR<0.1. Gene enrichment analyses were performed on DEGs (nominal p<0.05) using the canonical pathways function in Ingenuity Pathway Analysis (34) (IPA, version 01-20-04 (01-20-04)) with FDR<0.05 (Benjamini and Hochberg method).

Two-way hierarchical clustering of log2 transformed VST-normalized gene expression values was performed in JMP Pro (version 17.0.0 (623769), SAS Institute Inc., Cary, NC, 1989–2021). Gene functions and biological processes were determined using GeneCards (35) (version 5.11.0, build 656) and UniProt (36) (release 2023_01).

2.4 Metabolomics

Untargeted profiling of whole liver metabolites was performed using mass spectrometry (Metabolon, Inc (37).) as previously described (38). In brief, four distinct mass spectrometry platforms were utilized in order to capture metabolites across a wide range of varying chemical properties: 1) RP/UPLC-MS/MS with acidic positive ion mode electrospray ionization (ESI) for hydrophilic compounds, 2) RP/UPLC-MS/MS with acidic positive ion mode ESI for hydrophobic compounds, 3) RP/UPLC- MS/MS with basic negative ion mode ESI, and 4) HILIC/UPLC-MS/MS with negative ion mode ESI (38). Metabolites were identified based on their spectrum properties, which included criteria such as m/z and retention time/index. All data underwent extensive quality control by Metabolon, Inc. to account for sample and instrument variability, including interday differences in instrument tuning. To account for these differences while also maintaining appropriate sample variation, raw data were median normalized (38). Samples with measurements that fell below the limit of quantitation were imputed to the limit of quantitation.

Diet and POG effects on metabolite concentrations were analyzed by Principal Component Analysis (PCA) and linear regression modeling. PCA was performed on the normalized and imputed metabolite concentrations using the R package “ropls (ver 1.16.0).” Significant differences in relative concentrations of metabolites were identified in two ways using JMP Pro (version 17.0.0 (623769), SAS Institute Inc., Cary, NC, 1989–2021): (1) Main diet effects independent of POG - samples from both POGs were combined and data was modeled for diet effects adjusted for POG (~ POG + diet); and (2) POG x diet interactive effects were modeled (~ POG + diet + “POG x diet”). Model assumptions for normality (Shapiro-Wilk test) and homoscedasticity (Bartlett’s test) were tested for each model, and metabolites that failed either test (p<0.05) were transformed (log10 or Box-Cox) to best meet assumptions (Supplementary Tables 1D, E). In total, after transformation, we performed effect testing on 643 metabolites, including 35 metabolites that did not meet assumptions in the main diet model (1 with significant main diet effect, p<0.05) and 50 metabolites that did not meet assumptions in the interactive model (7 with significant interactive effects, p<0.05). All significant metabolites were reported, but metabolites that did not meet assumptions were not included in the pathway analyses. For pathway analyses, metabolites that met assumptions for linear regression and had significant diet or POG x diet effects (nominal p<0.05) were imported into Metaboanalyst (pathway enrichment function, version 5.0) (39, 40), and pathways with an FDR<0.1 were considered significant. Twenty-five metabolites, predominantly classified as lipids, had no discernable Metaboanalyst-matched name and could not be included in pathway analyses.

2.5 Bisulfite-seq

Genomic DNA was isolated using phenol-chloroform extraction as previously described (23). Quantity and quality of DNA were assessed by NanoDrop 2000 spectrometer (Thermo Scientific, DE) and by the Quant-it PicoGreen dsDNA assay (#P7589, Life Technologies, NY). 3 μg dsDNA was used as input, and library preparation was carried out using the SureSelectXT Mouse Methyl-Seq system (Agilent Technologies, CA) following the manufacturer’s protocol. Libraries were multiplexed and sequenced on the Illumina HiSeq2500 (Illumina Inc., CA) at the David H. Murdock Research Institute (DHMRI). Illumina HiSeq single-end 100 bp long reads were generated and trimmed as previously described (41), and the resulting reads were mapped to the UCSC (42) GRCm38/mm10 genome. Only CpGs with ≥10X coverage and a quality score ≥20 were included in the analyses. Single Nucleotide Polymorphisms (SNPs) and other sequences that would potentially interfere with CpG methylation calls were removed from the dataset. ‘Blacklisted’ regions (42) were also removed from the dataset as previously described (41).

MethylKit (version 0.9.5) (43) was used to identify differentially methylated CpGs (DMCs) by logistic regression, as previously described (41). Regression analyses were performed on data after normalizing the coverage between loci using a scaling factor based on the differences in median coverage distributions among CpGs. P-values from the regression analyses were adjusted for multiple comparisons using the SLIM method (FDR<0.05 cutoff) within the MethyKit (43) package. A false discovery rate threshold of 0.05 was applied to identify DMCs. Gene annotations for DMCs (q<0.05) were determined using the UCSC table browser tool (GRCm38/mm10 genome) (44) (Supplementary Tables 1F, G).

Differentially methylated genes (DMGs) with a q<0.05 were imported into PANTHER (32). DMGs that mapped to the PANTHER mus musculus reference genome (version 2021_03) were assigned to pathways and then tested for pathway overrepresentation using a Fisher’s exact test (PANTHER (32, 33), version 17.0, released 2022-10-13) with FDR<0.1.

3 Results3.1 Developmental vitamin D deficiency (DVD) causes persistent transcriptional dysregulation of liver cholesterol biosynthesis, energy metabolism, growth & development, inflammation, and detoxification pathways

To assess the role of DVD in transcriptional programming of molecular pathways of liver disease, we performed an ancillary study on a previously phenotypically characterized DVD cohort of mice (23). Total RNA-seq was performed on livers from 8-week-old adult male offspring (n=6/diet/POG, 24 samples total) treated with either a vitamin D depleted (VDD) or vitamin D sufficient (VDS) diet only from conception until weaning then repleted until adulthood (Figure 1A). Adult offspring were assessed for persistent effects of DVD on the liver transcriptome. Bodyweight, fat mass, and vitamin D exposure levels have been previously published for these mice (23). This model incorporates offspring from reciprocal crosses that are genetically identical except for mitochondria and sex chromosomes [POG1-(CC011xCC001)F1 and POG2-(CC001xCC011)F1 males] to characterize interindividual differences related to parental origin of the genome (POG) (Figure 1B).

Figure 1 Dietary treatment scheme and POG model. (A) Dams from both CC001 and CC011 strains were simultaneously treated with one of two diets (vitamin D sufficient (VDS) or vitamin D depleted (VDD)) for 5 weeks prior to the start of breeding and throughout breeding, gestation, and lactation. All offspring were weaned onto standard rodent chow and remained until adulthood (8-9 weeks of age) when they were euthanized and livers were collected. Figure created with BioRender.com (B) A reciprocal cross approach was used to generate F1 offspring (POG1 (CC011xCC001) & POG2 (CC001xCC011)) that differed genetically only for mitochondrial (mt) and sex chromosomes (csbio.unc.edu).

To investigate the effect of DVD independent of POG, we combined POG1 and POG2 datasets and adjusted for POG using linear regression analyses. Of the 53,569 genes queried across the liver genome, 29,471 had detectable expression levels above baseline, and 1,685 were identified as DVD-induced differentially expressed genes (DEGs, nominal p<0.05). Three pathways were significantly overrepresented and downregulated by DVD (PANTHER (32, 33), FDR<0.1): cholesterol biosynthesis (mevalonate pathway), energy metabolism (pentose phosphate pathway - converts glucose to pyruvate), and inflammation (20S proteasome subunits that regulate protein deubiquitination and are triggered by oxidative stress (45)) (Figure 2A, Supplementary Table 2). Despite relatively small fold changes in gene expression, we observed consistent downregulation of DEGs within these pathways (Supplementary Table 2).

Figure 2 Overrepresented and enriched gene pathways transcriptionally dysregulated by DVD independent of POG. Differentially expressed genes (DEGs) contributing to these pathways are listed in Supplementary Table 1A. (A) Gene pathway overrepresentation analysis (PANTHER) on 1,338 PANTHER-annotated DEGs with p<0.05 (detected after adjustment for POG). Pathways with FDR<0.1 are shown and classified by biological processes. PANTHER pathway ID numbers are shown in parentheses. (B, C) Gene pathway enrichment analysis (IPA) on 1,549 IPA-annotated DEGs with p<0.05 (detected after adjustment for POG). (B) Pathways with -log(adjp-value)<0.05 (after Benjamini and Hochberg (BH) correction) are shown and classified by biological process (colored boxes). Pathways are ordered from lowest (bottom) to highest (top) –log(adjp-value). Bubble size indicates the ratio of overlapping genes to the total number of genes within that pathway. Numbers overlaying each pathway bubble represent the number of DVD-altered genes within that pathway. Bubbles are colored by the z-score predictor of pathway activation or inactivation: activated pathways (light to dark orange), unchanged pathways (grey), and deactivated pathways (light to dark blue). (C) The 32 DVD-induced DEGs enriched in the oxidative phosphorylation pathway and their functional locations in each of the five electron transport chain complexes. Gene names are shown in green bubbles with the direction of change indicated (green=downregulated; red=upregulated). Data were analyzed through the use of IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis).

We expanded this finding using Ingenuity Pathway Analyses (IPA) to show the enrichment of DEGs in 39 mostly downregulated canonical pathways (FDR<0.05) that regulate cholesterol biosynthesis, energy metabolism, growth & development, inflammation, and liver detoxification (Figure 2B). Most notably, this included: (i) predicted downregulation of cholesterol biosynthesis; EIF2 and mTOR signaling, which can inhibit protein synthesis and respond to cellular stress (46); the NRF2-mediated oxidative stress response and subsequent glutathione redox reactions and glutathione-mediated detoxification; and xenobiotic sensing pathways regulated by the nuclear receptors CAR, PXR, and AHR; and (ii) predicted upregulation of the sirtuin signaling pathway, which regulates hepatic lipid and glucose metabolism, stress response, DNA repair mechanisms, and inflammation (47) (Figure 2B). The xenobiotic sensing pathways included differentially expressed glutathione-s-transferases (Gsts), which are key phase II liver detoxification enzymes. We also found enrichment of 32 downregulated DEGs in oxidative phosphorylation across all five electron transport chain (ETC) complexes (Figure 2C). This finding is consistent with VDR knockdown cell culture models that showed reduced ATP production and oxidative phosphorylation capacity (48–50). Unsupervised hierarchical clustering of the expression profiles for these pathways showed that the differences in gene expression mostly clustered by diet (Supplementary Figures 1A-H). However, although both strains exhibited similar responses to DVD on average, the extent of heterogeneity in response differed by strain. For example, across biological replicates, we observed more consistent diet profiles in POG1 for genes regulating oxidative phosphorylation, xenobiotic metabolism, and inflammation (including LPS/IL-1 mediated inflammation, neutrophil trap signaling, granzyme A, and protein ubiquitination). Meanwhile, POG2 had a more consistent diet profile for cholesterol biosynthesis and EIF2 signaling (Supplementary Figures 1A-H).

Seven DEGs (including 2 predicted VDR targets, Mkrn2 and Trip13) (51) remained significant after correction for multiple testing (FDR<0.1) and were similarly affected for both POGs (Supplementary Figure 2A). These included growth & development and inflammation genes involved in rRNA processing and DNA repair, protein trafficking, and chromosome recombination (downregulated by DVD) (35), and genes involved in protein modification, flagellated motility, and cell cycle regulation and apoptosis (upregulated by DVD) (35) (Supplementary Figures 2A, B).

3.2 DVD induces unique transcriptional responses on different genomic backgrounds

To define the effect of DVD that differed between POG, we analyzed the liver RNA-seq data after stratifying by POG. As expected, the combined POG analyses detected the greatest overall transcriptional response to DVD (1,685 DEGs, Figure 3A). However, it failed to find a substantial number of DVD-induced transcriptional effects that were unique to each POG (Figure 3A). POG1 mice exhibited 1,454 DEGs (760 new), while POG2 mice had only half as many (766, 400 new) (Figure 3A). In addition, each POG exhibited a distinct transcriptional response, with only 13% (101) of the POG2 DEGs overlapping with POG1 (Figure 3A).

Figure 3 POG-specific enriched gene pathways transcriptionally dysregulated by DVD. (A) Venn diagram showing the overlap of DEGs across all three datasets (1) DVD effects in POG1, (2) DVD effects in POG2, and (3) DVD effects for both POGs combined (adjusted for POG). Number of DEGs before (nominal p<0.05) and after (FDR<0.1) correction for multiple testing. (B) Gene pathway enrichment analysis (IPA) on 1,360 IPA-annotated DEGs (p<0.05) detected in POG1. (C) Gene pathway enrichment analysis (IPA) on 662 IPA annotated DEGs (p<0.05) detected in POG2. (B, C) Enriched pathways (-log(adjp-value)<0.05) are shown and classified by biological process (colored boxes). Bubble size indicates the ratio of overlapping genes to the total number of genes within that pathway. Numbers overlaying each pathway bubble represent the number of DVD-altered genes within that pathway. Bubbles are colored by the z-score predictor of pathway activation or inactivation: activated pathways (light to dark orange), unchanged pathways (grey), and deactivated pathways (light to dark blue). Single asterisk (*) indicates pathways also enriched in the DVD independent of POG dataset (Figure 2). Double asterisks (**) indicate pathways enriched in all three datasets. DEGs contributing to these pathways are listed in Supplementary Tables 1B, C. Data were analyzed through the use of IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis).

Despite a greater number of DEGs, POG1 had no significantly overrepresented pathways. POG2 had a significant overrepresentation (FDR<0.1) of the mevalonate pathway of cholesterol biosynthesis, which was previously detected in the combined POG analyses (Supplementary Table 3).

DEGs were significantly enriched (FDR<0.1, IPA (34)) in eight pathways for POG1 (Figure 3B) and in almost twice as many (13 pathways) for POG2 (Figure 3C). Although enrichment and downregulation of cholesterol biosynthesis were common to both POGs, POG2 had a greater number of relevant cholesterol biosynthesis DEGs, including enrichment of zymosterol biosynthesis upstream of cholesterol (Figure 3C). POG1 exhibited unique enrichment of downregulated genes involved in pathways that regulate lipid metabolism for energy production (triacylglycerol degradation and stearate biosynthesis), growth & development (reelin signaling), and inflammation (the LPS/IL-1 mediated inhibition of RXR function pathway and a SPINK1 cancer pathway) (Figure 3B). POG1 mice also exhibited unique enrichment of upregulated genes involved in acyl-CoA hydrolysis and γ-glutamyl cycle, which are energy metabolism pathways that increase acetyl-CoA and amino acid availability, respectively (Figure 3B). In contrast, POG2 mice exhibited unique enrichment of downregulated genes involved in inflammation (NRF2-mediated oxidative stress) and liver detoxification (glutathione-mediated detoxification) (Figure 3C).

For POG1, nine DEGs remained significant after correction for multiple testing (FDR<0.1, Supplementary Figure 3A): Mkrn2 and Commd10 (previously detected), and seven genes that regulate protein modification and metabolism, transcription, membrane trafficking, and inflammation (Supplementary Figures 3A-C). For POG2, only two DVD-induced DEGs remained significant after correction for multiple testing (FDR<0.1): Ptar1 (protein prenylation, previously identified) and Slc7a15 (amino acid transport) (Supplementary Figures 4A-C).

3.3 DVD leads to persistent disruption of liver metabolic processes regulating cholesterol biosynthesis and energy metabolism

We used untargeted metabolomics profiling on livers from a subset of transcriptionally profiled samples (n=3/diet/POG from 3 litters, 12 samples total) to investigate the impact of DVD-induced transcriptional changes on the liver metabolic processes. Principal component analysis (PCA) of the 651 metabolites detected revealed a POG-specific diet effect. PC1 (25%) primarily implicated an effect of DVD on POG1. PC2 (16%), in part, implicated an effect of DVD on POG1 and POG2 (Figure 4A). Linear regression analyses to identify the main diet effects (effects of DVD after adjustment for POG) identified 82 metabolites (nominal p<0.05, 0 significant with FDR<0.1) that exhibited mostly reduced abundance on both backgrounds and classified primarily as lipids (~67%), carbohydrates (~16%), and amino acids (~11%) (Figure 4B). These metabolites were not significantly enriched for any pathways (Metaboanalyst (39, 40)). However, consistent with the downregulation of cholesterol biosynthesis and pentose phosphate pathway genes in the liver, DVD-treated mice exhibited decreased (nominal p<0.05) levels of liver cholesterol metabolites (Figure 4C) and pentose phosphate pathway metabolites (Figure 4D).

Figure 4 Metabolic pathways dysregulated by DVD independent of POG. (A) Principal Component Analysis (PCA) run on all 654 metabolites after normalization and data imputation to the limit of quantitation. Top two principal components (PC1 & PC2) are shown with the percent of variation in metabolite levels they each explain. Each dot in the graph represents an individual mouse liver and colors indicate the diet and POG of that sample. (B) Bar graph showing the direction of change for 82 metabolites with significant DVD-induced changes (p<0.05) after adjustment for POG. Bar graphs represent the number of metabolites exhibiting up- or downregulation within each metabolite class and function calculated separately for both POGs combined, for POG1 only, and for POG2 only. Metabolites contributing to these pathways are listed in Supplementary Table 1D. (C) Metabolite levels for four sterols with known roles in cholesterol biosynthesis are plotted for VDS and DVD groups. POG is indicated by open (POG1) and closed (POG2) circles. Pathway illustration shows the differentially regulated cholesterol metabolites and genes detected within our datasets. Blue indicates downregulation of gene expression by DVD. Red indicates upregulation of gene expression by DVD. (D) Metabolite levels for four carbohydrates with known roles in the pentose phosphate pathway are plotted for VDS and DVD groups. POG is indicated by open (POG1) and closed (POG2) circles. Pathway illustration shows the differentially regulated pentose phosphate metabolites and genes detected within our datasets. Blue indicates downregulation of gene expression by DVD. Red indicates upregulation of gene expression by DVD. Illustrations were created with BioRender.com.

3.4 DVD disrupts distinct metabolic processes on each genomic background

Based on the POG-specific diet effect observed in the PCA analysis (Figure 4A), we used an interactive diet x POG linear regression model to identify DVD effects that differed between POGs. We identified 94 metabolites with significant DVD x POG interactive effects (nominal p<0.05, 16 with FDR<0.1 (Supplementary Table 4) that were mainly classified as amino acids (~35%), lipids (~33%), and nucleotide intermediates (13%) (Figure 5). Strikingly, the two POGs exhibited inverse responses to DVD for most of these metabolites (Figure 5). POG1 exhibited increased concentrations of macronutrient metabolites (amino acids, lipids, carbohydrates, etc.) indicative of increased synthesis and increased concentrations of micronutrient cofactors of energy production (e.g., for vitamin B5 & B6), while POG2 concentrations were decreased (Figure 5, Supplementary Figures 5A-E). This includes effects on branched-chain amino acids ((BCAAs) L-leucine, L-isoleucine, L-valine), energy cofactors (e.g., FAD, flavin mononucleotide, and vitamin B6 metabolites) and substrates (e.g., pyruvic acid, succinic acid, fructose), and purine and pyrimidine metabolites (Supplementary Figures 5A-E). Metabolites with DVD x POG interactive effects were enriched (FDR<0.1) in 15 pathways (Figures 6A, B, Supplementary Table 5). We observed inverse effects on metabolites in these pathways, with increased metabolite abundance for POG1 and decreased abundance for POG2 (Figure 6B).

Figure 5 POG-specific effects of DVD on metabolite profiles. Bar graph showing the direction of change for 94 metabolites with significant DVD x POG interactive effects (p<0.05). Bars represent the number of metabolites exhibiting up- or downregulation within each metabolite class and function calculated separately for POG1 and POG2. Metabolites contributing to these pathways are listed in Supplementary Table 1E.

Figure 6 POG-specific enriched metabolite pathways dysregulated by DVD. (A) Pathway enrichment analysis (Metaboanalyst) on 71 Metaboanalyst-annotated metabolites with significant DVD x POG interactive effects (p<0.05). Single asterisk (*) indicates pathways with significant enrichment (-log10(FDR)<0.1). Enrichment ratios indicate the ratio of overlapping metabolites over the total number of metabolites within that pathway. (B) Fold change (DVD/VDS) of metabolites from Metaboanalyst enrichment analysis. X-axis values >1 mean DVD increased abundance, while values <1 mean DVD decreased abundance.

3.5 DVD programs distinct liver DNA methylation profiles on each genomic background that explain only a small portion of the transcriptional changes

We previously showed that DVD causes loss of methylation (LOM) in sperm (41). To determine whether the persistent effects of DVD on liver transcription and metabolic processes could be epigenetically programmed by DNA methylation, we used genome-wide bisulfite-sequencing on the metabolically profiled liver samples (n=3/diet/POG from 3 litters, 12 samples total).

Given the distinct POG response to DVD observed for the transcriptome and metabolome, we identified DVD-induced differentially methylated CpGs (DMCs, q<0.05) for each POG. After removal of CpGs with ≤10X coverage, a quality score ≤20, SNPs, or other sequences that would potentially interfere with methylation calls, 493,989 CpGs were queried across the liver genome. Despite POG1 exhibiting a greater transcriptional response, POG1 exhibited a less robust epigenetic response than POG2 with fewer DMCs (77 DMCs annotated to 40 genes (DMGs)) compared to POG2 (459 DMCs annotated to 260 DMGs) (Figure 7A). For both POGs, DMCs were distributed across the genome (Figure 7B) and primarily LOM (Figures 7A, B). Only one DMC (Chr10:19088363-19088364) and one gene (Prrc2b, different DMCs for each POG) were shared by both POGs (Supplementary Figure 6). To determine whether the methylation profile induced by DVD was similar between the two POGs (despite poor overlap in “significantly” altered genes), we performed two-way hierarchical clustering of DMCs specific to POG1 (Figure 7C) and POG2 (Figure 7D). This analysis confirmed that each POG exhibited a distinct DNA methylation profile induced by DVD. Differentially methylated genes (DMGs) mapped to 10 pathways for POG1 (Supplementary Table 6) and 59 pathways in POG2 with the top 15 listed in Supplementary Table 7. No pathways were overrepresented for POG1. For POG2, we identified significant overrepresentation (FDR<0.1) of growth & development pathways (cadherin signaling pathway and Wnt signaling) (Figure 7E, Supplementary Table 8) with gene methylation profiles showing strong clustering of LOM by diet (Figure 7F).

Figure 7 POG-specific differentially methylated cytosines dysregulated by DVD. Blue indicates loss of methylation (LOM). Red indicates gain of methylation (GOM). (A) Total numbers of differentially methylated CpGs (DMCs, q<0.05) from genome-wide bisulfite sequencing analysis are shown for POG1 and POG2. (B) The distribution of DMCs across chromosomes is shown for POG1 and POG2. (C, D) Two-way hierarchical clustering heat maps generated with (C) All DMCs detected in POG1 (77, q<0.05) and (D) All DMCs detected in POG2 (459, q<0.05). (E) Pathway overrepresentation analysis (PANTHER) for POG2. (F) Two-way hierarchical clustering heat map for DMGs contributing to significant overrepresentation of Cadherin and Wnt signaling pathways in POG2. DMCs and gene annotations are listed in Supplementary Tables 1F, G.

There was little overlap between DMGs and DEGs (Figure 8, Supplementary Table 9). For POG1, only two of the DMGs were also DEGs (Figure 8A, Supplementary Figure 7A). Both were growth & development genes, Igf1r and Slc1a5, that exhibited LOM (Figure 8B). For POG2, there were 13 growth & development DMGs (including two predicted Vdr-targeted genes (Med27, Tfpt) (52)) with significantly altered liver expression: Aust2, Gata2, Hs6st3, Htatip2, Pdzm3, Tfpt (GOM) and Dpysl4, Mast4, Med27, Nrg2, Ped10a, Tial1, Ttn (LOM) (Figures 8C, D). All, except Hs6st3 and Htatip2, had increased expression (Supplementary Figure 7B).

Figure 8 POG-specific overlap between DVD-induced differentially methylated and differentially expressed genes. (A) Venn diagram showing the overlap between DEGs (p<0.05) and DMGs (q<0.05) detected in POG1. (B) Methylation status is shown for differentially methylated DEGs detected in POG1. (C) Venn diagram showing the overlap between DEGs (p<0.05) and DMGs (q<0.05) detected in POG2. (D) Methylation status is shown for differentially methylated DEGs detected in POG2. Blue indicates DEG with LOM. Red indicates DEG with GOM. Single asterisk (*) indicates significant differential methylation (q<0.05). N.S. = Not Significant.

To determine whether the POG-specific effects of DVD on DNA methylation were driven by differences in regulators of de novo methylation (53), we assessed gene expression levels for DNA methyltransferases Dnmt1, Dnmt3a, and Dnmt3b. No significant changes in gene expression were detected for either POG (Supplementary Figures 8A, B). Furthermore, no significant changes were detected amongst folate (54) and methionine-related metabolites (51) involved in the establishment and maintenance of DNA methylation (Supplementary Figure 8C). Individual SAM and SAH levels were not significantly altered for either POG, however, SAM/SAH ratios were significantly reduced by DVD only in POG1.

4 Discussion

In this study, we demonstrate that developmental vitamin D deficiency programs molecular pathways regulating essential liver functions (Figure 9A). This is the first in vivo study to define the persistent liver transcriptional profile induced by DVD, to connect these to actively dysregulated liver metabolic processes (metabolome), and to investigate the role of epigenetic mechanisms (DNA methylome) in driving these persistent effects. We found broad transcriptional suppression of gene pathways that are essential to the liver’s ability to respond to environmental stressors, including xenobiotic metabolism, inflammation, and cholesterol biosynthesis. This implicates a potential mechanism by which DVD exposure increases the risk of liver damage and disease. In support of our previous phenotypic findings (23), the Collaborative Cross POG + DVD model demonstrated interindividual differences in liver response to DVD. Although DVD altered similar molecular pathways for both POGs, they differed substantially in penetrance such that the proportion of responders vs. non-responders differed. This was especially striking for the epigenomic and metabolomic responses, where the latter showed inverse effects of DVD on each POG background. This finding shows that the genomic context controls an individual’s response to DVD at the molecular level, demonstrating that a gene x environment mechanism drives the deleterious effects of DVD.

Figure 9 Proposed model for adult liver pathways programmed by early life exposure to DVD. (A) Key POG-specific DVD-induced gene expression changes are grouped by major pathway classifications. Red arrows indicate downregulation by DVD. Green arrows indicate upregulation by DVD. Illustration created in part by BioRender.com. (B) Key DVD-induced metabolite changes are shown for POG1 (open circles) and POG2 (closed circles) and grouped by metabolic process. Solid arrows indicate a direct relationship between metabolites. Dashed arrows indicate an indirect relationship between metabolites. Red indicates increased by DVD. Blue indicates decreased by DVD. N.S., Not Significant; N.D., Not Detected.

The downregulation of genes controlling cholesterol biosynthesis in the liver was an effect of DVD for both POGs, although penetrance was greater in POG2. Given the essential role of liver-generated cholesterol in cellular structure and function and even vitamin D production, prolonged impaired cholesterol biosynthesis is likely to have deleterious consequences on health. We found subtle but consistent downregulation of genes in the mevalonate pathway of cholesterol biosynthesis and a reduced abundance of liver cholesterol metabolites from this pathway. This effect of DVD is likely a “programmed” effect specific to developmental exposure because we have previously shown that chronic vitamin D deficiency in adult mice from the same strains did not exhibit changes in liver cholesterol levels (38). In humans, vitamin D deficiency has actually been implicated in increased cholesterol (55, 56), but critical longitudinal studies have not been performed to show the time course of deficiency and how cholesterol levels and regulatory pathways respond over time. In fact, prolonged high intracellular cholesterol levels are known to induce mitochondrial dysfunction, leading to inhibition of the mevalonate pathway as a compensatory mechanism for decreasing cholesterol production (57, 58). A similar phenomenon has been demonstrated in VDR KO mice that exhibit decreased serum cholesterol levels at 6 months of age, suggesting that loss of VDR function can also induce attenuation of cholesterol biosynthesis (59). Despite the downregulation of the mevalonate pathway, a downstream target and regulator of protein prenylation, Ptar1, was unexplainably upregulated by DVD. Ptar1 functionality as a subunit of the prenyltransferase, GGTase3, has only recently been recognized and minimally described (60). Therefore, our finding implicates a previously unknown role for Ptar1 function that is independent of GGTase3 that could explain the upregulation of Ptar1 expression in the absence of upregulation of the mevalonate pathway. Regardless, our finding of DVD-induced downregulation of the mevalonate pathway has provided a previously undefined role of vitamin D in developmentally programming the molecular mechanisms regulating cholesterol production and, potentially, subsequent protein prenylation.

We found that DVD also primarily downregulated genes controlling the liver’s response to injury and infection. Downregulation of liver detoxification pathways is likely to weaken the liver’s defense system, making it more susceptible to oxidative stress and less effective at deactivating xenobiotic toxicants, all of which would increase susceptibility to liver damage and disease. Here, the penetrance differed by genomic context, with POG1 exhibiting a more consistent liver detoxification and inflammatory response across DVD-treated mice. Both POGs exhibited downregulation of the xenobiotic sensing pathways, CAR, PXR, and AHR, which included five glutathione-s-transferases (Gst) (liver detoxification genes); however, this effect was more consistently observed amongst POG1 mice. Similarly, POG1 mice exhibited more consistent downregulation of inflammatory pathways, including neutrophil trap signaling, granzyme A, protein ubiquitination, and LPS/IL-1 mediated inhibition of RXR function. POG2 demonstrated a similar response across granzyme A and protein ubiquitination pathways. Downregulation of these pathways could suggest a suppressed inflammatory response system and subsequent impaired response to injury or infection. In contrast, POG1 mice also exhibited potential pro-inflammatory signatures, including increased expression of Saa2, an inflammatory response gene (61, 62); Mkrn2, a negative regulator of inflammation (63, 64); and Spink1, which has been shown to regulate redox homeostasis (50) and is overexpressed in ~70% of HCC patients (65). Two stress-response pathways that regulate protein translation, EIF2 signaling and mTOR signaling (46), were downregulated by DVD across both POGs, with EIF2 signaling showing the most robust response. Programmed downregulation of EIF2 and mTOR signaling could suggest lower protein translation activity, or it could suggest a suppressed response to cellular stress. However, further data is required to differentiate these scenarios.

DVD altered genes controlling liver development and energy availability. For both POGs, DVD decreased intermediates of glycolysis and the pentose phosphate pathway (Figure 9B). However, POG1 demonstrated a unique upregulation of glycerol-3-phosphate, suggesting higher rates of glycerol synthesis via glyceroneogenesis (Figure 9B). Other key distinctions between POG1 and POG2 were defined for lipid and amino acid profiles. Lipid synthesis is an energy-demanding processes fueled by lipolysis (66). In POG1, which exhibited a lean phenotype (23), increased levels of long-chain fatty acids (lipids) were detected and supported by decreased levels of energy substrates (carbohydrates & TCA cycle metabolites) (Figure 9B). Furthermore, POG1 livers had increased levels of aromatic and BCAAs (Figure 9B). These responses were not seen in the adiposity-susceptible POG2 (23), where lipid and amino acid abundance were mainly decreased. One possible explanation for these POG differences in metabolic response could be the upregulation of acyl-CoA hydrolysis and the γ-glutamyl cycle in POG1, which can be upregulated in response to disrupted oxidative phosphorylation as a way to increase cellular availability of acetyl-CoA and amino acid substrates (cysteine, glycine, glutamate), respectively (67). Surprisingly, DVD downregulated oxidative phosphorylation in both POGs, but only POG1 exhibited this compensatory response. Overall, the data indicate that the POG1 response is indicative of increased energy utilization, while the POG2 response indicates decreased energy utilization. Future work is required to determine what characteristic(s) of the parental genome is causing these differences. However, given that most of the changes observed are linked to mitochondrial function, it is likely that genetic differences in the mitochondrial (mt) genes contributed by different founder strains (POG1=NZO/HILtJ; POG2=C57BL/6J) are at least in part responsible.

Our finding that DVD altered the adult liver DNA methylation profile with primarily loss of methylation epimutations supports our previous findings implicating vitamin D in developmental epigenetic programming (23, 41). Strikingly, DVD not only induced a completely different liver DNA methylation profile in each genomic context, but POG2 exhibited substantially greater genomic epigenetic perturbation (~6X as many DMCs and DMGs distributed across the genome). We previously observed a similar pattern of differences in POG epigenetic response at imprinting control regions (ICRs) (23). The overrepresentation of differentially methylated Wnt and Cadherin signaling genes in the liver also mirrors what we found in the germline (sperm) of DVD-treated males from POG2 (41). Surprisingly, the only overlap between DMGs and DEGs