Precision pharmacological reversal of strain-specific diet-induced metabolic syndrome in mice informed by epigenetic and transcriptional regulation

Diet-related metabolic syndrome is the largest contributor to adverse health in the United States. However, the study of gene-environment interactions and their epigenomic and transcriptomic integration is complicated by the lack of environmental and genetic control in humans that is possible in mouse models. Here we exposed three mouse strains, C57BL/6J (BL6), A/J, and NOD/ShiLtJ (NOD), to a high-fat, high-carbohydrate diet, leading to varying degrees of metabolic syndrome. We then performed transcriptomic and genome-wide DNA methylation analyses for each strain and found overlapping but also highly divergent changes in gene expression and methylation upstream of the discordant metabolic phenotypes. Strain-specific pathway analysis of dietary effects revealed a dysregulation of cholesterol biosynthesis common to all three strains but distinct regulatory networks driving this dysregulation. This suggests a strategy for strain-specific targeted pharmacologic intervention of these upstream regulators informed by epigenetic and transcriptional regulation. As a pilot study, we administered the drug GW4064 to target one of these genotype-dependent networks, the farnesoid X receptor pathway, and found that GW4064 exerts strain-specific protection against dietary effects in BL6, as predicted by our transcriptomic analysis. Furthermore, GW4064 treatment induced inflammatory-related gene expression changes in NOD, indicating a strain-specific effect in its associated toxicities as well as its therapeutic efficacy. This pilot study demonstrates the potential efficacy of precision therapeutics for genotype-informed dietary metabolic intervention and a mouse platform for guiding this approach.


Introduction
The advancement of personalized medicine, an emerging medical paradigm in which therapeutic regimens are configured on an individual basis, will be critical for addressing public health issues, especially those related to environmental exposure.The importance of genotype and the epigenome in mediating phenotypic responses to environmental factors, and the overwhelming diversity of individual responses compared to population-level measurements, has become increasingly clear in recent years [1][2][3][4].However, the role of the epigenome, including DNA methylation and transcriptional regulation, has not been well characterized in this context.One such example is the known role of genotype in modulating how diet contributes to obesity and metabolic syndrome [5,6].Such "gene-by-diet", or GxD, interactions have been shown to explain, in great part, why a dietary recommendation that is beneficial for one demographic may be ineffective or even deleterious in another [7][8][9].Personalized nutritional guidelines are therefore one promising solution to addressing obesity, though constructing such guidelines will require a continued effort towards establishing an in-depth understanding of GxD effects as well as their epigenetic and transcriptomic bases.
Despite the clear importance of identifying GxD interactions governing health effects in humans, there are many limitations to studying GxD in human cohorts [10].For example, dietary backgrounds and habits of individuals are highly diverse and variable over a lifetime, translating to a lack of control over prior exposures and the need for extremely large cohort sizes in GxD studies to compensate for substantial confounding factors and noise.Conversely, controlled experiments on volunteers involving strict dietary regimens are likely to suffer from compliance issues as well as noise from variation in other environmental factors.Moreover, many GxD effects may be specifically mediated through the epigenome of metabolically-relevant tissues such as the liver and pancreas [11,12], but obtaining samples of these tissues from human patients is often difficult or infeasible.Because human studies face these numerous challenges, animal models represent a powerful alternative for studying GxD with the crucial advantage of allowing for extensive control and reproducibility of the experimental design.Furthermore, we have previously shown both that epigenetic analysis of mice reveals patterns of genetic susceptibility that are conserved in humans [13] and that genetically diverse mice have highly disparate phenotypic responses to diet [7].
The principles of GxD interactions, in which genotype is a key determinant of response to a particular diet or nutrient, also extend to other environmental exposures.Notably, they may even apply to the efficacy of therapeutic interventions such as the administration of drugs to alleviate metabolic disease.Currently, mouse studies that test the efficacy and safety of therapeutic regimens generally use only one laboratory strain such as C57BL/6J (BL6), which restricts their applicability to genetically diverse human patient populations.That is, a drug that happens to benefit the one specific mouse strain tested may nevertheless appear to fail efficacy tests in a larger clinical trial if other genetic backgrounds are insensitive to the drug; conversely, "rare" adverse effects that are actually genotype-dependent may be undetectable in one mouse strain [14].Therefore, the use of multiple genetically distinct mouse strains in drug trials has evident advantages in identifying genotype-dependent responsiveness and adverse effects and has been advocated for despite such experimental designs remaining rare [15].Consequently, GxD experiments have the potential to vastly improve the design of therapeutic trials, as the identification of disease-associated pathways altered in only a subset of genetic backgrounds can be used to predict that those genotypes, and those genotypes alone, will benefit from a drug targeting that pathway.
Here we have designed a new experimental paradigm in which we identify GxD changes in gene expression and DNA methylation in a cohort of genetically diverse mouse strains, applied the pathway analysis from those studies to inform a literature-based selection of a predicted strain-specific drug intervention, and assessed the molecular and phenotypic consequences of this intervention in a strain-specific manner (Fig 1A and 1B).Specifically, we identified a Farnesoid X receptor (FXR) agonist, GW4064 [16], as a candidate strain-specific therapy based on the observation of strain-specific modulation of the FXR pathway in response to a high-fat, high-carbohydrate "American" diet, and found that this drug not only has beneficial effects on the strain predicted to be responsive, but also induces strain-specific adverse responses in the strain predicted to be insensitive.In so doing, we have identified an example of both strainspecific disease phenotype mitigation and strain-specific deleterious side effects.This integrative approach thus opens the door to harnessing mouse genetics to improve preclinical assessment of both benefit and risk for genotype-and diet-specific drug candidates.

Characteristics of the mouse-diet experimental model
To elucidate how diverse genetic backgrounds differentially mediate phenotypic response to the environment, we designed epigenetic and transcriptomic analyses of three mouse genotypes exposed to two different diets.We selected three inbred founder strains from the Collaborative Cross project [17] and two diets with relevance to human nutrition: control laboratory (A) Genetically divergent mouse strains BL6, A/J, and NOD were exposed to either the American diet or a standard mouse chow.Phenotypic measurements were collected and WGBS and RNA-seq were performed on DNA and RNA extracted from liver.Arrows represent differences from the American diet vs. standard diet comparison, as reported in [7].(B) BL6 and NOD mouse strains were exposed to an American diet and were treated with the FXR agonist GW4064 or control vehicle.Phenotypic measurements were collected and WGBS and RNA-seq were performed on DNA and RNA extracted from liver.Arrows represent differences from the American diet + GW4064 vs. standard diet + vehicle comparison.Graphics were created with BioRender.com.https://doi.org/10.1371/journal.pgen.1010997.g001

PLOS GENETICS
chow (termed throughout as "standard") and a high-fat, high-carbohydrate "American" diet designed to match typical Western nutritional intakes which have been found to be detrimental to metabolic health in both animal and human studies [7,18].Mice were fed their assigned diets for 6 months.Macronutrient and lipid compositions of these diets are provided in S1 Table, and full ingredient, nutrient, and lipid compositions are detailed in [7].As we have previously reported, each of the three selected mouse strains has a different degree of phenotypic response to the American diet when compared to the standard diet [7].BL6, the most commonly used laboratory strain and the source of the standard mouse reference genome, displays strong negative changes in metabolic phenotypes on the American diet, including large increases in body fat, hepatic triglycerides, and total cholesterol [7].These phenotypic changes are in line with previous literature showing that the BL6 strain is adversely affected by high-fat diets [19,20].In contrast, the A/J and NOD/ShiLtJ (NOD) strains are more resistant to the American diet, with A/J mice displaying only mild changes in metabolic measurements and NOD showing moderate changes [7].This outcome is also in line with previous literature observations that the A/J strain is resistant to high-fat diets [21] and that NOD is a non-obese diabetic model [22].
The wide range of phenotypic responses in our chosen combination of strains and diets provides a robust experimental framework for identifying epigenetic GxD interactions (Fig 1A).This enables the identification of differentially expressed genes (DEGs) and differentially methylated regions (DMRs) denoting diet-dependent expression and methylation changes for each strain individually.Subsequently, by assessing the extent to which expression or methylation changes overlap between strains, "strain-specific" DEGs and DMRs can be identified that are altered on a genotype-by-diet basis.

Strain-specific gene expression analysis reveals common and distinct targets and pathways
First, we performed RNA sequencing (RNA-seq) on liver tissue from BL6, A/J, and NOD mice on the standard and American diets (n = 5 per group; n = 4 per group for NOD).The liver was chosen for its relative tissue homogeneity as well as its physiological relevance in metabolism and disease, including previously observed changes to hepatic triglyceride levels in all three of these strains in response to the American diet [7].Pseudo-alignment was performed for each strain relative to strain-specific reference transcriptomes in order to minimize strain and alignment biases in expression quantification, as described further in the methods section.Expression data was subsequently used to identify diet DEGs with a Benjamini-Hochberg (BH) false discovery rate (FDR)-adjusted p-value < 0.05 between the American vs. standard diet groups for each strain.
Our analysis identified 1,307 genes that were differentially expressed in at least one strain (Fig 2A and 2B and S2-S4 Tables).Overall, BL6 had by far the most diet-induced changes in expression (1,129 DEGs), while A/J displayed minimal expression changes (75 DEGs), and NOD had an intermediate number of changes (245 DEGs).This result closely mirrors the degree of phenotypic diet response in these strains [7].
Notably, nearly all diet DEGs are strain-specific, with only 26 common across all three strains (Fig 2B and S5 Table).Of these 26 strain-agnostic diet DEGs, a total of 21 are associated with metabolic pathways, with 18 involved in the regulation of or response to steroid and/or cholesterol levels as indicated by their annotations within the KEGG pathway database and the Ingenuity Knowledge Base (S5 Table) [23,24].These genes include the cholesterol transporter Abcg5, which helps to maintain cholesterol homeostasis [25]; Pcsk9, which functions to regulate plasma LDL cholesterol levels [26]; and several key enzymes and kinases in the cholesterol biosynthesis pathway (Acat2, Sqle, Sc5d, Hmgcs1, Mvk, Idi1, Cyp51, Dhcr7, Fdps, and Nsdhl) [24,27].
To determine processes or pathways where gene expression is uniquely modulated by diet in each strain, we performed Ingenuity Pathway Analysis (IPA) on the diet DEGs identified in each strain [24]  the mevalonate pathway I, and the superpathway of geranylgeranyl diphosphate biosynthesis I (via mevalonate) in particular.Overall, cholesterol biosynthesis pathways are highly enriched in all strains, in line with the significant changes in cholesterol levels observed to varying degrees in all three strains as well as the large proportion of strain-agnostic diet DEGs associated with cholesterol regulation.The remaining pathways are primarily unique to a subset of strains.For example, BL6 diet DEGs are uniquely enriched for atherosclerosis signaling, possibly reflecting the strain's highly adverse phenotypic responses to the American diet.Furthermore, while A/J and NOD had relatively few diet DEGs compared to BL6, the DEGs from these two strains are enriched for the estrogen biosynthesis pathway, which has been shown to play a role in the regulation of hepatic lipid metabolism [28] and is not significantly enriched in BL6 diet DEGs.Thus, these gene ontology results reveal potential epigenetic mechanisms by which strain-specific regulation of gene expression confers protective or deleterious responses to the same diet in different strains.

Applying GxD pathway analysis predicts strain-specific efficacy of metabolic drugs
The large GxD effects we observed in phenotypic and transcriptomic responses to the American diet indicate that genotype is a major factor in determining an individual's sensitivity to environmental challenges.Notably, these findings significantly reinforced our hypothesis that an individual strain's responsiveness to drug treatment, specifically one aimed at protecting against deleterious diet effects, could be equally genotype-dependent.Thus, we sought to determine if GxD pathways identified through our transcriptomic analysis as strain-specific could be utilized to predict certain drugs as having genotype-dependent efficacy.
To identify druggable gene targets which could elicit strain-specific protection against dietinduced obesity-associated phenotypes in our mouse strains, we performed a regulatory network analysis utilizing IPA, further described in our methods section, to identify upstream transcriptional regulators which could explain observed changes in gene expression [24].Of the master upstream regulators identified by this analysis, some are ubiquitously dysregulated by diet across all three strains, while others show a large degree of strain-specificity (Fig 3A).We then searched for commercially available metabolic drugs with established efficacy and characterized mechanism of function in BL6, namely a protective effect against a high-fat, high-carbohydrate diet, and asked whether any of them affected BL6-specific master upstream regulators.
Based on these criteria, we selected GW4064, a well-documented FXR agonist previously shown to prevent diet-induced obesity in BL6 mice [16].FXR (Nr1h4) is a master upstream regulator predicted via our IPA regulatory network analysis to be downregulated in BL6 and unchanged in A/J and NOD (Fig 3A and S4 Fig).We sought to further validate IPA's predicted BL6-specific dysregulation of FXR utilizing gene set enrichment analysis of each strain's diet DEGs using Enrichr [29].In concordance with IPA, this analysis reported significant enrichment (BH FDR-adjusted p-value < 0.05) of the NR1H4 gene set within the "ARCHS4 Transcription Factors Co-Expression" gene set library in the BL6 diet DEGs, but not in the A/J or NOD diet DEGs.As examples of individual genes affected by the activation of FXR, G6pc, an enzyme critical for the maintenance of glucose homeostasis, was downregulated upon American diet exposure in BL6, but not in A/J or NOD (Fig 3B).G6pc is predicted by IPA to be activated by FXR, and thus this expression pattern indicates a BL6-specific downregulation of the FXR pathway, in concordance with the IPA master regulatory analysis (Fig 3A).Additionally, the cholesterol transporter Abcg8 was significantly upregulated upon American diet exposure in A/J and NOD mice but was unchanged in BL6 (Fig 3C).Abcg8 is an indirect target of FXR and has been shown to upregulated by the FXR agonist GW4064 [30,31], suggesting that the BL6-specific downregulation of FXR removes this strain's ability to upregulate Abcg8 expression in response to the American diet, and further indicates that this response may be restored by GW4064 treatment.Overall, these strain-specific gene expression patterns observed within the FXR pathway suggest that it is specifically downregulated in response to the American diet in BL6 but not A/J or NOD mice.We thus reasoned that GW4064-induced activation of FXR might protect BL6 mice, but not A/J or NOD mice, from the deleterious effects of an American diet.

FXR agonist GW4064 exhibits strain-specific efficacy and toxicity
Due to A/J's relative phenotypic and transcriptomic insensitivity to the effects of the American diet in comparison to BL6 and NOD, we hypothesized that it would be difficult to distinguish the cause of A/J's responsiveness to GW4064 treatment.This is due to the potential ambiguity between a minimal A/J drug response stemming from either a lack of the initial dietary response or strain-specificity of the drug's efficacy.Thus, we tested the strain-specificity of the GW4064 treatment using BL6 and NOD only (Fig 1B).To do so, we exposed a cohort of 15-week-old BL6 and NOD mice (n = 5 per group) to the standard or American diet for 6 weeks.Mice on the American diet were treated with either GW4064 or vehicle via intraperitoneal injection; standard diet mice were treated with vehicle only.To characterize the health effects of diet and GW4064 treatment on these mice, we took a variety of metabolic measurements throughout the 6-week period, including lean and fat body mass, cholesterol levels, and hepatic triglyceride levels (S6 Table ).Four NOD mice were euthanized early due to health issues during testing and were removed from subsequent analyses.
Both the American diet and GW4064 treatment exerted strongly strain-specific effects on several metabolic phenotypes (Fig 4A).BL6 mice exhibited a large increase in body fat on the American diet that GW4064 treatment almost entirely reversed, consistent with previous reports and our earlier results (Fig 4A , left).In contrast, NOD mice did not exhibit any such increase in body fat and thus had no phenotype to revert with GW4064 treatment.Even more notable were trends in cholesterol and hepatic triglycerides which were both increased in both strains upon exposure to the American diet but reverted upon GW4064 treatment only in BL6 mice (Fig 4A, center and right).That is, treatment of American-diet mice with GW4064 caused statistically significant decreases in cholesterol and hepatic triglycerides in BL6 alone, with little to no change observed in NOD mice.
Hepatic histopathology also revealed strain-specific responses to GW4064 (S5 Fig and S7 Table ), though to a lesser degree.All mice from both strains on the American diet developed moderate to severe levels of hepatic lesions that were not seen in liver samples from mice of either strain when on the standard diet.This hepatic damage was prevented in all five BL6 mice on the American diet when also given GW4064.In contrast, one of three NOD mice on the American diet treated with GW4064 developed these hepatic lesions, suggesting a moderate strain-specific prevention of hepatic damage in response to GW4064 treatment.Observed hepatic lesions mainly consisted of increased glycogen and lipid deposition, which was low in all mice on the standard diet but increased significantly when mice were fed the American diet, irrespective of strain.The increase in hepatic glycogen deposition observed in the livers of both BL6 and NOD in response to the American diet is counterintuitive given the observed induction of insulin resistance in response to the American diet [7,32,33], which would typically increase glycogenolysis and result in decreased glycogen deposition [34,35].However, there have been several reported exceptions to this paradigm including cases of glycogen storage diseases contributing to hepatomegaly [36], as well as during hyperglycemic diet exposure [37].Further studies are needed to understand the mechanism of increased glycogen and lipid deposition in affected livers.Average levels of hepatic glycogen and lipid deposition were significantly lower in both strains fed the American diet while also receiving GW4064, albeit levels remained slightly higher in NOD compared to BL6, especially for lipid deposition which was consistent with hepatic triglyceride differences described earlier.Overall, this hepatic histopathology data further supports that the FXR agonist GW4064 has a strain-specific effect on mice fed an American diet, namely beneficial responses in BL6 mice compared to minimal responses in NOD mice.
To explore potential mechanisms underlying the strain-specific response to GW4064, we performed RNA-seq on livers of BL6 and NOD mice on the standard diet given the vehicle and on the American diet, under the vehicle or GW4064 treatment conditions.We first sought to validate the effect of GW4064 on the FXR pathway.To do so, we compared the expression of genes in IPA's Nr1h4 master regulatory network between the vehicle and GW4064 treatment conditions for each strain on the American diet.28 and 20 genes within this regulatory network were nominally differentially expressed (p-value < 0.05 and |log2FC| � 0.5) in BL6 and NOD, respectively (S8 and S9 Tables).12 of these genes are common to both BL6 and NOD, including Hmgcr, Fabp2, and Srebf2.The remaining strain-specific genes include Ppargc1a and Sc5d specific to BL6 as well as Fabp1 and Mvk specific to NOD.This differential expression observed in the Nr1h4 master regulatory network indicates successful modulation of FXR activity by GW4064 treatment.
We next performed a transcriptome-wide search to identify drug DEGs between the vehicle and GW4064 treatment conditions for each strain on the American diet (S10 and S11 Tables; BH FDR-adjusted p-value < 0.05).BL6 mice had a larger number of drug DEGs (237) compared to NOD (64), consistent with its stronger phenotypic response to GW4064 (Fig 4B).Furthermore, drug DEGs in BL6 and NOD are almost entirely strain-specific, with only 6 common to both strains.This suggests highly disparate responses to GW4064 between these genotypes and mirrors our earlier observations of strongly strain-specific responses to diet (Fig 2B).
To elucidate whether unique biological processes were affected by GW4064 treatment in BL6 compared to NOD, we used IPA to examine gene ontology of the drug DEGs.The top 10 most significantly enriched pathways (BH FDR-adjusted p-value < 0.1) for each strain are shown in Fig 4C and all significantly enriched pathways are shown in S6 and S7 Figs.BL6 drug DEGs are enriched in pathways highly relevant to metabolism, showing very strong enrichment in cholesterol biosynthesis in particular.By contrast, NOD drug DEGs show no enrichment for these same metabolic pathways, but instead have strong enrichments for immunerelated pathways including interferon signaling and activation of antiviral response pathways.This consistent enrichment in immune-related pathways among the NOD drug DEGs serves to potentially highlight that drug treatments, in addition to having strain-specific benefits resulting from the return of dysregulated gene networks to normal function, may also have strain-specific toxicities resulting from the augmentation of gene networks which were not dysregulated to begin with.
Due to the observed diet and drug pathway enrichment overlaps in BL6 and not NOD, we sought to further analyze the relationship between the effects of GW4064 treatment and the American diet on gene expression.To do so, we identified diet DEGs between the American vs. standard diets given the vehicle for each strain within the drug cohort (S12 and S13 Tables) and overlapped these DEGs with the previously defined drug DEGs (S14 and S15 Tables).126 of the 237 drug DEGs (53.2%) in BL6 are reversals of differential expression observed from the change in diet: i.e., 60 genes that were upregulated in BL6 mice on the American diet with the vehicle were downregulated upon GW4064 treatment, while 66 genes downregulated by the American diet with the vehicle were upregulated by GW4064.These include many functionally relevant genes such as the aforementioned Abcg5; Hmgcr, the rate-limiting enzyme in cholesterol biosynthesis [27]; and Scap, which regulates the SREBP family of transcription factors involved in the regulation of cholesterol and lipid metabolism [38].This suggests that, in BL6 mice, GW4064 is able to successfully reverse or protect from some of the dysregulated gene expression that arises upon exposure to the American diet.Conversely, only 4 of the 64 NOD drug DEGs (6.3%) exhibit this reversed differential expression.This result is consistent with the hypothesis that GW4064 does not significantly alter NOD's response to the American diet, whereas its effect on NOD mice leads to the dysregulation of pathways which were not affected by the change in diet, likely contributing to the observed toxicities.

Common and distinct targets and pathways in diet-associated DNA methylation changes
Given the capability of the epigenome, and in particular DNA methylation, to mediate phenotypic responses to both genetic and environmental factors, we were interested to see how methylation was altered in response to dietary changes in each of the three strains tested.To examine this relationship, we performed whole-genome bisulfite sequencing (WGBS) to assess DNA methylation in liver tissue of the same BL6, A/J, and NOD mouse individuals on which we performed RNA-seq (n = 5 per group; n = 3 per group for NOD).In order to accurately measure methylation in different mouse genotypes, which can be biased by alignment to incorrect reference genomes, we aligned WGBS reads to personalized reference genomes for each strain [39].We then performed a permutation-based analysis for each strain to identify diet-associated candidate DMRs and further filtered these DMRs by removing those with less than a 10% mean methylation difference between the standard and American diets.Due to the exploratory nature of this analysis, we chose to further analyze all nominally significant DMRs which also pass the post-hoc methylation difference filtering.1,316 DMRs were identified with a diet-associated methylation change in at least one strain (Fig 5A and S16-S18 Tables).BL6 contains the highest number of diet DMRs (728), which may reflect its high phenotypic responsiveness to the American diet relative to the other two strains, while NOD exhibits an intermediate number of diet DMRs (367), and A/J has the fewest (299).As observed with dietassociated expression changes, this trend perfectly matches the degrees of phenotypic changes associated with each strain.
Notably, of the 1,316 identified diet DMRs, only 17 regions (1.3%) exhibit diet-associated methylation changes in the same direction in all three strains.The vast remainder of the genome is instead characterized by strain-specific DMRs, with 1,258 (95.6%) unique to only one strain (Fig 5A ).One example of a strain-specific diet DMR is a region in the Srebf2 gene, which is involved in cholesterol homeostasis [40] and became hypomethylated only in BL6 mice, but not A/J or NOD mice, upon American diet exposure (Fig 5B ).As an example of a diet DMR that is common to all three strains, a region within the Adam11 gene became hypermethylated in all strains on the American diet (Fig 5C).Notably, Adam11 was also significantly overexpressed in all three strains on the American diet (S2-S4 Tables), potentially indicating a link between the methylation of this DMR and Adam11 gene expression.
Using the ChIPSeeker R package, which performs annotation of genomic regions [41], we next examined the gene associations of these diet DMRs.DMRs were associated with the promoter or gene body of 887 Ensembl genes; naturally, due to the strain-specific nature of most diet DMRs, each strain was associated with a distinct set of genes, with only 19 DMR-associated genes common to all three strains.We used IPA to identify enriched pathways in diet DMR-associated genes for each strain (Fig 5D and S8-S10 Figs).As we observed with dietassociated DEGs, cholesterol biosynthesis pathways and activation of the LXR/RXR pathway, which involves nuclear receptors with roles in lipid metabolism and transport [42,43], are enriched in all three strains, while other pathways are unique to certain strains; for example, the superpathway of geranylgeranyl diphosphate biosynthesis, which has been implicated in high-fat diet-associated non-alcoholic fatty liver disease fibrosis [44], is only enriched in A/J and NOD, while BL6 is uniquely enriched for the complement system, which has been previously associated with obesity and may also play a role in the development of insulin resistance and diabetes mellitus [45].Overall, our results indicate that, while certain dysregulated pathways appear to be universal or strain-agnostic effects of dietary intervention, there are many strain-specific methylation changes in distinct pathways which underlie the differential response of BL6, A/J, and NOD mouse strains to the same American-diet exposure.
We next compared the methylation and gene expression changes between the American and standard diets and observed that there is little overlap between diet DMRs and diet DEGs within each strain (S19-S21 Tables).In BL6, only 68 of 519 diet DMR-associated genes overlap with the 1129 diet DEGs, while in A/J, 17 of the 215 diet DMR-associated genes overlap with the 75 diet DEGs, and in NOD, 18 of the 250 diet DMR-associated genes overlap with the 245 diet DEGs.This occurs despite many enriched pathways, including cholesterol biosynthesis and FXR/RXR activation, being shared between the diet DEG and DMR gene sets within each strain (S1-S3 and S7-S9 Figs).This seemingly counterintuitive observation is partially driven by a subset of genes within the small overlapping set which are involved in these commonly enriched ontologies (10 of 68, 11 of 17, and 7 of 18 genes for BL6, A/J, and NOD, respectively; S19-S21 Tables).These include the previously discussed Abcg5 and Sc5d, which are each involved in the synthesis or regulation of cholesterol and are identified as both diet DEGs and diet DMR-associated genes in all three strains, as well as strain-specific overlaps such as the BL6-specific Hmgcr and the A/J-specific Mvd.The prevalence of these overlapping genes within these enriched ontologies common to the expression and methylation analyses suggests that the limited degree of concordance is functionally relevant.However, a portion of the genes within these commonly enriched pathways are indeed unique to either the diet DMRassociated gene set or the diet DEG gene set.For example, within the superpathway of cholesterol biosynthesis, Cyp51, Dhcr7, Hmgcs1, Hsd17b7, Mvd, Mvk, Idi1, Lss, Msmo1, Nsdhl, Pmvk, and Sqle were identified as BL6 diet DEGs with no observed differential methylation, while Dhcr24 was identified as a BL6 diet DMR-associated gene with no observed differential expression (S2 and S16 Tables).Overall, this suggests that diet exerts both overlapping and distinct effects on methylation and gene expression which contribute towards a common deleterious metabolic phenotype.

Strain-specific methylation responses to diet-protective drug treatment
Finally, to examine methylation patterns that might contribute to the observed strain-specific drug response, we performed WGBS on livers of the same BL6 and NOD mice from all three diet and drug treatment groups (standard diet with vehicle, American diet with vehicle, and American diet with GW4064).Using this exploratory small-sample WGBS study (n = 2 mice per group), we sought to confirm the presence of phenotypically relevant methylation patterns over obesity-and metabolism-related genes.We first called 355 and 309 "drug DMRs" between vehicle-and GW4064-treated mice on the American diet for BL6 and NOD, respectively, which represent regions where GW4064 treatment elicits a significant methylation change in each strain (S22 and S23 Tables).Of these DMRs, only 26 (7.32%) in BL6 and 15 (4.85%) in NOD also overlap a significant diet-induced DMR between standard and American diets given the vehicle (S24 and S25 Tables).This suggests that reversion of methylation changes caused by the American diet in BL6 accounts for only a subset of GW4064's effects, and that, instead, GW4064 promotes many methylation changes at regions where BL6 would not otherwise respond to American diet exposure.These GW4064-dependent changes, in turn, could represent protective metabolic effects that BL6 mice cannot normally produce on their own.We next examined NOD methylation over BL6 drug DMRs to determine how NOD mice might respond differently at these loci that undergo drug-induced changes in BL6.For 73 (20.6%) of the BL6 drug DMRs, BL6 but not NOD displayed diet-induced methylation change of at least 10% which is reverted with GW4064 treatment.For example, at the Crhr2 gene, which plays roles in lipid and cholesterol metabolism and has been implicated in obesity [46,47], BL6 had significant changes on the American diet which were restored to standard-diet levels with GW4064 treatment, while NOD had no change in any diet or treatment group (Fig  5E).Other patterns of strain-specific methylation differences were also observed.For 44 (12.4%) of the BL6 drug DMRs, BL6 mice were hyper/hypomethylated relative to NOD mice on the standard diet and this BL6 methylation level did not change upon American diet exposure, whereas GW4064 treatment induced a change in methylation to reach levels similar to those in NOD.This suggests that NOD's intrinsic methylation level over these regions could be protective against the American diet.An example of this pattern is observed in the intron of Ncoa2, a gene with important roles in adipogenesis and lipid metabolism [48][49][50] (Fig 5F).At another DMR, over the ncRNA 9530059O14Rik, NOD but not BL6 mice display hypermethylation upon American diet exposure (though the observed 8.8% methylation change is slightly under the 10% cutoff), while BL6 mice only gain methylation when treated with GW4064 (Fig 5G).This type of DMR could represent protective responses to the American diet in the NOD genotype that are not activated in BL6 mice, and that must instead be compensated for by GW4064 treatment.Furthermore, this DMR overlaps quantitative trait loci influencing obesity, hepatic cholesterol accumulation, and triglyceride concentrations [51,52].
Lastly, we searched for genomic regions in which methylation is only changed in GW4064-treated NOD mice.To do so, we identified NOD drug DMRs over which methylation changes were not observed between the American and standard diets given the vehicle in either BL6 or NOD mice, nor between the vehicle-and GW4064-treated BL6 mice on the American diet.These regions could represent GW4064-induced methylation changes that contribute to the NOD-specific immune activation predicted from the pathway analysis of the previously described transcriptomic data.288 (93.2%) of NOD drug DMRs fall within this category, and several can be functionally related back to the predicted drug-related immune response.As an example, a DMR following this pattern is observed over Map3k5 (also known as Ask1), a gene which plays a role in apoptotic signaling and has been shown to be induced by inflammatory cytokines [53] (Fig 5H).Overall, these DNA methylation and gene expression studies on GW4064-treated mice point toward a model in which strain-specific methylation and expression, including both intrinsic levels as well as disparate responses to an environmental exposure, can predetermine, in part, the consequences of a diet or the efficacy of a drug treatment for a given genotype.

Discussion
In this study, we coupled diet exposure with a set of genetically diverse mouse strains to ascertain the effect of genetic variation on epigenetic and transcriptomic responses to the environment, as measured by DNA methylation and gene expression.While a small number of metabolic genes were commonly activated or repressed across all mouse strains, such as in cholesterol biosynthesis pathways, a far larger set of genes falling into unique metabolic regulatory networks were altered on a strain-unique basis.We also applied our observations on gene-environment interactions to test strain-specific responses to the metabolic drug GW4064.We demonstrate that treatment with GW4064, while protective against the consequences of the high-fat, high-carbohydrate American diet in BL6 mice, had limited to no effect on dysregulated metabolic phenotypes in NOD mice and may even induce NOD-specific toxicities.Together, these results suggest that any given diet or drug does not induce consistent pathways of response across a population; rather, each individual's response is likely highly specific to that individual and is governed by the complex interactions between that individual's unique set of genetic variants.
FXR activation is not the only potential candidate strain-specific drug which has been identified through this analysis (Fig 3A).Inhibitors which target Slc13a1 and agonists which target Nr1i2 or Ago2 are predicted to have a similar BL6-specific protective effect against diet-induced obesity.Additionally, inhibitors of Insig1, Insig2, Acacb, or Mapk14 and agonists of Plau, Casp2, Egr1, or Eif6 are predicted to confer an A/J-specific protective effect.Finally, inhibitors of Rptor could have a NOD-specific protective effect.Perhaps most importantly, our integrative analysis also identifies pathways that may be druggable to confer a strain-agnostic protective effect across all three tested strains.These include inhibitors which target Dbp and agonists which target Scap, Insr, Atp7b, or Ncoa2.
It should be noted that "strain-specific" is a term necessary only for our pilot study, which is limited to three strains.With more genotypes-for example, the Collaborative Cross (CC) mouse panel in which alleles can be mapped with high resolution [54]-strain-specific DMRs and DEGs resolve into "non-conserved" or "genotype-specific" DMRs and DEGs present in a proportion of strains, which can be mapped to a genetic variant.In addition, these DMRs and DEGs can, in a sense, be treated as a more granular substitute for physiological phenotypes, which are few and broadly controlled.In this case, it is plausible that one large mapping panel could resolve a multitude of high-resolution genome-epigenome effects at once.This is an exciting prospect for any field, including dietetics, which would otherwise require the analysis of many complex interactions across heavily intertwined gene networks.Hypothetically, diets could even be reduced to their individual components, allowing for an even more detailed association of genotype with one nutrient of interest.Furthermore, a genetically diverse mouse panel such as the CC would allow for the use of model-based methods for the identification and analysis of GxD interactions [55,56].Such methods provide for a more statistically rigorous analysis but would require a higher degree of genetic diversity than the three strains analyzed in the current study.
In addition, this pilot study is limited in sample size which subsequently reduces statistical power, particularly for the identification of diet-and drug-associated DEGs and DMRs.Genome-wide corrections were applied when identifying DEGs, though only nominal p-values were utilized during DMR finding in order to broaden results.Future studies attempting to utilize a similar experimental and analytical methodology to identify transcriptomic and epigenetic drivers of strain-or genotype-specific phenotypic responses to environmental perturbations should seek to have additional biological replicates to overcome the stringency of these genome-wide corrections.
GW4064 has not been extensively explored in a clinical setting, largely due to its limited bioavailability and concerns that the presence of a stilbene group may cause hepatic toxicity as demonstrated previously in rats [57,58].In the present study, GW4064 treatment led to abnormal upregulation of immune and inflammatory response genes only in NOD mice, while BL6 mice demonstrated no adverse effects.These results have implications for how drug screening and preclinical animal trials could be performed.Specifically, they suggest that testing a diverse panel of mouse strains will prove far more valuable for identifying both genotype-specific responses and adverse effects compared to common experimental designs where compounds are tested on a single strain of laboratory mouse.Such an approach could identify drug candidates that, while appearing ineffective across the general population, have high efficacy for a subset of individuals, and can explain genetic risk factors underlying rare adverse events.We propose that utilizing such experimental designs will become the paradigm going forward, as we expect that doing so will greatly expand the pool of gene targets for clinical study.

Ethics statement
All animal protocols were approved by the University of North Carolina, North Carolina State University, and Texas A&M University Institution Animal Care and Use Committees.

Sample and animal information
For the diet study, 4-week-old A/J, C57BL/6J, and NOD/ShiLtJ mice were obtained from The Jackson Laboratory (Bar Harbor, ME) and acclimated for 2 weeks on a standard laboratory diet, then randomly assigned to standard or American diet groups, with five mice per strain and diet across two equally-sized cohorts studied in two locations: North Carolina (NC cohort) and Texas A&M University (TAMU cohort).Mice were fed their assigned diet in powdered form for roughly 6 months (24 weeks).Mice from the NC cohort were housed at the University of North Carolina during the first 4 months for analysis of body composition, metabolic rate, and physical activity.Mice were transferred to North Carolina State University for the final 2 months of diet exposure, where they also underwent necropsy, and tissue collection.Mice from the TAMU cohort were housed at TAMU for the duration of the analysis.Mice were housed five per cage and maintained at 22˚C under a 12-hr light cycle; they were maintained, and protocols were followed in accordance with the University of North Carolina, North Carolina State University, and Texas A&M University Institution Animal Care and Use Committee guidelines.Mice were sacrificed with carbon dioxide, and tissues were flash-frozen in liquid nitrogen or fixed in formalin.The complete protocols for animal handling and all phenotypic data are reported in [7].Only male mice from the NC cohort were used for RNA sequencing (RNA-seq) and whole-genome bisulfite sequencing (WGBS).
For the GW4064 drug study, 9-week-old male C57BL/6J and NOD/ShiLtJ mice were obtained from The Jackson Laboratory and acclimated for 6 weeks on a standard laboratory diet, then randomly assigned to standard-vehicle, American-vehicle, or American-GW4064 treatment groups, with five mice per strain and treatment.Four NOD mice, one on the standard diet, one on the American diet with vehicle, and 2 on the American diet with GW4064, were euthanized early due to health issues during testing.Mice were fed their assigned diet in pelleted form for roughly 6 weeks and treated with vehicle or GW4064 as described below (see "GW4064 formulation and treatment" section).Mice were housed at Texas A&M University at five per cage and maintained at 22˚C under a 12-hr light cycle; protocols were followed in accordance with Texas A&M University Institution Animal Care and Use Committee guidelines.Mice were killed with carbon dioxide, and tissues were flash frozen in liquid nitrogen or fixed in formalin.

Diet composition
Powdered and pelleted diets were designed in collaboration with Research Diets (New Brunswick, NJ); the American diet (D12052705) was based on the US Department of Agriculture's 2008 Dietary Assessment of Major Food Trends, as described in [7].A purified control mouse diet (D12052701, powder; D17031601, pellet; Research Diets) was used as a standard diet for comparison to the American diet.Diets were designed to recapitulate human diets as closely as possible, matching macronutrient ratio, fiber content, types of ingredients, and fatty acid ratios to the human diets.Accordingly, nutrient sources were selected to match intakes of human diets, e.g.beef protein to match red meat intake in the American diet.

Animal phenotyping
Total body weight was measured weekly for all individuals.Fat and lean mass were measured using echo magnetic resonance spectroscopy (EchoMRI, Houston, TX, USA).Fat mass and lean mass were recorded before and after the feeding trial.Body fat percentage is defined as the percentage of total body fat mass relative to body weight at the time of measurement.The percentage of body fat gained during the feeding trial reflects the calculated difference between the percentage of body fat at the end of the feeding trial and the percentage of body fat at the beginning of the feeding trial.Total cholesterol, HDL, and LDL measurements were performed as described in [7].Briefly, samples were measured in duplicate using the EnzyChrom AF HDL and LDL/VLDL Assay kit (BioAssay Systems, Hayward, CA, USA).Liver triglyceride levels were determined using the Folch extraction method, as previously described in [59].In brief, after performing necropsy the liver of the mouse was collected, frozen, and stored at -80˚C for subsequent analysis.50mg of liver tissue was homogenized in a 2:1 chloroformmethanol solution and allowed to equilibrate at room temperature for 15 minutes.After adding 100μL of 0.9% w/v NaCl, the samples were vortexed for 1 min and centrifuged at 2000 × g for 15 minutes at 4˚C.The lower organic phase was collected, evaporated using a nitrogen stream, and then resuspended in 500μL of a 0.5% Triton X-100/PBS solution.The samples were sonicated for 5 minutes using a Bioruptor and placed in a drying bath at 55˚C for 5 minutes.Infinity Triglyceride reagent (Thermo Scientific, USA) was added, the samples were incubated for 5 minutes at 37˚C, and absorbance at 500/660nm was measured to quantify triglyceride concentration according to the manufacturer's instructions.

GW4064 formulation and treatment
GW4064 was purchased from Cayman Chemical Company (Catalog Number 10006611) and stored at -20˚C in 20mg aliquots.Preparation of GW4064 solution for animal administration was performed on the day of injection.20mg of GW4064 was first dissolved in 1000μL of 99.5% DMSO, then diluted with 1985μL of water to reduce DMSO concentration.16μL of TWEEN 80 was then added to return GW4064 to solution.Vehicle solution was prepared with DMSO, water, and TWEEN 80 without GW4064.50mg/kg of GW4064 solution, and equivalent volumes of vehicle, were administered to mice via intraperitoneal injection twice a week.

DNA extraction and whole-genome bisulfite sequencing
Genomic DNA was extracted from liver samples using the Qiagen DNEasy kit, with an additional RNase incubation step (50μg/sample, 30 minutes) prior to column application to remove RNA.For the American vs. standard diet comparison, WGBS single indexed libraries were generated using NEBNext Ultra DNA Library Prep kit for Illumina (New England Bio-Labs) according to the manufacturer's instructions with modifications.500ng of input gDNA was quantified using the Invitrogen Qubit dsDNA BR assay and spiked with 1% unmethylated Lambda DNA (Promega, cat # D1521) to monitor bisulfite conversion efficiency.Input gDNA was fragmented by Covaris S220 or LE220 Focused-ultrasonicator to an average insert size of 350bp.Size selection to isolate insert sizes of 300-400bp was performed using AMPure XP beads.The EZ DNA Methylation-Gold Kit or EZ DNA Methylation-Lightning Kit (Zymo cat#D5005, cat#D5030) were used to bisulfite convert samples after size selection following the manufacturer's instructions.Amplification was performed after the bisulfite conversion using Kapa Hifi Uracil+ (Kapa Biosystems, cat# KK282) polymerase using the following cycling conditions: 98˚C 45s / 8cycles: 98˚C 15s, 65˚C 30s, 72˚C 30s / 72˚C 1 min.AMPure cleaned-up libraries were run on 2100 Bioanalyzer (Agilent) High-Sensitivity DNA assay, samples were also run on Bioanalyzer after shearing and size selection for quality control purposes.Quantification of libraries was performed by qPCR using the Library Quantification Kit for Illumina sequencing platforms (KAPA Biosystems, cat#KK4824) and the 7900HT Real-Time PCR System (Applied Biosystems).WGBS libraries were sequenced on an Illumina HiSeq2000 or HiSeq2500 instrument using 100bp paired-end indexed reads (v3 chemistry, BL6 and A/J samples) or 125bp paired-end indexed reads (v4 chemistry, NOD samples) with 10% PhiX spikein.For the GW4064 drug study, the above protocol was followed with the following modifications: libraries were dual indexed, size selection was performed using SPRIselect beads, qPCR quantification was performed using the CFX384 Real-time system (BioRad), WGBS libraries were sequenced on an Illumina NovaSeq6000 instrument using 150bp paired-end dual indexed reads (S4 flowcell, version 1.5 reagents) with 5% PhiX spike-in.

WGBS read alignment
TrimGalore (v0.6.6) was used to perform adapter removal and quality trimming of sequencing reads.In order to accurately estimate methylation while accounting for strain differences in genomic sequence, samples from BL6, A/J, and NOD were aligned to their respective reference genomes, obtained from UNC Systems Genetics (build 37), as we have previously described [39].Alignment was performed using Bismark (v0.23.0) and Bowtie2 (v2.9.2).Reference genomes were combined with the λ phage genome for measurement of conversion efficiency.The Bismark function deduplicate_bismark was then used to remove PCR duplicates.M-bias plots were generated using bismark_methylation_extractor with the-mbias_only flag in order to identify the positions of biased CpG sites most commonly resulting from library prep end-repair.bismark_methylation_extractor was then used to extract methylation values.For BL6 and A/J samples, 8 and 3 bp from the 5' and 3' ends of read 1, respectively, and 12 and 5 bp from the 5' and 3' ends of read 2 were removed using the-ignore,-ignore_3prime,-ignore_r2, and-ignore_3prime_r2 flags based on M-bias results.For NOD samples, 5 and 3 bp from the 5' and 3' ends of read 1 and 15 and 5 bp from the 5' and 3' ends of read 2 were removed.CpG positions from A/J and NOD were converted into the BL6 (mm9) reference coordinate system using modmap [60].These methylation values were used as input in the subsequent differential methylation analysis.

Differential methylation analysis
Raw CpG methylation data from the cytosine reports output by bismark_methylation_extractor were imported into R version 3.6.1 using the read.bismarkfunction of bsseq (v1.22.0)[61].DMR identification was performed using dmrseq (v1.6.0)[62] with the default DMR-finding parameters except for the following: minNumRegion = 3, maxPerms = 100.For diet DMRs, the raw BSmooth object was subset to CpGs where coverage was greater than 2x in 4 out of 5 samples from each strain/diet, except for NOD samples in which coverage had to be greater than 2x in 2 out of the 3 samples from each diet.dmrseq was run using diet as the test covariate separately for each strain.Significant DMRs were defined as those with a dmrseq nominal pvalue < 0.05 and a smoothed methylation difference between the groups, averaged across the entire DMR, of greater than 10%.To perform this calculation, CpG data were smoothed using the BSmooth/bsseq package (v1.22.0) with the default DMR-finding parameters (ns = 70, h = 1,000, maxGap = 10e8).Smoothing was performed over common and strain-unique CpGs to allow comparison of imputed methylation values across such sites, and then subset by coverage as previously described.To identify DMRs in the drug study, the same import, testing, and significance procedures and parameters as employed for the diet study were utilized, with a few differences.The raw BSmooth object was subset to CpGs where coverage was greater than 2x in all 12 samples.Two separate tests were run using dmrseq for each strain, one between the standard and American diet both given the vehicle and one between the American diet given the vehicle and the American diet given GW4064.Note that one mouse which was euthanized early was included in the drug DMR analysis (NOD sample on the American diet with GW4064 treatment).For the diet comparison, analysis was run between the American vs. standard diets within each strain.For the drug comparison, analysis was run between the American diet + GW4064 vs.American diet + vehicle within each strain, as well as between the American diet + vehicle vs. standard diet + vehicle within each strain.

DMR annotation
DMR annotation was performed with ChIPSeeker version 1.22.1 [41].Regions were annotated using the mm9 transcript database (TxDb.Mmusculus.UCSC.mm9.knownGeneR package version 3.2.2) and the genome-wide annotation for mouse (org.Mm.eg.dbR package version 3.8.2) with the promoter region defined as 3-kb upstream or downstream of the transcription start site.DMR-associated genes for each strain were defined as genes annotated for a nonintergenic DMR with at least 10% standard-American methylation difference in that strain.CpG islands were obtained via the UCSC Genome Browser.CpG shores were defined as the 2-kb regions upstream and downstream of a CpG island.

RNA extraction and RNA-seq
RNA was isolated from mouse liver using a Maxwell 16 LEV simplyRNA kit (Promega).For the American vs. standard diet comparison strand-specific mRNA libraries were generated using the TruSeq Stranded mRNA protocol (Illumina, cat# RS-122-2101).Libraries were performed following the manufacturer's protocol (Illumina, Part#15031050) with minor modifications.The input was 500ng (BL6 and NOD samples) or 2μg (A/J samples) and samples were fragmented for 6 min.The following PCR cycling conditions were used: 98˚C 30s / 13 (BL6 and NOD samples) or 12 (A/J samples) cycles: 98˚C 10s, 60˚C 30s, 72˚C 30s / 72˚C 5 min.Stranded mRNA libraries were sequenced on an Illumina HiSeq4000 (BL6 and NOD samples) or HiSeq2500 (v4 chemistry; A/J samples) instrument using 75bp (BL6 and NOD samples) or 70bp (A/J samples) paired-end indexed reads and 1% of PhiX control.For the GW4064 drug study, strand-specific mRNA libraries were generated using the NEBNext Ultra II Directional RNA library prep Kit for Illumina (New England BioLabs #E7760), and mRNA was isolated using Poly(A) mRNA magnetic isolation module (New England BioLabs #E7490).The preparation of libraries followed the manufacturer's protocol (Version 2.2 05/19).The input was 500ng and samples were fragmented for 15 min for an RNA insert size of ~200 bp.The following PCR cycling conditions were used: 98˚C 30s / 8 cycles: 98˚C 10s, 65˚C 75s / 65˚C 5 min.Stranded mRNA libraries were sequenced on an Illumina Hiseq4000 instrument using 48bp paired-end dual-indexed reads and 1% PhiX control.

RNA-seq read alignment, quantification, and analysis
RNA-seq reads were quantified using the kallisto program (v0.46.1) [63], which uses pseudoalignment to match reads with target genes.cDNA FASTA files for BL6, A/J, and NOD genomes were obtained from The Jackson Laboratory.Each file was used to generate a strainspecific reference index to which reads from the corresponding strain's samples were pseudoaligned and gene abundances estimated in order to minimize alignment bias from sequencing differences.RNA-seq data were analyzed in R version 3.6 using the edgeR (v3.28.1) and limma (v3.42.2) packages [64,65].Gene and transcript IDs were obtained from Ensembl release 101 and used as a target-mapping key to summarize kallisto abundance data at the gene level.Genes were filtered to those with a CPM (counts per million) greater than 1 in all 28 samples for the diet comparison and all 26 samples for the drug comparison; note that all mice which were euthanized early were removed prior to analysis.The normalization factors for library sizes were determined with edgeR using the TMM method.A contrast matrix was designed to look for differential expression in each strain individually.For the diet comparison, analysis was run between American vs. standard diets within each strain.For the drug comparison, analyses were run between American diet + GW4064 vs.American diet + vehicle within each strain, as well as between the American diet + vehicle vs. standard diet + vehicle within each strain.Raw counts were transformed to log-CPM values using the voom function from limma, then linear modeling was performed according to the contrast matrix to identify differentially expressed genes.Differentially expressed genes were defined as those with BH FDR-adjusted p-values less than 0.05.

Enrichment analyses
QIAGEN's Ingenuity Pathway Analysis software was utilized to perform pathway enrichment analyses of RNA-seq differential gene expression and WGBS differential methylation data [24].This enrichment analysis was performed on eight separate gene sets originating from the following three analyses: (1) differentially methylated genes associated with the change from standard mouse chow to the American diet from each of BL6, A/J, and NOD strains (three gene sets); (2) differentially expressed genes associated with the change from standard mouse chow to the American diet from each of BL6, A/J, and NOD strains (three gene sets); (3) differentially expressed genes associated with GW4064 treatment while on the American diet from the BL6 and NOD strains (two gene sets).For the gene sets of diet-induced differential methylation, expression analyses were performed using all genes associated with significant DMRs via ChIPSeeker.Default analysis settings were used except for the following: species was set to mouse only and tissues and cell lines were set to hepatocytes and liver.For the gene sets of diet-induced differential expression, an expression analysis was performed using the log2 fold change (log2FC) values and the BH FDR-adjusted p-values as inputs.Cutoffs of |log2FC| � 1 and BH FDR-adjusted p-value < 0.05 were used.Default analysis settings were used except for the following: species was set to mouse only and tissues and cell lines were set to hepatocytes and liver.For the gene sets of drug-induced differential expression, an expression analysis was performed using the log2FC values and the BH FDR-adjusted p-values as inputs.Cutoffs of | log2FC| � 0.5 and BH FDR-adjusted p-value < 0.05 were used.Default analysis settings were used except for the following: species was set to mouse only and tissues and cell lines were set to hepatocytes and liver.Gene set enrichment analysis was performed using Enrichr on differentially expressed genes (BH FDR-adjusted p-value < 0.05) associated with the change from standard mouse chow to the American diet from each of BL6, A/J, and NOD strains (three gene sets).

Ingenuity Pathway Analysis-Causal Network Analysis
QIAGEN's Ingenuity Pathway Analysis software was also utilized to perform Causal Network Analysis on differentially expressed genes associated with the change from standard mouse chow to the American diet from each of BL6, A/J, and NOD strains.Separately, an expression analysis with the settings and inputs described above was performed on the gene sets from each of these strains.A comparison analysis was then run between the three expression analyses.Information from the Causal Networks section of the Upstream Analysis tab is presented in this paper.A detailed description of the Causal Network Analysis tool is provided in [24].

Liver pathology
Liver samples were fresh fixed in 10% neutral buffered formalin before being processed in a Leica ASP300 tissue processor for paraffin embedding by the Texas A&M Rodent Preclinical Phenotyping Core.After embedding, 5μm sections were cut on a Leica 2165 rotary microtome and sections were H&E stained on a Leica HistoCore SPECTRA ST Stainer.After cover-slipping, slides were scored blinded by a board-certified veterinary pathologist.The severity of increased glycogen and lipid deposition was scored on a 0-4 scale with 0 = normal and 4 = severe.Note that none of the four mice which were euthanized early were included in this analysis.Hepatic intracellular glycogen was confirmed by staining with Periodic-acid Schiff

Fig 1 .
Fig 1.Schematic of experimental design.(A)Genetically divergent mouse strains BL6, A/J, and NOD were exposed to either the American diet or a standard mouse chow.Phenotypic measurements were collected and WGBS and RNA-seq were performed on DNA and RNA extracted from liver.Arrows represent differences from the American diet vs. standard diet comparison, as reported in[7].(B) BL6 and NOD mouse strains were exposed to an American diet and were treated with the FXR agonist GW4064 or control vehicle.Phenotypic measurements were collected and WGBS and RNA-seq were performed on DNA and RNA extracted from liver.Arrows represent differences from the American diet + GW4064 vs. standard diet + vehicle comparison.Graphics were created with BioRender.com.
. The top 10 most significantly enriched pathways (BH FDR-adjusted pvalue < 0.1) for each strain are shown in Fig 2C and all significantly enriched pathways are shown in S1-S3 Figs.Note that substantial overlap exists among several of these pathways, and thus similar gene sets drive the enrichments of each of the cholesterol biosynthesis pathways,

Fig 2 .
Fig 2. Gene-by-environment interactions of hepatic gene expression in BL6, A/J, and NOD mice.(A) Volcano plots of liver gene expression indicating significantly upregulated (red) and downregulated (blue) genes from the American vs. standard diet comparison of BL6, A/J, and NOD, as determined by the BH FDR-adjusted p-value.(B) Venn diagram of overlap between strains of differentially expressed genes (DEGs) upon American diet exposure.Marked categories are strain-specific, highlighting potential gene-by-diet interactions.(C) IPA results showing enriched pathways among diet DEGs (American vs. standard diet) for each strain; the top 10 most significantly enriched pathways with BH FDR-adjusted p-values < 0.1 are shown.https://doi.org/10.1371/journal.pgen.1010997.g002

Fig 3 .
Fig 3. Regulatory network analysis identifies ubiquitous and strain-specific transcriptional regulators of diet response.(A) Output of IPA regulatory network analysis identifying both common and strain-specific master upstream transcriptional regulators predicted based on diet DEGs (American vs. standard diet) in BL6, A/J, and NOD.The top 30 master upstream regulators ranked by z-score are shown.FXR (Nr1h4) is highlighted.(B) Plot of G6pc gene expression in BL6, A/J, and NOD mice on the standard and American diets.(C) Plot of Abcg8 gene expression in BL6, A/J, and NOD mice on the standard and American diets.The y-axes represent counts per million (CPM) of each gene.Error bars represent SE. * p < 0.05, ** p < 0.01, *** p < 0.001 from RNA-seq BH FDR-adjusted p-values.https://doi.org/10.1371/journal.pgen.1010997.g003

Fig 4 .
Fig 4. FXR agonist GW4064 elicits strain-specific responses in BL6 and NOD mice on the American diet.(A) Phenotypic measurements of BL6 and NOD mice under three conditions tested (standard diet + vehicle, American diet + vehicle, American diet + GW4064).Phenotypes shown are percent body fat gain (left), hepatic triglyceride levels (middle), and total cholesterol (right).Error bars represent SE. * p < 0.05, ** p < 0.01, *** p < 0.001 from ANOVA between means within strains, with post-hoc Tukey HSD test.(B) Volcano plots of liver gene expression indicating significantly upregulated (red) and downregulated (blue) genes from the American diet + GW4064 vs.American diet + vehicle comparison of BL6 and NOD, as determined by the BH FDR-adjusted p-value.(C) IPA results showing enriched pathways among drug DEGs (American diet + GW4064 vs.American diet + vehicle) for each strain; the top 10 most significantly enriched pathways with BH FDR-adjusted p-values < 0.1 are shown.https://doi.org/10.1371/journal.pgen.1010997.g004

Fig 5 .
Fig 5. Gene-by-environment interactions in hepatic DNA methylation of BL6, A/J, and NOD mice.(A) Heatmap of mean methylation differences between diets across a 3-kb window centered around the 1,316 diet DMRs, with each row corresponding to a single DMR.The x-axis represents the genomic distance of each locus relative to the center of the DMR and the color gradient represents the mean methylation difference between the American vs. standard diet at each locus.DMRs are further grouped by strain specificity, with group numbers shown to the right of the heatmap.DMRs in group 1 are nominally significant and have a dietary mean difference > 10% in all three strains; groups 2 through 7 denote DMRs which are nominally significant and have a mean difference > 10% in only one or two strains.(B) Example of a strain-specific diet DMR, showing the BL6-specific hypomethylation of Srebf2 and Mir33 upon exposure to the American diet.(C) Example of a strain-agnostic diet DMR, showing the ubiquitous hypermethylation of Adam11 upon exposure to the American diet.(D) IPA results showing enriched pathways among diet DMR- associated genes (American vs. standard diet) for each strain; the top 10 most enriched pathways with BH FDR-adjusted p-values < 0.1 are shown.(E-H) Examples of four distinct categories of diet-and drug-associated DMRs overlapping phenotypically relevant genes functionally implicated in various metabolic (E-G) or inflammatory (H) pathways.https://doi.org/10.1371/journal.pgen.1010997.g005