Introduction

A long-term lipid-rich diet is associated with multiple metabolic disorders, such as obesity (Hasegawa et al., 2020), cardiovascular disease (Lutsey, Steffen and Stevens, 2008; Maurya et al., 2023), and systemic low-grade inflammation (Duan et al., 2018; Christ, Lauterbach and Latz, 2019). The gastro-intestinal tract is the primary site of adaptation to dietary challenge, due to its roles in nutrient absorption, immunity and metabolism (Enriquez et al., 2022). Dietary challenges or other environmental or genetic factors can lead to prolonged inflammation and eventually damage the gastro-intestinal tract (Huang et al., 2017; Enriquez et al., 2022). IBD encompasses several chronic inflammatory gut disorders, including ulcerative colitis (UC) and Crohn’s disease (Chang, 2020; Adolph et al., 2022). Patients with IBD, have a higher risk of developing colorectal cancer (CRC), one of the most lethal cancers (Kim and Chang, 2014; Shah and Itzkowitz, 2022). The incidence of IBD has increased worldwide (Alatab et al., 2020; Freeman et al., 2021) during the last decade, in part due to increased consumption of lipid-rich diets (Maconi et al., 2010; Hou, Abraham and El-Serag, 2011). Furthermore, mouse studies show that HFD leads to more inflammation in the dextran sulfate sodium (DSS)-induced UC models compared to normal diet (Zhao et al., 2020). However, the response to HFD is variable across individuals (Zeevi et al., 2015) and the association between the lipid-rich diet and the risk of IBD in clinical studies is inconclusive (Kreuter et al., 2019), possibly due to genetic factors underlying inter-individual variability in gut inflammation and dysbiosis (Baumgart and Sandborn, 2012). More than 200 risk genes associated with IBD were identified through human genome-wide association studies (GWAS) (Huang et al., 2017), which have implicated epithelial function, microbe sensing and restriction, and adaptive immune response as drivers (Graham and Xavier, 2020; Kong et al., 2023). However, there is still no effective treatment for IBD. Current therapies, such as anti-tumor necrosis factor alpha (TNF-α) antibodies (Rutgeerts et al., 2005) and integrin α4β7 antibodies, blocking leukocyte migration (Feagan et al., 2013), can temporarily alleviate inflammation in a subset of patients (Rutgeerts et al., 2005; Feagan et al., 2013) but cause adverse effects (Harbord et al., 2017) and fail to prevent relapses (Doherty et al., 2018). Therefore, it is important to understand the gene-by-environment (GxE) interactions underpinning pre-clinical gut inflammation that eventually evolves into IBD, to aid in designing novel preventive and therapeutic strategies for intestinal inflammatory disorders.

Heterogeneity in clinical presentations as well as diversity in diet and lifestyle among human IBD patients render human genetic studies challenging (Molodecky et al., 2011). Experiments in laboratory mice allow to control several environmental factors, such as temperature and diet, when exploring the genetic modulators of IBD and also enable the collection of several relevant tissues to help elucidate tissue-specific mechanisms (Nadeau and Auwerx, 2019; Li and Auwerx, 2020). In addition, to mirror the heterogeneity of human populations, genetically diverse populations, such as mouse genetic reference populations (GRPs), can be used in a systems genetics paradigm (Nadeau and Auwerx, 2019; Li and Auwerx, 2020). This not only allows the mapping of clinically relevant traits in controlled environments but also the characterization of intermediate molecular phenotypes from tissues that cannot easily be obtained in humans (Williams et al., 2016; Li and Auwerx, 2020). For example, studies on the molecular basis of non-alcoholic fatty liver disease in the Collaborative Cross founder strains illustrated the importance of the genetic background in determining susceptibility to steatosis, hepatic inflammation and fibrosis (Benegiamo et al., 2023). Moreover, the BXD GRP was used to identify genetic variants associated with metabolic phenotype variation, such as bile acid homeostasis (Li et al., 2022), lipid metabolism in plasma (Jha, McDevitt, Halilbasic, et al., 2018) and liver (Jha, McDevitt, Gupta, et al., 2018) as well as mitochondrial dysregulation (Williams et al., 2016), using GWAS or quantitative trait locus (QTL) mapping (Wu et al., 2014; Williams et al., 2016). Thus, large mouse GRPs are useful tools for identifying the tissue-specific mechanisms of complex diseases.

In order to decipher the genetic and environmental contributions to the development of intestinal inflammation, we measured the colon transcriptome of 52 BXD strains fed with CD or HFD (Williams et al., 2016). HFD feeding from 8 to 29 weeks of age induced an IBD-like transcriptomic signature in colons of some, but not all, BXD strains, uncovering a subset of BXD strains that could be susceptible to HFD-induced IBD-like state. Gene co-expression analysis revealed two IBD-related modules in the colons of HFD-fed mice, one of which is likely under the control of a ModQTL. Through a systems genetics prioritization of genes under this ModQTL, we identified candidate IBD-related genes that we validated using GWAS in the UK Biobank (UKBB) for human IBD.

Results

HFD feeding leads to highly variable transcriptomic adaptations in the colon of BXD strains

For this study, we used an extensively characterized BXD mouse panel of 52 BXD strains fed with a chow diet (CD) or high-fat diet (HFD) from 8 to 29 weeks of age (Williams et al., 2016; Jha, McDevitt, Gupta, et al., 2018; Jha, McDevitt, Halilbasic, et al., 2018), in which we mapped genetic determinants of metabolic traits in the liver (Williams et al., 2016; Jha, McDevitt, Gupta, et al., 2018) and plasma (Jha, McDevitt, Halilbasic, et al., 2018). These mice underwent metabolic phenotyping, with many metabolic traits being altered by HFD (Figure 1A), and multiple organs were harvested and flash-frozen for future use (Williams et al., 2016). Here, we focused on proximal colon samples from this population and performed microarray-based transcriptome analysis of this tissue (Figure 1A).

Figure 1.

Principal component analysis (PCA) of all transcriptomes (Figure 1—figure supplement 1A) showed that the first principal component (PC1) separated mice by diet, indicative of a global diet effect in the population. Nevertheless, transcriptomes of several strains (such as BXD12, BXD84 and BXD81) on HFD had very similar PC1 values to their CD counterparts (Figure 1—figure supplement 1B), suggesting that they were resistant to dietary changes. Similarly, BXD strains did not cluster completely by diet based on hierarchical clustering analysis indicating that the genetic differences can override the impact of diet on the transcriptome in the colon (Figure 1—figure supplement 1C). To obtain a global, strain-independent, view of the HFD effect, we performed a differential expression analysis and identified 115 up-and 295 down-regulated differentially expressed genes (DEGs, absolute Log₂(Fold Change) > 0.5 and Benjamini-Hochberg (BH)-adjusted P value < 0.05, Figure 1B). Of note, Cldn4, one of claudins implicated in intestinal permeability (Ahmad et al., 2017), was significantly down-regulated and serum amyloid A (Saa1 and Saa3), which have been involved in the inflammatory response (Ye and Sun, 2015; Tannock et al., 2018), were up-regulated upon HFD (Figure 1B). Furthermore, gene set enrichment analysis (GSEA) showed an upregulation of inflammation, cell proliferation and translation, mitochondrial respiration, and stress response-related pathways upon HFD, while genes involved in the intermediate filament - that contribute to maintaining intestinal barriers (Misiorek et al., 2016; Mun, Hur and Ku, 2022) - were down-regulated (Figure 1C). All in all, the transcriptome data are consistent with an HFD-induced downregulation of components of the intestinal barrier, enhanced permeability, induction of the unfolded protein response (UPR) and increased inflammation in BXD colons, much like HFD does in humans (Bischoff et al., 2014). However, as in humans, not every strain exhibited the same response to dietary challenges. GSEA analyses applied individually to the diet effect in each strain showed a high degree of diversity in the inflammatory response (Figure 1— figure supplement 1D). For example, BXD44, 45 and 55, highlighted in red, were the 3 most susceptible strains to gut inflammation upon HFD, whereas BXD1, 67, and 85, colored in green, showed no significant enrichment in gut inflammation. This diversity in responses provided the basis for a systems genetics investigation of HFD-driven gut inflammation determinants in the BXD.

The transcriptomic response to HFD of a subset of BXD strains resembles DSS-induced ulcerative colitis (UC)

IBD is characterized by increasing inflammation in the gastro-intestinal tract (Adolph et al., 2022). To investigate the disease relevance of the chronic inflammation seen in BXD colons upon HFD, we extracted the transcriptomic signatures from DSS-induced mouse UC models (Czarnewski et al., 2019) and two IBD human studies (GSE16879 (Arijs et al., 2009) and GSE83687 (Peters et al., 2017), Materials and Methods) and used these signatures as custom gene sets in GSEA on the global HFD effect. DSS is widely used to induce UC in mouse models and disease severity increases over time (Czarnewski et al., 2019). GSEA analyses showed that DSS-induced genes from days 4 (early inflammatory phase), 6 and 7 (acute inflammatory phase) were significantly enriched in genes upregulated by HFD, especially the dysregulated genes in the later stage of DSS-induced UC (Figure 1D, bottom panel). Similarly, genes involved in human IBD (UC and Crohn’s disease (CDs)) were also enriched in those same genes (Figure 1D, top panel). The same trend was observed for downregulated genes in mouse and human IBD, which were negatively enriched (Figure 1D), illustrating that HFD induced an IBD-like transcriptomic signature in BXD colons.

While the average response across all BXDs shared features of mouse and human IBD, we assessed the strain-specificity of this response by measuring each strain’s response to IBD using GSEA (Figure 2A). Hierarchical clustering of the normalized enrichment scores (NES) in mouse IBD datasets classified the BXDs into three groups: susceptible strains highlighted in red (19 strains), intermediate strains represented in blue (11 strains), and resistant strains colored in green (17 strains) (Figure 2A, top panel). Of note, in line with colon histological lesions comparison of DSS-induced colitis mouse models in the literature (Mähler et al., 1998), the C57BL/6J strain, one of the parental strains of the BXDs, was classified as one of the susceptible strains while the other parental strain DBA/2J belonged to the resistant group (Figure 2A, top panel), suggesting that genetic determinants inherited from the parental strains may determine the susceptibility of BXD strains to HFD-induced IBD-like inflammation in the colons.

Figure 2.

To establish the functional relevance of this transcriptome-based classification on systemic inflammation, we compared plasma cytokine levels of these three groups under HFD (Williams et al., 2016). Interestingly, the susceptible group have significantly lower levels of the anti-inflammatory cytokine - Interleukin (IL)-10 (Figure 2B, two-tailed t-test p < 0.01) and increased the proinflammatory cytokine - IL-15 (Figure 2C, two-tailed t-test p < 0.0001) compared to the resistant strains. IL10 itself has been identified as an IBD-related candidate gene using GWAS in humans (Franke et al., 2008) and IL-10-deficient mice are also well-known mouse model for IBD research (Keubler et al., 2015). IL-15 is another important cytokine involved in intestinal inflammation and is elevated in the human guts with IBD (Liu et al., 2000). IL-15 knock-out mice are also reported to have less severe symptoms, such as weight loss and histological scores, following DSS administration (Yoshihara et al., 2006). In summary, susceptibility to HFD-induced IBD-like inflammation in the colon, as assessed by changes in levels of genes associated with IBD, correlates with markers of the general inflammatory status of mice.

Identifying IBD-related gene modules in BXD colons

Since different BXD strains seem to exhibit different susceptibility to IBD, we set out to explore gene expression signatures underlying these differences. For that, we used Weighted Gene Co-expression Analysis (WGCNA) to construct CD-and HFD-specific gene co-expression networks to identify modules of co-expressed genes (Figure 3A, Appendix 1 - Table 1). Disease-associated modules were then defined as modules under HFD are significantly enriched in mouse DSS-induced UC signatures by an over-representation analysis (ORA, BH-adjusted P value < 0.05 and number of enriched genes > 5, Figure 3A). The HFD co-expression network consisted of 39 modules ranging in size from 34 to 1,853 genes and containing a total of 14,723 genes (Appendix 1 - Table 1). We visualized this network using Uniform Manifold Approximation and Projection (UMAP) (Figure 3B), reflecting that the majority modules were closely connected in the co-expression network.

Figure 3.

Enrichment analyses indicated that modules HFD_M9 (484 genes), HFD_M16 (328 genes), and HFD_M28 (123 genes) were enriched with genes that are upregulated by DSS-induced colitis, while HFD_M15 (368 genes), HFD_M24 (159 genes), and HFD_M26 (135 genes) were significantly enriched with downregulated genes (Figure 3C). Of note, more than 20% of genes involved in HFD_M9 and HFD_M28 were part of the dysregulated genes of the acute phase of mouse UC (day6 and day7) (Figure 3C). Interestingly, genes perturbed during IBD pathogenesis in humans were also enriched in HFD_M9 and HFD_M28 (Figure 3C).

While IBD-related genes were predominantly found in HFD modules, we also found that two modules, CD_M28 (185 genes) and CD_M32 (142 genes), in CD-fed mouse colons were associated with IBD (Figure 3—figure supplement 1A). These two-modules significantly overlapped with the IBD-related HFD_M9 and HFD_M28 modules, respectively (BH-adjusted P value < 0.05) (Figure 3—figure supplement 1B). Moreover, the molecular signatures underlying human UC and Crohn’s disease were also clustered in these two modules (CD_M28 and CD_M32) under CD (Figure 3—figure supplement 1C). Collectively, the co-expression and enrichment analyses identify HFD_M9 and HFD_M28 as IBD-related modules on which we focus our subsequent investigation.

Identifying biological roles and transcriptional regulation of the IBD-related modules

To identify the biological function of the IBD-related modules, we performed enrichment analyses using the Hallmark database and the cell-type gene signatures (Kong et al., 2023) (Materials and Methods). Genes in HFD_M9 were enriched in KRAS signaling and inflammation-related pathways, while HFD_M28 was enriched in IFN-α/γ responses (BH-adjusted P value < 0.05) (Figure 4A). Both modules were enriched in IFN-γ response genes (Figure 4A). IFN-γ is an essential cytokine for innate and adaptive intestinal immune responses (Brasseit et al., 2018). It has been reported to play a key role in mouse (Ito et al., 2006) and human (Tilg et al., 2002) IBD pathogenesis, and was identified as a potential therapeutic target to alleviate inflammatory response in IBD (Li et al., 2021). In addition, genes that are dysregulated in immune cells of Crohn’s disease patients (Macrophages, B cell and immune cycling cells) were enriched in HFD_M9 (Figure 4B). In contrast, genes of HFD_M28 were not only enriched for genes that are dysregulated in immune cells, but also in intestinal epithelial cells of diseased individuals, such as Goblet and stem cells (Figure 4B). Overall, HFD_M9 and HFD_M28 are both involved in inflammatory response, while genes involved in HFD_M28 also potentially influence intestinal epithelial barrier.

Figure 4.

To identify transcriptional drivers of the two IBD-related modules, we performed a transcription factor (TF) enrichment analysis (Materials and Methods) and found that ZIC2, SMAD3, REL, FOSL1, and BATF are the top enriched transcription factors for the genes in HFD_M9 (Figure 4C), while the expression of genes in module HFD_M28 may be regulated by Interferon regulatory factors (IRFs, IRF1, IRF2, IRF7, and IRF9) and the signal transducer and activator of transcription families (STAT, STAT2) (Figure 4D). In fact, most of these TFs have been reported to be involved in gut inflammation. For example, Smad3 mutant mice were more susceptible to intestinal inflammation (Yang et al., 1999). Moreover, the IFN-STAT axis is essential to initiate the type-I IFN induction that is critical for human immune defense, such as IBD diseases (Stolzer et al., 2021) and primary immunodeficiency diseases (Mogensen, 2019) as well as for disease tolerance (Mottis et al., 2022). Collectively, we have identified TFs that likely control the expression of the two IBD-related modules to play an essential role in gut inflammation regulation.

Identifying ModQTLs for IBD-related modules and filtering of candidate genes

To analyze how the genotype impacts the IBD-like inflammatory response associated to HFD, we performed module QTL mapping analysis (ModQTL) for both IBD-related modules (HFD_M9 and HFD_28) (Figure 5A). We found a suggestive QTL for HFD_M28 (P value < 0.1), on chromosome 16 containing 552 protein-coding genes (Figure 5A, Appendix 1 - Table 2). The ModQTL analysis was also performed on the modules that are significantly enriched in IBD-downregulated genes (HFD_M15, HFD_M24, and HFD_M26), but no significant or suggestive QTLs were detected. Therefore, we focused on the QTL for IBD-induced genes in HFD_M28 and annotated its candidate genes based on three criteria (Figure 5B): (1) presence of high-impact genetic variants (such as missense and frameshift variants) in BXDs, (2) association with inflammation based on literature mining (Materials and methods), (3) presence of cis-expression QTLs (eQTLs), that is, whether the expression of the gene is controlled by the QTL. The 27 genes satisfying at least two of the above criteria were considered as candidate genes driving the expression of module HFD_M28 (Figure 5C).

Figure 5.

To further prioritize candidate genes regulating module HFD_M28, we applied GWAS to detect Crohn’s disease-and UC-associated genetic variants using whole genome sequence (WGS) dataset in UKBB (Figure 5C, Materials and Methods). Interestingly, the genetic variants of two genes under the QTL peak, i.e, ephrin type A receptor 6 (EPHA6, P value = 2.3E-06) (Figure 5C, Figure 5—figure supplement 1A, Appendix 1 - Table 3) and Mucin 4 (MUC4, P value = 1.2E-06) (Figure 5C, Figure 5—figure supplement 1B, Appendix 1 - Table 4) were also associated with UC in humans. EPHA6 belongs to Eph/Ephrin Signaling and this pathway has been associated with gut inflammation (Coulthard et al., 2012) and proposed as a potential target to alleviate the inflammatory response in IBD (Grandi et al., 2019), but the association between EPHA6 and IBD is not explored yet. The Gene-Module Association Determination (G-MAD) (Li et al., 2019) (https://systems-genetics.org/gmad) also revealed that expression of Epha6 in mouse gastro-intestinal tract correlates with genes involved in inflammation-related pathways, such as IL-6 production and regulation of inflammatory response (Figure 5D, Appendix 1 - Table 5). MUC4 is a transmembrane mucin (Gao et al., 2021) and highly expressed in gastro-intestinal tract according to the human protein atlas (Uhlén et al., 2015) (https://www.proteinatlas.org/humanproteome/tissue/intestine) (Figure 5C). The expression of MUC4 in the human gastro-intestinal tract correlates with genes that are enriched for CRC and O-linked glycosylation based on G-MAD (Li et al., 2019) (Figure 5E, Appendix 1 - Table 5). O-linked glycans are expressed by the intestinal epithelium to maintain barrier function, especially mucin type O-glycans, and gut disorders can be affected by dysfunction of O-linked glycosylation (Brazil and Parkos, 2022). Moreover, MUC4 is upregulated in enterocytes and Goblet cells in colons of Crohn’s disease patients (Figure 5F). MUC4 hence is a strong candidate because of its role in maintaining the intestinal epithelium and controlling the gut inflammatory response (McGuckin et al., 2011) and EPHA6 might be a novel candidate gene to impact gut inflammation. Based on the results of our QTL mapping, human GWAS in UKBB, and existing literature, we hypothesize that MUC4 and EPHA6 impact on colon integrity and inflammation and may be important players in gut inflammation or IBD triggered by an unhealthy, lipid-rich diet.

However, it is unclear through what mechanisms the genetic variants in the candidate genes affect IBD susceptibility. One possibility is that genetic variation leads to altered levels of expression of the gene, ultimately affecting disease susceptibility. To test this possibility, we examined the GTEx resource (GTEx Consortium, 2013) and found that MUC4, but not EPHA6, has cis-eQTLs in the sigmoid and transverse colon. To establish likely causal links with IBD incidence, we used these associations as instruments in a two-sample Mendelian randomization (MR) (Hemani, Tilling and Smith, 2017; Hemani et al., 2018) analysis. Using publicly available GWAS summary statistics for IBD, Crohn’s disease, and ulcerative colitis (Liu et al., 2015; Elsworth et al., 2020) as outcomes, we found suggestive evidence that increased expression of MUC4 in the sigmoid, but not transverse, colon may increase the risk of IBD (nominal P value = 0.033, Appendix 1 - Table 6). No eQTLs were reported for EPHA6 in the colon, precluding us from investigating the potential consequences of changes in its expression in these tissues.

Discussion

Dietary, environmental and genetic factors have all been reported to influence intestinal inflammation (Adolph et al., 2022). Indeed, HFD can impair the intestinal epithelial barrier and trigger pre-clinical inflammation in the gastro-intestinal tract, eventually leading to inflammatory disorders of the gut (Enriquez et al., 2022). In addition, genetic factors identified by GWAS can also predispose to IBD. For example, the interleukin-1 and −7 receptors (IL-1R2 and IL-7R) were identified as candidate genes that regulate the immune response in IBD (Khor, Gardet and Xavier, 2011). However, the heterogeneity of diet and other environmental factors in human studies limits our ability to identify GxE interactions and pinpoint the genes and pathways involved in diet-induced gut inflammation. Studies in model organisms such as the mouse, where the environment can be carefully controlled, provide a valuable complement to human genetics studies that by nature are mainly observational (Nadeau and Auwerx, 2019; Li and Auwerx, 2020). Unfortunately, most mouse studies only evaluate mice from a single genetic background, limiting their generalizability and translatability to humans (Nadeau and Auwerx, 2019; Li and Auwerx, 2020). Conversely, GRPs such as the BXDs can mimic at least in part the heterogeneity of human populations and allow us to estimate the effect of GxE interactions on complex diseases (Jha, McDevitt, Gupta, et al., 2018; Li and Auwerx, 2020).

Here, we utilized a panel of 52 BXD genetically diverse mouse strains fed with either HFD or CD to explore the genetic and dietary modulators of inflammation seen in the colon transcriptomes using systems genetics approaches. The colon transcriptomic response to HFD in this mouse population recapitulated several of the general features observed in DSS-induced UC mouse models and human IBD patients. In particular, we identified the upregulation of inflammation-related genes and the UPR as well as the downregulation of intercellular adhesion-related genesets as common signatures induced by HFD (Kreuter et al., 2019). Moreover, our dataset not only was informative about the transcript changes of IBD at the population level, but also unveiled extensive strain-specific effects that allowed us to classify strains based on their propensity to develop IBD-like signatures. The fact that these susceptibility groups also differed in anti-and pro-inflammatory plasma cytokine levels (IL- 10 and IL-15, respectively) suggests a relation between these tissue-specific transcriptional signatures and systemic low-grade inflammation. Since gene interactions determine cellular processes and the molecular functions of correlated genes are often similar (Nayak et al., 2009), we attempted to elucidate the mechanism underlying the diversity of IBD-like signatures and chronic inflammation in BXD colons using gene co-expression analyses. This led us to identify two IBD-related gene modules (HFD_M9 and HFD_M28).

As most differentially expressed genes are likely to be driven by and not be a cause of disease (Porcu et al., 2021), we attempted to understand whether the signatures in the colon are causes or consequences of chronic inflammation. A first step was to characterize possible transcriptional and genetic regulators of IBD-related modules. Enrichment analyses showed that both IBD-associated modules largely consisted of immune response-related genes. Specifically, genes involved in HFD_M9 and HFD_M28 are both differentially expressed in immune cells in inflamed tissues of Crohn’s disease patients (Kong et al., 2023). Moreover, the HFD_M28 module was enriched for TF motifs of STAT2 and IRF family, and HFD_M9 for SMAD3 and REL, which were illustrated to control the expression of these gut inflammation-related genes, and influence the inflammatory response triggered by HFD in the colon.

While we found IBD-related gene modules and the TFs driving their expression, the genetic drivers of the diversity of gut inflammatory responses observed across the BXDs remained elusive. To find candidate genes causing gut inflammation upon HFD, we then performed Module QTL (ModQTL) analysis and allocated a suggestive ModQTL that may be controlling one of IBD-related module (HFD_M28) under HFD. Importantly, through our prioritization scheme for the genes under the ModQTL, we identify two plausible candidates, Epha6 and Muc4, that have high-impact variants in the BXDs, are related to inflammation, and harbor variants in humans that are associated with IBD based on UKBB GWAS result. Mendelian randomization analysis suggests that higher expression of MUC4 in the sigmoid colon may increase the risk of IBD. Furthermore, Muc4 knock-out mice have been shown to be more resistant to DSS-induced UC through upregulating the expression of Muc2 (mucin secretion) and Muc3 (transmembrane mucin) (Das et al., 2016). A GWAS study also indicated that mutations in EPHA6 increase risk for CRC (Guda et al., 2015), but its potential association with IBD is a new finding. Therefore, these results point to important potential roles of Muc4 and Epha6 in gut chronic inflammation leading to inflammatory gut disorders.

Although studies in the BXD cohort are limited to variants present in the parental strains, C57BL/6J and DBA/2J, our analysis nevertheless shows how genetic diversity in this population allows us to detect the genetic modulators of chronic intestinal inflammation, that are more difficult to identify in widely used IBD mouse models on a single genetic background. In support of the generalizability of our data, the identified candidate genes in our mouse models were also associated to human UC, demonstrating that chronic inflammation induced upon HFD feeding may indeed be a prelude to human UC.

In conclusion, our systems genetics investigation of the colon in a controlled GRP, complemented with human GWAS studies, enabled the prioritization of modulators of IBD susceptibility that were generalizable to the human situation and may have clinical value.

Materials and methods

Population handling

Mice were studied as previously described (Williams et al., 2016) and multiple organs were harvested for further analysis. Briefly, in groups of 3-5 animals from the same strain and diet, in isolator cages with individual air filtration (500 cm2, GM500, Tecniplast) and provided water ad libitum. Mice were fed CD ad libitum until 8 weeks of age. From 8 weeks to 29 weeks, half of the cohort was fed ad libitum HFD and the rest continued to be fed a CD (Figure 1A). CD composition: 18% kCal fat, 24% kCal protein and 58% kCal of carbohydrates (Teklad Global 18% Protein Rodent Diet 2018 chow diet, Envigo, Indianapolis, USA). HFD composition: 60.3% kCal fat, 18.4% kCal protein and 27.3% kCal of carbohydrates (Teklad Custom Diet TD.06414, Envigo, Indianapolis, USA). All mice were fasted overnight (from 6pm to 9am) prior to euthanasia. All procedures were approved by the veterinary office of canton Vaud under animal experimentation license number VD2257. In this work, proximal colons were extracted from the bio-banked samples and we did not use any new animals.

Transcriptome of the proximal colon in BXDs

A ∼1 cm portion of the proximal half part of the colon was excised following euthanasia, washed in PBS and immediately stored in liquid nitrogen. Approximately 5 animals of the same strain fed the same diet were pooled at equal mass concentration for further RNA extraction. Total RNA was extracted using Direct-zol (Zymo Research) including the DNase digestion step. 100ng of total RNA was amplified using the Ambion® WT Expression Kit from Life Technologies (part number 4411974) and 5,500ng of cDNA was fragmented and labeled using the Affymetrix WT terminal labeling kit (part number 900671) all following manufacturers protocols. Labeled cDNA was hybridized on an Affymetrix Clariom S Assay microarray platform (GPL23038) in ∼16 hours of incubation, then washed and stained using an Affymetrix 450 Fluidics Station according to Affymetrix protocols. Finally, arrays were scanned on Affymetrix GSC3000 7G Scanner. Microarray data preprocessing was performed using apt-probeset-summarize from the Array Power Tool (APT) suite (v2.11.3) with the gc-sst-rma-sketch standard method and resulting expression values were log-transformed. Microarray probes targeting polymorphic regions in the BXD population were ignored in the process. For probesets targeting a same transcript, only the probeset with the highest value was considered.

Differential gene expression analysis

General differences in mRNA expression profiles between diets was assessed using Principal Component Analysis (PCA). Differential expression of individual transcripts between diets was assessed using the limma R Bioconductor package (version 3.48.3) (Ritchie et al., 2015). Briefly, statistical significance was assessed using an empirical Bayes method (eBayes function) with an additive linear model accounting for diet and strain effect and adjusted P values were calculated by the Benjamini-Hochberg (BH) approach. Transcripts showing BH-adjusted P value below 0.05 and absolute Log₂ (Fold Change) above 0.5 were considered significantly associated with the effect of the diet.

Gene set enrichment analysis (GSEA) and Over-representation analysis (ORA)

Gene sets used in GSEA and ORA consisted of two parts: (1) the gene sets from the GO, KEGG, Hallmark, and Reactome databases were retrieved through the msigdbr R package (version 7.2.1) (Liberzon et al., 2011). (2) the gene signatures of mouse and human IBD were used as custom gene sets (Table 1).

GSEA was performed using clusterProfiler R package (version 3.10.1) (Yu et al., 2012) based on the log₂(Fold Change) ranking using parameters (nPerm = 100000, minGSSize = 30, maxGSSize = 5000, pvalueCutoff = 1). The gene sets with absolute NES higher than 1 and BH-adjusted P value lower than 0.05 were identified as the significantly enriched gene sets.

ORA analysis was also performed using clusterProfiler R package (version 3.10.1) (Yu et al., 2012) using parameters (minGSSize = 30, maxGSSize = 800). The gene sets with adjusted P value calculated by BH lower than 0.05 were identified as the significantly enriched gene sets.

Weighted gene correlation network analysis (WGCNA)

We used WGCNA R package (v1.51) (Langfelder and Horvath, 2008) to construct co-expression networks under CD and HFD, respectively. Firstly, the correlations between all pairs of gene across all BXDs fed with CD or HFD were calculated by Pearson correlation. Then, a best soft-thresholding power of 4 and 3 was chosen using pickSoftThreshold function with parameters (networkType = "signed hybrid", blockSize = 25000, corFnc = "bicor") for CD and HFD datasets in BXD colons separately. According to the calculated correlation coefficients, a network was constructed using parameters (networkType = "signed hybrid", minModuleSize = 30, reassignThreshold = 1e-6, mergeCutHeight = 0.15, maxBlockSize = 25000). The constructed co-expression gene modules were assigned color names and the module eigengenes were also identified for further analyses. To detect the preserved CD-modules in the co-expression modules under HFD, we defined gene modules under CD as custom genesets and performed ORA on each HFD-modules.

Transcription factor (TF) enrichment analysis

We first constructed a lognormal background distribution using the sequences of + 5kb region around the transcription starting site (TSS) of all genes and then downloaded the mouse HOCOMOCO-v10 (Kulakovskiy et al., 2018) motifs from R package motifDB to perform TF enrichment analyses using R package PWMenrich. The significantly enriched motifs (P value < 0.001) were selected and then ranked based on the percentage of enriched promoters.

Module Quantitative Trait Locus (ModQTL) mapping in the BXDs

We first downloaded genotype information of each BXD mice from GeneNetwork (https://gn1.genenetwork.org/webqtl/main.py?FormID=sharinginfo&GN_AccessionId=600) and generated the kinship matrix of BXD mice using the leave-one-chromosome-out (LOCO) method. We then used the eigengenes of each module as phenotype input to perform Module QTL (ModQTL) with the R package qtl2 (version 0.28) (Broman et al., 2019) and the threshold of each QTL mapping analysis was obtained from a permutation test with 10,000 repeats. The peaks of QTL were calculated by find_peaks function with parameter: prob=0.95.

The same methods were also applied to gene expression QTL mapping (eQTL) and the significance threshold of each gene was obtained from a permutation test with 1,000 repeats. The significant peaks overlapped with the location of their corresponding gene were identified as cis-eQTL.

Literature mining

To explore the inflammation related genes, we first used candidate gene names and keywords (“IBD”, “inflammatory bowel disease”, “Ulcerative colitis”, “Inflammation”, “Inflammatory”, “Crohn’s disease”) to search the title or abstract of associated literature using R package easyPubMed (version 2.13). Then, the genes involved in inflammation were confirmed by manual curation.

Genome-wide association study (GWAS) in UKBB

The phenotype data of inflamed Ulcerative colitis (Data-Field 131629, n = 6,459) and Crohn’s disease (Data-Field 131627, n = 3,358) were firstly downloaded from UKBB (Bycroft et al., 2018). 200,030 individuals with whole genome sequencing (WGS) (Halldorsson et al., 2022) in UK Biobank were selected and then the population of European descent (including with 1,173 patients with Crohn’s disease and 2,295 patients with UC) was extracted for further GWAS analyses. Control individuals (n = 143,194) were included based on the following criteria: (1) Individuals without non-inflamed colitis (Data-Field 131631), Crohn’s disease, and UC. (2) Individuals not taking any IBD-related medicine (Appendix 1 - Table 7).

WGS data provided by UK Biobank and used for GWAS were processed starting from pVCF files. We used REGENIE step1 to estimate population structure and then REGENIE step2 were applied to test associations between phenotypes and genetic variants and also included the following covariates in our model: the first 10 genetic principal components, age, sex, age:sex interaction, Body Mass Index (BMI), and smoking status. All data preparation and GWAS steps were run on DNAnexus.

Mendelian randomization (MR) analysis

eQTLs in sigmoid colon and transverse colon were selected in and their effect sizes obtained from the GTEx Portal on 2023-03-28 (v8, https://www.gtexportal.org/home/datasets, dbGaP Accession phs000424.v8.p2) (GTEx Consortium, 2013). No eQTLs were found for EPHA6 but 147 and 87 eQTLs were found for MUC4 in the sigmoid colon and transverse colon, respectively.

GWAS summary statistics for outcomes of interest, namely inflammatory bowel disease, ulcerative colitis, and Crohn’s disease, were obtained from the IEU OpenGWAS project (Elsworth et al., 2020), using the ieugwasr package (v0.1.5, https://mrcieu.github.io/ieugwasr/). As multiple GWAS exist for these diseases, we selected a representative one for each, prioritizing those with greater sample sizes and more cases, while still providing enough genetic variants for a large overlap with GTEx data (> 1 M SNPs). The selected GWAS were inflammatory bowel disease (ieu-a-31) (Liu et al., 2015), ulcerative colitis (ieu-a-32) (Liu et al., 2015), and non-cancer illness code, self-reported: Crohn’s disease (ukb-b-8210) (Elsworth et al., 2020).

For each outcome, MUC4 eQTLs which were present in the outcome GWAS were pruned for independence (more than 10 kbp away or r² < 0.01) using Plink v1.90b6.21 (Purcell et al., 2007) with the GTEx LD reference panel. In most cases this resulted in only a single eQTL being retained, with the exception of the sigmoid colon - Crohn’s disease combination which resulted in two independent eQTLs.

Mendelian randomization (MR) was performed using the TwoSampleMR R package (version 0.5.6) (Hemani, Tilling and Smith, 2017; Hemani et al., 2018). In the case with more than one eQTL, we used inverse-variance weighted MR, otherwise the Wald ratio. Because the magnitude of the normalized effect sizes provided by GTEx have no direct biological interpretation (https://gtexportal.org/home/faq#interpretEffectSize), the resulting causal effect estimates do not have an associated unit and cannot be translated into direct biological consequences. The direction (sign) of the effect remains interpretable.

Acknowledgements

We thank the Schoonjans’ and Auwerx’s lab members for technical assistance and discussions and Giacomo von Alvensleben for providing the GWAS analysis pipeline in human UKBB. The work in the JA laboratory was supported by grants from the Ecole Polytechnique Fédérale de Lausanne (EPFL), the European Research Council (ERC-AdG-787702), the Swiss National Science Foundation (SNSF 31003A_179435) and the Global Research Laboratory (GRL) National Research Foundation of Korea (NRF 2017K1A1A2013124). XL was supported by the China Scholarship Council (201906050019).

Author contributions

The study was conceived by XL, MBS and JA. EW, MBS and AB performed laboratory experiments. Data analyses were carried out by XL, AB, AR, JS and JP. XL and JA wrote the original manuscript. XL, MBS, JDM, GB, AP, KS and JA reviewed and edited the manuscript with contributions from all co-authors.

Competing interests

Authors declare no conflict of interest related to the work reported.

Data availability

The data that support the findings are available upon request to the corresponding authors (MBS and JA). The microarray data are available under the GEO numbers GSE225791. To review this dataset, please use this link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE225791 and the review token is exwtuucwpdqlzsr. Methods, materials and external resources are included in the Materials and Methods.

Table Supplements and Legends

Appendix 1 - Table 1. Co-expression networks under CD or HFD. Related to Figure 3

Appendix 1 - Table 2. Genes under QTL peak of module HFD_M28. Related to Figure 5

Appendix 1 - Table 3. Associations between genetic variants of EPHA6 and IBD. Related to Figure 5 and its figure supplement 1

Appendix 1 - Table 4. Associations between genetic variants of MUC4 and IBD. Related to Figure 5 and its figure supplement 1

Appendix 1 - Table 5. G-MAD result for Epha6 in mice and MUC4 in humans. Related to Figure 5

Appendix 1 - Table 6. MR result for MUC4.

Appendix 1 - Table 7. Medicine for human IBD.