The Effects of Sex and Diet on Physiology and Liver Gene Expression in Diversity Outbred Mice

Diversity Outbred, Nutrigenomics, QTL, Complex Traits 13 14 Karen L. Svenson 15 The Jackson Laboratory 16 600 Main St. 17 Bar Harbor, ME 04609 18 207-288-6342 19 [email protected] 20 21 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 3

22 Inter-individual variation in metabolic health and adiposity is driven by many factors. 23 Diet composition and genetic background and the interactions between these two 24 factors affect adiposity and related traits such as circulating cholesterol levels. In this 25 study, we fed 850 Diversity Outbred mice, half females and half males, with either a 26 standard chow diet or a high fat, high sucrose diet beginning at weaning and aged them 27 to 26 weeks. We measured clinical chemistry and body composition at early and late 28 time points during the study, and liver transcription at euthanasia. Males weighed more 29 than females and mice on a high fat diet generally weighed more than those on chow. 30 Many traits showed sex- or diet-specific changes as well as more complex sex by diet 31 interactions. We mapped both the physiological and molecular traits and found that the 32 genetic architecture of the physiological traits is complex, with many single locus 33 associations potentially being driven by more than one polymorphism. For liver 34 transcription, we find that local polymorphisms affect constitutive and sex-specific 35 transcription, but that the response to diet is not affected by local polymorphisms. We 36 identified two loci for circulating cholesterol levels. We performed mediation analysis by 37 mapping the physiological traits, given liver transcript abundance and propose several 38 genes that may be modifiers of the physiological traits. By including both physiological 39 and molecular traits in our analyses, we have created deeper phenotypic profiles to 40 identify additional significant contributors to complex metabolic outcomes such as 41 polygenic obesity. We make the phenotype, liver transcript and genotype data publicly 42 available as a resource for the research community. 43 44 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 4 45

46 47 Many factors affect the physiology and transcriptional landscape of individuals. Intrinsic 48 factors, such as sex and genetic background, play a role in shaping individuals. External 49 factors, such as diet and other environmental stimuli, also play a role in determining the 50 health and well-being of each person. However, wide variation in individual response to 51 diet is an enormous challenge to developing prevention and treatment strategies aimed 52 at reducing the incidence of obesity and metabolic disorders. The response to diet is 53 affected by both sex (G RIFFIN et al. 2016) and individual genetic background (S UHRE 54 and GIEGER 2012). Disparity in obesity prevalence among sexes has been ascribed to 55 both biological and sociocultural factors (KANTER and CABALLERO 2012). Sex differences 56 in fat gain and storage can lead to considerable differences in health outcome among 57 women and men (P OWER and S CHULKIN 2008). Additionally, sex-dependent single 58 nucleotide variants have been reported that underlie differential contributions to the 59 development of obesity (KVALOY et al. 2013; SALDANA-ALVAREZ et al. 2016). 60 61 In human populations, it is difficult to dissect the genetic basis for differential responses 62 to diet between the sexes due to differences in lifestyle and uncontrolled covariates. 63 While epidemiological models can estimate correlations between different traits, the 64 controlled conditions in mouse models allow us to apply randomization and factorial 65 designs to detect causal associations between obesity and physiological traits. In 66 mouse models, we can also control the genetic background of the mice, thereby 67 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 5 reducing another uncontrolled variable that influences the response to diets between 68 the sexes. Most mouse models of obesity use a single inbred strain genetically 69 engineered or experimentally manipulated to become obese (reviewed in (H ARIRI and 70 THIBAULT 2010)). However, genetic background is known to influence the response of 71 individuals to dietary fat (WANG et al. 2002; STOEHR et al. 2004; SU et al. 2008; LIN et al. 72 2013). In order to generalize our results from mice to humans, it is critical to include 73 structured genetic diversity in mouse models of dietary response. 74 75 Multi-parent advanced intercross (MAGIC) populations are powerful models for mapping 76 genetic modifiers of complex traits due to their high minor allele frequency and fine 77 mapping resolution (C HURCHILL et al. 2004; R AKSHIT et al. 2012; R AT GENOME et al. 78 2013; G ATTI et al. 2014). In MAGIC populations derived from known founders, the 79 haplotype structure of each sample genome can be reconstructed in terms of the 80 founder genomes (MOTT et al. 2000; GATTI et al. 2014). When the founders have been 81 fully sequenced, the founder sequences can be imputed onto the MAGIC genomes to 82 allow for whole genome association mapping (Y ALCIN et al. 2005), which improves the 83 ability to identify candidate genes that influence traits. Transcript profiling in a relevant 84 tissue adds another important dimension to genetic mapping studies and can be used to 85 perform mediation analysis on each significant genomic locus (C HICK et al. 2016). 86 There are several mouse MAGIC populations available, including the Northport 87 Heterogeneous Stock (V ALDAR et al. 2006), the Collaborative Cross (T HREADGILL and 88 CHURCHILL 2012), the Heterogeneous Stock/Collaborative Cross (IANCU et al. 2010) and 89 the Diversity Outbred (SVENSON et al. 2012). 90 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 6 91 In this study, we fed Diversity Outbred mice of both sexes either standard chow or a 92 high-fat/high-sucrose diet from weaning until approximately 6 months of age. We 93 measured a variety of physiological traits throughout the study and performed liver 94 transcriptional profiling at the end of the study. Here, we report on the differential effects 95 of diet on each sex in terms of physiological and transcriptional trait, provide interactive 96 viewers for the results and release the entire data set to the public through 97 supplemental materials and at http://do.jax.org. 98 99

100 101 Mice and husbandry: Diversity Outbred mice were obtained from The Jackson 102 Laboratory (Bar Harbor, ME). This study used five independent cohorts of 100-200 non-103 sibling DO mice from generations 4 to 11 (G4-G11) for a total of 850 animals, which 104 builds on an initial study that has previously been reported (S VENSON et al. 2012). In 105 each cohort, half the animals were from first litters in the respective generation and half 106 were from second litters. An equal number of females and males were included in each 107 set of animals received. Mice were housed at a density of five same-sex mice per pen in 108 pressurized, individually ventilated cages (Thoren #11 Duplex II; Thor en Caging 109 Systems, Hazelton, PA) with pine bedding (Crobb Box, Ellsworth, ME) and free access 110 to food (diets described below) and acidified water. Light cycle was 12h:12h light:dark, 111 beginning at 0600. All animal procedures were approved by the Animal Care and Use 112 Committee at The Jackson Laboratory (Animal Use Summary # 06006). 113 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 7 114 Phenotyping : Upon receipt, when mice were 3 weeks of age, equal numbers of each 115 sex were randomly assigned to chow (LabDiet 5K52, LabDiet, Scott Distributing, 116 Hudson, NH) or high fat, high sucrose feeding (Envigo Teklad TD.08811, Envigo, 117 Madison, WI) for the duration of the study protocol (26 weeks). Caloric content of the 118 high fat diet (HFD) was 45% fat, 40% carbohydrates and 15% protein. Tail biopsies 119 were taken at wean for DNA preparation. Weight was monitored weekly throughout the 120 study. At age 8 weeks mice began a pipeline of noninvasive phenotyping assays to 121 profile metabolic health (Table 1). Some modifications to the pipeline were made as the 122 number of cohorts progressed, such that all parameters were not measured in all mice. 123 Table 1 lists each phenotypic measurement and the number of mice tested. Data was 124 obtained from 846 mice and 154 traits were used for analysis. Clinical chemistries, 125 urinalysis and body composition assessments were performed at two time points in the 126 study to evaluate stability of traits under prolonged HFD. Hence, calculated traits 127 comparing first and second measures were generated as derived traits. Details about 128 blood collection and analysis and body composition by dual-energy x-ray 129 absorptiometry (DEXA) have been described previously for this pipeline (Svenson 2012 130 Genetics). Additional tests include body composition by qNMR (EchoMRI), 131 electrocardiogram, intraperitoneal glucose tolerance test (ipGTT), and evaluation of 132 chemokines by electrochemiluminescence. To quantitate lean and fat tissue and free 133 and total water, EchoMRI (EchoMRI, Houston, TX) without anesthesia was used, 134 providing three time points for evaluation of tissue composition during the study and 135 minimizing the need for anesthesia during the pipeline. Electrocardiography was 136 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 8 performed using the ECGenie™ (Mouse Specifics, Quincy, MA) system, whereby 137 unanesthetized mice are placed on a platform raised 18” above the laboratory bench 138 containing a lead plate. When animals contact the plate with any three paws the trace 139 begins. Fast Fourier analysis (AnonyMouse™ software v2.2; Mouse Specifics) defines 140 interval durations from which heart rate, variability and other features of cardiac 141 conduction can be assessed. To evaluate glucose clearance, mice were fasted 142 overnight (15 hours) and in the morning mice were weighed and a small blood sample 143 from a tail tip incision was used in the Abbott glucometer system to measure fasted 144 glucose (GTT time 0; t0). A glucose solution was then injected intraperitoneally at 2 mg 145 glucose/gram body weight and tail tip blood samples were obtained at 15, 30, 60, 120 146 and 180 minutes after injection. Plasma leptin, insulin, ghrelin and adiponectin were 147 measured from nonfasted mice using the Meso Scale Discovery™ 148 electrochemiluminescent assay detection system according to the manufacturer’s 149 protocols. We transformed all traits to ranked Z-scores before performing statistical 150 analyses. 151 152 Genotyping and Diplotype Reconstruction: DNA was prepared from tail biopsies and 153 genotyped using two versions of the Mouse Universal Genotyping Array (MUGA) 154 (MORGAN et al. 2016). We genotyped 531 samples on the MUGA and 293 samples on 155 the Megamuga (GeneSeek, Lincoln, NE). We used the intensities from each array to 156 infer the haplotype blocks in each DO genome using a hidden Markov model (HMM) 157 (GATTI et al. 2014). 158 159 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 9 Genotyping by RNA-sequence: Genotyping by RNA-sequence (GBRS) is a set of 160 software tools that reconstruct individual genomes of each sample in multi-parent 161 population (MPP) models by decoding known polymorphisms of founder strains from 162 RNA-Seq data without resorting to genotyping arrays. The new method is efficient since 163 it avoids maintaining hundreds of individualized genome indexes by aligning RNA-Seq 164 reads to a common pooled transcriptome of all founder strains a single time. Since our 165 reusable model parameters can be easily estimated from separate RNA-Seq data of 166 inbred founder strains or from simulations, we can quickly process each MPP sample 167 independently. The software package implements our alignment strategy and statistical 168 models and is freely available at https://github.com/churchill-lab/gbrs. We used the 169 GBRS haplotype reconstructions to fill in samples that failed to genotype due to low call 170 rates on the MUGA or Megamuga. 171 172 Merging Haplotype Reconstructions from Different Methods: The MUGA and 173 Megamuga have different numbers of markers (MUGA: 7,854, Megamuga: 77,642) and 174 the HMM produced diplotype probabilities only at each marker. In contrast, the GBRS 175

produced diplotype probabilities at each gene that was expressed in the liver. 176 In order to merge diplotype probabilities from the data, we interpolated both markers 177 grids to an evenly spaced 64,000 marker grid (0.0238 cM between markers). After 178 merging the diplotype reconstructions, we had a total of 835 samples. 179 180 Principal Component Analysis of Physiological Traits and Liver Transcription: We 181 retained 129 out of 160 physiological traits with < 50% missing data across samples. 182 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 10 The 24 traits with > 50% missing data were ACR1, ACR2, Adiponectin, BW.3, BW.27, 183 BW.28, BW.29, BW.30, fat.g.mri, free.h20, FRUC1, Ghrelin, GTT.AUC, GTT.t0, 184 GTT.t15, GTT.t30, GTT.t60, GTT.120, GTT.180, Lipase1, non.fast.Calcium, TBIL1, 185 TBIL2 and tot.h20 (see File S1 for abbreviations). We used the Probabilistic PCA 186

of the pcaMethods software (S TACKLIES et al. 2007) to impute missing data in 187 the remaining traits and calculated the first 10 principal components. 188 189 Physiological Traits Correlations: We regressed out the outbreeding generation, sex 190 and diet effects from each of the physiological traits and calculated the pairwise 191 Pearson correlation between all physiological traits. 192 193 Alignment, Quantification and Normalization of Liver Transcription Data: We 194 aligned reads from the DO liver data to pooled transcriptomes derived from the eight 195 DO founder strains by incorporating strain-specific SNPs, insertions and deletions into 196 the reference genome sequence. We quantified expected read counts using an 197 expectation maximization algorithm (EMASE, https://github.com/churchill-lab/emase) 198 (CHICK et al. 2016). We retained 12,067 genes with mean transcripts per million across 199 all samples greater than one. We normalized effective counts to the upper quartile value 200 and transformed them to rank normal scores. 201 202 Differential Expression and Gene Set Enrichment Analysis of Liver Transcript 203 Data: 204 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 11 We performed analysis of variance (ANOVA) on the normalized liver transcription data 205 to identify genes that were differentially expressed between sexes, diets or that had a 206 sex by diet interaction. We regressed the expression of each gene on generation and 207 litter, sex, diet and the sex by diet interaction. We adjusted the p-values using the 208 Benjamini & Hochberg false discovery rate (FDR)(BENJAMINI and HOCHBERG 1995). 209 210 We searched for Gene Ontology (GO) categories (A SHBURNER et al. 2000) that were 211 differentially expressed between sexes, diets or that had a sex by diet interaction using 212 the SAFE software package (B ARRY et al. 2005). We used the “t.Student” local statistic 213 to test for differential expression for each gene and the “Wilcoxon” global statistic to test 214 for differential enrichment between categories. For sex effects, we regressed diet from 215 each gene and then tested for the effect of sex. For diet effects, we regressed sex from 216 each gene and then tested for the effect of diet. For the sex by diet interaction, we 217 compared the reduced model with sex and diet to the full model containing sex, diet and 218 the sex by diet interaction. We determined the empirical p-value for each category using 219 10,000 permutations and retained GO categories with p-values ≤ 0.05. 220 221 Linkage Mapping: At each marker, we regressed each phenotype on generation, sex, 222 diet and the diplotype probabilities for each mouse and included an adjustment for 223 correlation between residuals due to kinship. We used the same model for liver 224 expression QTL mapping. We performed 5,000 permutations of a rankZ transformed 225 phenotype and selected significance thresholds from the empirical distribution of 226 maximum LOD scores. We estimated the founder allele effects using a Best Linear 227 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 12 Unbiased Predictor (BLUP) in which we fit a mixed-effects model at each marker that 228 shrinks the founder effects in proportion to the magnitude of the standard errors. We 229 used the qtl2 R package available at: https://github.com/rqtl . Full details of the linkage 230 mapping model are in (GATTI et al. 2014). 231 232 Association Mapping: We imputed the DO founder SNPs from the Sanger Mouse 233 Genomes Project (REL-1505) onto each founder haplotype block in the DO genomes. 234 We then regressed each phenotype on generation, sex, diet and the SNP probabilities 235 for each mouse and included an adjustment for correlation between residuals due to 236 kinship. Full details of the association mapping model are in (GATTI et al. 2014). 237 238 Mediation Analysis: For each physiological QTL peak with a genome-wide adjusted p-239 value above 0.05, we performed mediation analysis to identify candidate liver genes 240 that might be responsible for the peak (C HICK et al. 2016). We fit a null model by 241 regressing the phenotype on generation, sex, diet and the diplotype probabilities at the 242 markers with the maximum LOD score. We added the expression of each of the 12,067 243 genes to the model and recorded the drop in the LOD score compared to the null 244 model. We estimated the standard deviation of the LOD drops and report genes that 245 decreased the LOD score by more than 6 standard deviations. We refer to these 246 standardized LOD score drops as “Z-scores”. We used the intermediate R package 247 available at https://github.com/simecek/intermediate . 248 249 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 13 Data and Reagent Availability: J:DO mice are available for purchase from The 250 Jackson Laboratory (Strain # 009376) at https://www.jax.org/strain/009376. The liver 251 gene expression data is archived at the Short Read Archive under project number 252 PRJNA35625. The physiological phenotypes are described in File S1, the raw 253 phenotypes are in File S2 and the normalized phenotypes are in File S3. The genotype 254 data for all mice and the R data objects used in all analyses are available at 255 http://do.jax.org. We used the Sanger REL-1505 SNPs and structural variants (KEANE et 256 al. 2011) and the Ensembl build 82 transcripts (YATES et al. 2016). 257 258 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 14

259 260 Impact of Sex and Diet on Physiological Traits and Liver Transcription 261 262 We maintained mice of both sexes on either a chow diet (n = 449) or a high fat diet (n = 263 397) from wean to at least 26 weeks. We measured a range of physiological traits 264 throughout the study and measured several traits at two time points (File S1). We 265 calculated the principal components of the physiological traits and found that the mice 266 grouped by sex and diet (Figure 1A). Principal component (PC) 1 accounted for 29.2% 267 of the variance and is correlated with sex. PC2 accounts for 7.6% of the variance and is 268 correlated with diet. We also quantified liver transcription at 26 weeks in a subset of 478 269 mice. When we calculated the PCs for a subset of 12,067 genes, we found that the 270 samples also clustered by sex and diet (Figure 1B). PC1 and 2 accounted for 12.1% 271 and 10.9% of the variance, respectively. Mice on the chow diet formed tighter clusters 272 than mice on the high fat diet, reflecting larger variance in liver gene expression in mice 273 on the high fat diet. 274 275 Correlation between Physiological Traits: We calculated the pairwise Pearson 276 correlation of all traits after regressing out sex and diet and identified several clusters of 277 correlated traits (File S4). Most traits measured at two time points clustered near each 278 other, indicating that genetic background affects many traits throughout the mouse’s 279 lifespan. Body weight (BW) at all time points was positively correlated with other BW as 280 well as adiponectin, insulin, bone mineral density (BMD) and lean and fat tissue mass. 281 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 15 Cholesterol (CHOL) at 19 weeks was highly correlated with high-density lipoprotein 282 (HDLD, ρ = 0.95, p < 10-16),as expected for mice, triglycerides (TG, ρ = 0.40, p < 10-16), 283 glucose (ρ = 0.28, p = 1.7 x 10-16), non-esterified fatty acids (NEFA, ρ = 0.44, p < 10-16), 284 body weight ( ρ = 0.20, p < 6.6 x 10 -9) and circulating calcium (Ca, r = 0.50, p < 10 -16). 285 These correlations are in agreement with recent evidence that circulating calcium levels 286 are associated with worsening lipid profiles in humans (G ALLO et al. 2016) and that 287 coronary artery calcification is an independent risk factor for atherosclerosis and 288 cardiovascular disease (B UDOFF et al. 2007). The area under the curve of the glucose 289 tolerance test at 24 weeks (GTT) was positively correlated with BW at 24 weeks (ρ = 290 0.29, p = 4.70 x 10-5) and other time points, GLUC at 19 weeks (ρ = 0.30, p = 1.87 x 10-291 6), leptin at 8 weeks (ρ = 0.21, p = 2.62 x 10-3), indicating a connection between appetite 292 and circulating glucose levels. Leptin ( ρ = 0.80, p < 10 -16), insulin ( ρ = 0.42, p < 10 -16), 293 adiponectin (ρ = 0.40, p < 10-16) and % fat at both time points were positively correlated, 294 indicating a connection between appetite and adiposity. Glutamate dehydrogenase 295 (GLDH) at 19 weeks was positively correlated with and BW traits at ages over 15 weeks 296 (ρ = 0.212, p = 2.45 x 10 -9). While this observation may suggest that liver injury is 297 associated with increased weight, it is not correlated with increased fat tissue mass as 298 might be expected. It is likely that this correlation is an effect of aging and may be driven 299 by those animals that were fed HFD. The correlations are available as an interactive on-300 line tool at: http://churchill-lab.jax.org/www/Svenson850/corr.html. 301 302 Impact of Sex and Diet on Physiological Traits: We tested each trait for the effect of 303 sex, diet and a sex by diet interaction in order to identify the effects of each on the 304 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 16 physiological traits (File S5). This analysis stratified our results into four effect classes, 305 each demonstrated by examples in Figure 2. There were 12 traits for which no 306 difference was found between sexes or diet groups, including eosinophil counts (EOS) 307 and spleen weight (Figure 2A). Sex had an effect on 130 traits at a false discovery rate 308 (FDR) ≤ 0.05 (Figure 2B). Males had higher mean values for 101 of these traits, 309 including body weight, monocyte counts (MONO), neutrophil counts (NEUT), glucose 310 tolerance test (GTT.AUC), heart rate (HR), mean corpuscular volume (MCV), mean 311 platelet volume (MPV), phosphorous, platelet counts (PLT), red blood cell distribution 312 width (RDW) and cholesterol (CHOL). Females had higher mean values for 29 traits, 313 including mean corpuscular hemoglobin concentration (CHCM), hemoglobin (HGB), 314 ghrelin, %Fat and red blood cell counts (RBC). 315 316 Diet had an effect on 116 traits at an FDR ≤ 0.05. Mice on the HFD had higher mean 317 traits values for adiponectin, weight and fat traits, %Fat, GTT.AUC, HGB, insulin, leptin, 318 total bilirubin (TBIL), glutamate dehydrogenase (GLDH) and CHOL. Mice on the chow 319 diet had higher blood urea nitrogen (BUN), kidney weight, MONO, NEUT, PLT, 320 triglycerides (TG) and urine creatinine and glucose. 321 322 CHOL was higher in males as compared to females at both time points and on both 323 diets (File S5). At eight weeks, males on the chow diet (94.1 mg/dL) had CHOL levels 324 20% higher than females (78.4 mg/dL). At 19 weeks, CHOL levels on the chow diet 325 were similar to levels at eight weeks in males (94.2 mg/dL) and females (75.8 mg/dL). 326 The HFD increased CHOL levels at both time points. At eight weeks, males on the HFD 327 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 17 (126 mg/dL) had CHOL levels that were 20% higher than females (105 mg/dL). At 19 328 weeks, males (128 mg/dL) had CHOL levels that were 16.3% higher than females (110 329 mg/dL) (Figure 2C). At 19 weeks, the HFD increased CHOL by 45.1% compared to the 330 chow diet in females and by 35.9% in males. CHOL levels did not change greatly 331 between eight and 19 weeks. Female CHOL levels on the HFD increased by 4.7% from 332 105 mg/dL to 110 mg/dL and males increased by 1.6% from 126 mg/dL to 128 mg/dL. 333 Therefore, the increase in CHOL levels compared to chow values was established by 334 eight weeks in the HFD group and increased only minimally by the second time point. 335 336 There were 14 traits for which sex and diet showed a significant interaction (FDR ≤ 337 0.05, Figure 2D). BMD2 had one of the strongest sex by diet interactions, with mice on 338 the chow diet having similar BMD between the sexes, but males having higher BMD 339 than females on the HFD. This may be due to males gaining more weight and needing 340 stronger bones to carry the weight. This is consistent with higher BW in males and the 341 correlation of BW to BMD, such that male weight gain may require bone fortification to 342 support increased body mass. 343 344 Impact of Sex and Diet on Liver Transcription: We tested each of the 12,067 345 transcripts for the effects of sex, diet and a sex by diet interaction (File S6). We found 346 7,723 genes with sex effects (FDR ≤ 0.01) and 5,299 genes with significant diet effects 347 (FDR ≤ 0.01). We found 757 genes with a significant sex by diet interaction. In order to 348 interpret the functional relevance of these large gene lists, we searched for Gene 349 Ontology (GO) categories that were enriched in for each effect. We identified 212 GO 350 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 18 Biological Process (GO.BP) categories out of 2,570 in which genes were differentially 351 expressed between the sexes (p ≤ 0.05, File S7). Of those traits affected by sex, organ 352 regeneration (GO:00031100) was the most significant category and was higher in 353 males. This was followed by lipid metabolism (GO:0006629) and transport 354 (GO:0006869), which was higher in males. However, cholesterol metabolism 355 (GO:0008203, GO:0006695) was lower in males. Fatty acid metabolism (GO:0000038 & 356 GO:0070542) and beta-oxidation were higher in male mice along with catabolism of 357 triglycerides (GO:0019433). Of note, male mice had higher expression of genes 358 involved in unfolded protein responses (GO:1900103, GO:0072321), extracellular matrix 359 disassembly (GO:0022617), and fibroblast proliferation (GO:0048146). This suggests 360 that males experienced greater stress from unfolded protein responses, cellular 361 remodeling and proliferation. 362 363 We identified 212 GO.BP categories that contained genes that were differentially 364 expressed by diet (p ≤ 0.05, Supplemental File S11). For these, heat generation 365 (GO:0031649) and energy reserve metabolism (GO:0006112) were upregulated in mice 366 on HFD. Lipid metabolism (GO:0045834) was upregulated while lipid biosynthesis 367 (GO:0051055), including fatty acids (GO:0006633, GO:0042761) and phospholipids 368 (GO:0008654, GO:0015914) was down-regulated in mice fed the HFD. Lipid storage 369 (GO:0010888), cholesterol transport (GO:0030301) and gluconeogenesis 370 (GO:0045721) were all decreased in mice fed the HFD. Insulin secretion in response to 371 glucose stimulation was suppressed (GO:0061179) and glucose metabolism was 372 increased (GO:0010907). Overall, the HFD produced increases in lipid and energy 373 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 19 metabolism while decreasing lipid biosynthesis and storage, and perturbed glucose 374 homeostasis. 375 376 Physiological Trait Mapping 377 378 We mapped the physiological traits using each of three models: a model in which sex 379 and diet were additive covariates (additive model); on in which sex and diet were 380 additive covariates and sex interacted with genotype (sex-interactive model) and one in 381 which sex and diet were additive covariates and diet interacted with genotype (diet-382 interactive model). 383 384 Additive Model: We identified 82 additive QTL with genome-wide p-values ≤ 0.05 (File 385 S9). Of these, 39 were hematology traits, 23 were body weight or body composition 386 traits, 14 were clinical chemistry traits, 3 were electrocardiogram traits and 3 were 387 urinalysis traits. 388 389 Circulating cholesterol at 8 and 19 weeks (CHOL1 & CHOL2) had additive QTL on 390 chromosome 1 at 171.37 Mb with a LOD of 13.92 for CHOL1 and 13.57 for CHOL2 391 (Figure 3A). The pattern of founder allele effects at this locus was similar at both time 392 points (Figure 3B). The 129S1/SvImJ (129S1) and WSB/EiJ (WSB) alleles on the distal 393 end of chromosome 1 were associated with higher cholesterol levels. We performed 394 association mapping around the peak at 171.37 Mb and found that the most significant 395 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 20 SNPs (rs587286870 & rs580179709) had founder allele patterns for which the 129S1 396 and WSB strains carried the minor allele (Figure 3C). When we fit the association 397 mapping model again by including these SNPs as covariates, the maximum LOD score 398 in the region between 170 and 175 Mb decreased to 1.72, which is well below the 399 significance threshold. Both rs587286870 and rs580179709 are in an intron of prefoldin 400 2 ( Pfdn2), which is part of a molecular chaperone complex that stabilizes unfolded 401 proteins. It is unclear how Pfdn2 might impact circulating cholesterol levels. However, 402 we note that both of these SNPs are near a gene that is known to influence cholesterol 403 levels, apolipoprotein A-II ( Apoa2) located at 171.2 Mb. The 129S1 strain carries a 404 private alanine to valine substitution (rs8258226) that increases cholesterol levels 405 (WANG et al. 2004). The WSB strain carries a non-synonymous SNP (rs229811374) 406 that changes a serine to an asparagine and is located six nucleotides upstream of 407 rs8258226. If both of these SNPs increase cholesterol levels, this may explain the 408 increase in the 129S1 and WSB alleles effects at the Apoa2 locus (Figure 3D) and the 409 strong association with all SNPs where the minor allele occurs in both 129S1 and WSB. 410 The peak in this region is broad and may include more than one locus and allele. When 411 we regressed out the effects of the 129S1 and WSB alleles at the Apoa2 locus, the 412 maximum LOD of 5.35 occurred at 138.178 Mb on chromosome 1. We searched for 413 other genes expressed in the liver that might influence cholesterol levels by mediating 414 the peak with the expression of each gene (see Methods). We found six genes that 415 reduced the LOD score by greater than six standard deviations (i.e. had a Z-score < -6, 416 Figure 3E); inhibitor of kappaB kinase epsilon (Ikbke), peptidase M20 domain containing 417 1 ( Pm20d1), adenosine A1 receptor (Adora1), coagulation fa ctor X III, beta ( F13b), 418 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 21 complement factor H-related 1 (Cfhr1) and cathepsin E (Ctse). Of these, Cfhr1 had a Z-419 score of -18.6, which was lower than any other gene. The next lowest Z-score was -13 420 for Ctse. We tested whether Ctse would still reduce the LOD score by performing 421 mediation analysis with Cfhr1 in the model and found that Ctse still had a Z-score of -422 9.6. We found that F13b had a Z-score of -6.7 in this scan as well. The founder allele 423 effects for Cfhr1 (Figure 3F) show that mice carrying the PWK/PhJ (PWK) allele have 424 higher Cfhr1 levels and mice carrying the A/J allele have lower levels. For Ctse, the 425 129S1, WSB and NZO/HlLtJ (NZO) strains have higher Ctse expression (Figure 3G). 426 For F13b, the WSB allele is associated with lower F13b expression (Figure 3H). Cfhr1 427 is upregulated in the mouse retina in response to a different high fat diet (Z HENG et al. 428 2015). Ctse deficient mice fed a high fat diet showed hypercholesterolemia, reduced 429 body weight gain and impaired fat development compared to controls mice (K ADOWAKI 430 et al. 2014). F13b is part of the coagulation cascade and has not been previously 431 associated with cholesterol metabolism. 432 433 We found a peak for CHOL2 on chromosome 5 at 123.760 Mb for which the NZO allele 434 was associated with higher cholesterol levels (Figure 4A). When we mediated the QTL 435 peak at 123.76 Mb with the liver expression of each gene on chromosome 5, we found 436 that TRAF-type zinc finger domain containing 1 ( Trafd1) and scavenger receptor class 437 B1 (Scarb1) reduce the LOD score by more than six standard deviations (Figure 4B). 438 Trafd1 is associated with the regulation of innate immune responses and is not known 439 to have a function in cholesterol metabolism or transport. However, the pattern of allele 440 effects (Figure 4C) and the mediation score suggest that it may play an unknown role. 441 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 22 We repeated the mediation scan using Trafd1 expression as a covariate and found that 442 Scarb1 was the only gene that reduced the LOD score by more than 6 standard 443 deviations. DO mice carrying the NZO allele at the QTL had lower transcript levels of 444 Scarb1 (Figure 4D), which is consistent with the founder allele effects for CHOL2. 445 Scarb1 is the primary receptor for HDL-cholesterol uptake by the liver and steroidogenic 446 tissues and is vital for reverse cholesterol transport. Targeted mutations in Scarb1 lead 447 to abnormal lipoprotein metabolism and increased cholesterol levels (MOHR et al. 2004). 448 Liver-specific reduction in Scarb1 expression as a result of an ENU-induced point 449 mutation has also been reported, in which mice exhibit 70% higher plasma HDL-450 cholesterol levels due to reduced HDL selective uptake (S TYLIANOU et al. 2009). Scarb1 451 was proposed as a candidate gene for hypercholesterolemia in an intercross between 452 NZB/B1NJ and SM/J (PITMAN et al. 2002), but the authors found no sequence, mRNA or 453 protein differences. However, whole-genome sequencing of NZO has revealed a stop 454 gain mutation (rs233349576, Figure 4E & F) at residue 37 in ENSMUST00000137783. 455 This may produce an incomplete protein and may alter its function. 456 457 QTL that Interact with Diet: There were 12 QTL with p-values ≤ 0.05 for which 458 genotype interacted with diet (File S9), including 6 clinical chemistry traits, heart rate, 459 lymphocyte counts, urinary creatinine, body weight and body length. Cholesterol at 8 460 weeks (CHOL1) had a QTL that interacted with diet on chromosome 10 at 21.99 Mb 461 with a LOD of 10.6 (p ≤ 0.001, Figure 5A). Mice carrying the NOD allele on the HFD had 462 higher cholesterol levels. We mediated the QTL peak with all liver transcripts and found 463 that E030030I06Rik decreased the LOD by more than six standard deviations. 464 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 23 E030030I06Rik has a local eQTL on chromosome 10 at the same location. Association 465 mapping near the QTL peak produced significant associations with one gene, 466 Gm20125, which is a gene model with no known function. 467 468 Heart rate at 13 weeks (HR) had a QTL that interacted with diet on chromosome 6 at 469 125.63 MB with a LOD of 10.4 (p ≤ 0.001, Figure 5B). Mice carrying the A/J or 470 C57BL/6J allele on the HFD had higher HR than mice carrying the CAST allele. We did 471 not perform mediation analysis because we do not have heart transcript information on 472 these mice. Association mapping near the peak produced two SNPs (rs48596855, 473 rs38346309) in the introns of anoctamin 2 (Ano2), a calcium activated chloride channel, 474 a class of genes that may play a role in cardiac function (G UO et al. 2008; HARTZELL et 475 al. 2009). Another gene, potassium voltage-gated channel, shaker-related 1, (Kcna1), is 476 located 1 Mb distal from the peak SNPs and has been associated with changes in heart 477 rate (GLASSCOCK et al. 2010). 478 479 Glutamate dehydrogenase, a marker of liver injury, at 19 weeks (GLDH2) had a QTL 480 that interacted with diet on chromosome 9 at 92.19 Mb with a LOD of 11.86 (p ≤ 0.001, 481 Figure 5C). Mice carrying the A/J, NZO or PWK alleles on the HFD had higher GLDH 482 levels. There were no genes that had mediation Z-scores less than -6 near the QTL 483 peak. When we performed association mapping near the peak, we found four 484 transcripts with intronic SNPs and LOD scores over 7; Gm29478, 1700057G04Rik, 485 phospholipid scramblase 4 ( Plscr4), and procollagen lysine, 2-oxoglutarate 5-486 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 24 dioxygenase 2 ( Plod2). Both Plscr4 and Plod2 had local eQTL on chromosome 9. 487 Plscr4 is a membrane protein that is involved in the organization of phospholipids and 488 interacts with the CD4 receptor of T lymphocytes (PY et al. 2009). It is up-regulated with 489 HFD feeding (SONG et al. 2012). Plod2 hydroxylates lysine residues and is involved in 490 remodeling of the extracellular matrix (GILKES et al. 2013) and fibrosis (VAN DER SLOT et 491 al. 2003). Under hypoxic conditions, Plod2 is expressed in hepatocellular carcinoma 492 and is correlated with tumor size and macroscopic intrahepatic metastasis. It is a 493 prognostic factor for disease-free survival (NODA et al. 2012). 494 495 QTL that Interact with Sex: There were 9 QTL with p-values ≤ 0.05 for which genotype 496 interacted with sex (File S9), including 5 clinical chemistry traits and 4 hematology traits. 497 Blood urea nitrogen at 19 weeks (BUN2) had a QTL that interacted with sex on 498 chromosome 10 at 95.85 Mb with a LOD of 11.7 (p ≤ 0.001). Males carrying the 129S1, 499 C57BL/6J and WSB alleles were associated with higher BUN and females carrying the 500 PWK allele with lower BUN. We did not perform mediation analysis because we do not 501 have kidney transcript information on these mice. Total bilirubin at 19 weeks (TBIL2) 502 had a QTL that interacted with sex on chromosome 19 at 14.89 Mb with a LOD of 11.1 503 (p ≤ 0.001, Figure 5D). Females carrying the NZO allele had higher bilirubin and males 504 carrying the CAST allele had lower bilirubin. Mediation analys is using liver gene 505 expression did not reveal a candidate gene. The most significant SNPs in the 506 association mapping on chromosome 19 were near five transcripts: Gm8630, 507 Gm31441, Gm37997, Gm26026 and transducin-like enhancer of split 4 ( Tle4). None of 508 the gene models had a QTL. Tle4 is a transcriptional corepressor factor that regulates 509 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 25 mouse hematopoiesis and bone development (W HEAT et al. 2014), and has also been 510 used for histological application as a podocyte nuclear marker in glomeruli 511 (VENKATAREDDY et al. 2014). 512 513 Liver Expression QTL Mapping: We performed linkage mapping on 12,067 liver 514 genes and identified additive QTL (File S10), QTL that interact with sex (File S11) and 515 QTL that interact with diet (File S12). We mapped local and distant eQTL using an 516 additive linkage model, and two models in which sex or diet interacts with genotype. We 517 plotted the location of significant QTL versus gene location for each model and found 518 that the additive and sex-interactive models produced local eQTL (Figure 6A & B) and 519 the diet-interactive model did not (Figure 6C). At a significance threshold of 0.05, we 520 found that 8,127 local eQTL out of 9,754 total (83.3%) in the additive model, 332 out of 521 532 (62.4%) in the sex-interactive model and 23 out of 219 (10.5%) in the diet-522 interactive model. We have provided an on-line visualization tool at http://churchill-523 lab.jax.org/qtl/svenson/DO478/. 524 525 XO Females: XO females have been previously reported in high numbers in DO mice 526 (CHESLER et al. 2016). In order to search for XO females, we plotted the liver expression 527 of Xist as an indicator of X chromosome gene expression versus Ddx3y as a marker of 528 Y chromosome gene expression (Figure 7). As expected, females that were dizygous 529 for the X chromosome had high Xist expression and low Ddx3y expression while males 530 had high Ddx3y expression and lower Xist expression. There were two females (out of 531 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 26 244, 0.82%) that had low Xist expression, consistent with hemizygosity on the X 532 chromosome, and low Ddx3y expression and these samples are XO females. In 533 humans, Turner Syndrome describes females with the XO genotype and is 534 characterized by short stature, a propensity for ovarian dysfunction and infertility, and 535 heart defects. The two XO females in our study, one fed chow and the other fed the 536 HFD, had very different phenotypes and were not outliers for any particular trait. One of 537 them, however, was especially resistant to weight gain on HFD, gaining only 5 grams of 538 body weight over the course of the study, and died 4 weeks before the scheduled end of 539 the study. 540 541 542

543 544 Multi-parent populations are excellent tools for studying the effects of genetic diversity 545 on phenotypic variation because they offer increased genetic diversity, high minor allele 546 frequencies and fine recombination block structure. The large number of variants leads 547 to perturbation of genes throughout the genome and the high minor allele frequency 548 produces high power to detect the effects of these polymorphisms. The fine 549 recombination block structure leads to fine mapping resolution to identify loci that 550 contain a manageable number of candidate genes. These loci may contain genetic 551 variants that alter either protein structure or transcript levels or both. In fact, many loci in 552 multi-parent crosses may be caused by more than one genetic variant and 553 disentangling the signal from these different alleles is a complex process. For the 554 cholesterol loci on chromosomes 1 and 5, we combined association mapping using 555 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 27 imputed SNPs with mediation analysis using liver transcripts and we identified 556 candidate genes using both missense SNPs and transcript levels. This approach is 557 broadly applicable in the DO and other multi-parent populations, but requires 558 transcriptional profiling in a relevant tissue. 559 560 CHOL levels were associated with two loci on chromosomes 1 and 5. Variation at these 561 loci affected mice of both sexes and on both diets. In contrast, CHOL at eight weeks 562 had a QTL on chromosome 10 with effects that were modified by diet. The HFD 563 increased CHOL levels by 35.9% in males and 45.1% in females, indicating that diet 564 increases CHOL levels in most mice. However, DO mice carrying the NOD allele on 565 chromosome 10 have higher CHOL levels than mice carrying other alleles at the same 566 locus. The locus covered several Mb and may contain more than one polymorphism 567 that affects CHOL levels, independent of diet. 568 569 Association mapping on chromosome 1 produced a set of SNPs with high LOD scores 570 for which 129S1 and WSB contributed the minor allele. Initially, we expected to find 571 SNPs with this allele pattern that alter the protein structure or expression levels of some 572 gene. However, in this case, we believe that two closely located SNPs in Apoa2, each 573 of which is private to either 129S1 or WSB, are influencing CHOL levels. The SNP in 574 129S1 (rs8258226) has been shown to increase CHOL levels and we hypothesize that 575 the SNP in WSB (rs229811374), which is 2 residues away from rs8258226, may also 576 increase CHOL levels. This highlights the complexity of analyses in multi-parent 577 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 28 crosses. Had we not known of any candidate genes in this region, we may have been 578 led to consider other genes that are unrelated to CHOL metabolism. 579 580 Mediation analysis identified several candidate genes on chromosome 1. By its nature, 581 mediation analysis is a hypothesis generating analysis. All of the genes, Cfhr1, Ctse 582 and F13b are within 10 Mb of each other. Of these, Ctse is the only gene that had been 583 previously associated with hypercholesterolemia (Z HENG et al. 2015). The situation is 584 similar for the CHOL peak on chromosome 5. Scarb1 has been associated with 585 differences in CHOL levels, but Trafd1 is a new candidate gene that may have effects of 586 CHOL independent of Scarb1. These findings suggest that some strong associations 587 that appear in MAGIC populations may be due to multiple, tightly linked polymorphisms. 588 If this is true, more sophisticated, multi-locus approaches to candidate gene selection 589 will be needed to find causal genes. In this study, we suggest that both association 590 mapping and mediation analysis with transcript levels improve the reliability of candidate 591 gene selection because it allows investigators to search for polymorphisms that affect 592 protein structure or transcript levels. 593 594 When we performed liver eQTL mapping, we found local eQTL for both the additive 595 model and the sex-interactive model. This suggests that local polymorphisms in 596 regulatory elements modulate transcript levels constitutively in both sexes and 597 differentially by sex. However, when we mapped liver transcript levels using a model in 598 which diet interacts with genotype, we found very few local eQTL. This suggests that 599 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 29 the response to diet is less influenced by the interaction of diet with local polymorphisms 600 affecting transcript levels, we speculate that distant loci are acting through non-601 transcriptional mechanisms such as interactions between proteins (CHICK et al. 2016). 602 603 Our analysis of this large study using DO mice includes a novel multi-tiered approach 604 that has identified plausible candidate genes underlying physiological traits, has 605 extended to considering effects of transcrip tion, and provides compelling evidence for 606 further investigation of the role of novel candidates in driving metabolic traits. We 607 present an analysis of the complex interplay of sex and diet and how these factors 608 influence important inter-individual variation in outcome. We found that diet increases 609 many traits related to body size and composition. Interestingly, we observed that a high 610 fat, high sucrose diet increased the variance of liver gene expression, suggesting that 611 genotype influences the range of responses to diet. We mapped loci for multiple traits 612 but focused on CHOL to demonstrate the complexity of the underlying genetic loci. We 613 note that both association mapping and mediation analysis using liver transcript data 614 help to dissect the causal alleles underlying mapping peaks. Finally, we observe that 615 changes in liver transcription in response to diet are not primarily altered by local 616 genetic variants. This suggests that other molecular measurements, such as protein or 617 metabolite levels may be more useful in determining the effects of diet on organisms. 618 We have made the phenotype and genotype data fully available to the public and have 619 released an interactive viewer to allow the reader to explore the liver expression QTL 620 data. 621 622 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 30

623 This work was funded by P50GM076468 and R01GM070683 to G.A.C. from the 624 National Institute of General Medical Sciences. 625 626 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 31 FIGURES 627 628 Figure1. Principal component analysis (PCA) of physiological and liver transcription629 traits. (A) PCA plot of the first and second principal components of the physiological630 traits. Each point represents one mouse. Females are red; males are blue; mice on the631 chow diet are shown with open symbols; mice on the HFD with closed symbols. Dashed632 lines are ellipses from bivariate Gaussian di stributions fit over each of the four sex and633 diet groups. (B) PCA plot of the first and second principal components of the liver634 transcription traits. 635 636 637 on al he ed nd er certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 32 638 Figure 2. Physiological traits in the DO. Each plots shows female (circles ) and male639 (squares) mice on chow (open symbols) and HFD (solid symbols). The open symbols640 connected by dashed lines show the chow diet group means. The closed symbols641 connected by solid lines show the HFD group means. (A) Spleen weight did not vary642 between sexes or diet groups. (B) Lean tissue mass at 12 weeks differed by sex and643 was not affected by diet. (C) Cholesterol at 19 weeks had both a significant difference644 le ls ls ry nd ce certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 33 between sexes and diets. (D) Bone mineral density at 21 weeks had one of the most 645 significant sex by diet interactions, showing a greater response in males to HFD 646 647 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 34 648 Figure 3. Analysis of cholesterol QTL on chromosome 1. (A) Genome scan of 649 cholesterol at 19 weeks shows peaks on chromosomes 1 and 5. The horizontal axis 650 shows the mouse genome and the vertical axis shows the LOD score. Red and yellow 651 lines are the p = 0.05 and 0.2 significance levels, respectively. The horizontal axis in 652 panels B through H shows the location in Mb on chromosome 1. (B) DO founder allele 653 effects on chromosome 1 for cholesterol at 19 weeks (CHOL2). Each colored line 654 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 35 represents the estimated effect of one of the founder alleles along the chromosome. (C) 655 Association mapping near the peak on chromosome 1 shows that the SNPs for which 656 129S1 and WSB contribute the minor allele have the highest LOD scores (red 657 triangles). (D) Genes in the region one chromosome 1 shown in panel C. Apoa2 and 658 Pfdn2 are colored in red and are mentioned in the text. (E) Mediation analysis shows 659 that six genes reduce the LOD score by more than 6 standard deviations. The vertical 660 axis shows the Z-score (scaled LOD across all genes). Each point represent the Z-661 score (standardized LOD score reduction) for one gene in the mediation analysis. The 662 genes are located along the horizontal axis. The red line shows Z = -6. DO founder 663 allele effects for liver expression of (F) Cfhr1, (G) Ctse and (H) F13b on chromosome 1. 664 665 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 36 666 Figure 4. Cholesterol QTL at 19 weeks (CHOL2) on chromosome 5 suggest Trafd1 and 667 Scarb1 as a candidate gene for circulating cholesterol levels. (A) Founder allele effects 668 for CHOL2 show that the NZO allele at 123.7 Mb is associated with higher cholesterol 669 levels. Each colored line represents the estimated effects of one founder allele. (B) 670 Mediation analysis of the CHOL QTL on Chr 5 using liver transcripts as mediators. Each 671 dot represents the Z-score of the LOD decrease after including one gene in the 672 mapping model. The red line is the Z = -6 threshold. (C) Founder allele effects for liver 673 Trafd1 transcript levels have a similar pattern of allele effects as CHOL2. (D) Founder 674 allele effects for liver Scarb1 transcript levels show that DO mice with the NZO allele on 675 chromosome 5 at 123.7 Mb have lower levels of Scarb1. (E & F) Association mapping in 676 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 37 the interval near the QTL identifies a single SNP (rs233349576) that introduces a stop 677 codon into a Scarb1 transcript. 678 679 680 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 38 681 Figure 5. Sex- and diet-interactive QTL plots. Each plot has two panels. The top panel682 shows the additive LOD (dashed line) and the interactive LOD (solid line). The lower683 panel shows the difference between the interactive and additive founder allele effects,684 with standard errors in light shading. Each founder is represented by a separate color :685 A/J yellow, C57BL/6J grey, 129S1 /SvImJ pink, NOD/ShiLtJ cyan, NZO/HlLtJ blue,686 CAST/EiJ green, PWK/PhJ red and WSB/EiJ purple. (A) Cholesterol at eight weeks had687 el er ts, : e, ad certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 39 a diet-interactive QTL on chromosome 10. (B) Heart rate had a diet-interactive QTL on 688 chromosome 6. (C) Glutamate dehydrogenase at 19 weeks had a diet-interactive QTL 689 on chromosome 9. (D) Bilirubin at 19 weeks had a sex-interactive QTL on chromosome 690 10. 691 692 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 40 693 Figure 6. Liver expression QTL maps generated fr om an (A) additive model (File S10) ,694 (B) a model in which sex interacts with genotype (File S11), or (C) a model in which diet695 interacts with genotype (File S12). Each dot represents the location of a QTL peak for696 one gene above the p = 0.05 threshold. Each panel plots the QTL position on the697 horizontal axis and the transcript position on the vertical axis. 698 699 , iet for he certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 41 700 Figure 7. XO females in the DO population. For each mouse, we plotted the 701 untransformed expression of Xist versus Ddx3y and found two XO females (in red). 702 Females have high Xist expression and low Ddx3y expression. Males have low Xist 703 expression and high Ddx3y expression. 704 705 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 42

706 A SHBURNER , M ., C. A. B ALL , J. A . B LAKE , D. B OTSTEIN , H. B UTLER e t a l . , 2000 Gen e ontology: too l for the 707 unification of biology. The G ene On tolog y Consortium. Na t Gen et 25: 25-29. 708 B ARRY , W . T., A. B . N OBEL a n d F . A . W RIGHT , 2005 Significance analysis of functi onal categori es in gen e 709 expr ession studi es: a str uctur ed permu ta tion appr oach. Bi oinformatics 21: 1943-1949. 710 B ENJAMINI , Y., and Y. H OCHBERG , 1995 Controlling the F alse Discovery Ra te - a Practical and Powerfu l 711 Approach to M ultipl e Testing . Jou rnal o f the Royal S ta tistical S ociety S eries B- Metho dological 712 57: 289-300. 713 B UDOFF , M. J ., L. J . S HAW , S. T. L IU , S. R. W EINSTEIN , T. P. M OSLER e t a l . , 2007 Long-term prognosis 714 associated with coron ary calcificatio n: o bservations from a r egistry of 25,253 pa t ients. J Am Coll 715 Cardiol 49: 1860-1870. 716 C HESLER , E. J ., D. M. G ATTI , A. P. M ORGAN , M . S TROBEL , L. T REPANIER et al. , 2016 Diversity Out bred Mice a t 717 21: Maintai ning Allelic Varia tion in the Fa ce of Selection . G3 (Bet h esda) 6: 3893-3902. 718 C HICK , J . M . , S . C . M UNGER , P. S IMECEK , E . L. H UTTLIN , K. C HOI e t a l . , 2016 Definin g the conse quences of 719 genetic varia tion on a p rot eome-wide scale. N atu re 534: 500-505. 720 C HURCHILL , G . A . , D . C . A IREY , H . A LLAYEE , J. M . A NGEL , A. D. A TTIE et a l. , 2004 The Collabora tive Cross, a 721 community resource for the gen etic anal ysis of complex traits. Nat ure gen etics 36: 1133-1137. 722 G ALLO , L. , M. C. F ANIELLO , G . C ANINO , C. T RIPOLINO , A. G NASSO e t a l . , 2016 Serum Calcium Increas e 723 Correla tes Wit h Worsen ing of Lipid Pro file: An Observa tional S tudy on a Large Cohort Fro m 724 South I taly. Me dicine (Baltimo re) 95: e2774. 725 G ATTI , D. M ., K. L. S VENSON , A. S HABALIN , L . Y. W U , W. V ALDAR et al. , 2014 Quan tita t ive trait locus mapping 726

for diversi ty outbr ed mice. G3 ( Beth esda) 4: 1623-1633. 727 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 43 G ILKES , D . M . , S . B AJPAI , P. C HATURVEDI , D. W IRTZ a n d G . L . S EMENZA , 2013 Hypoxia-i nducible factor 1 (HIF-1) 728 promot es e xtr acellula r mat rix remod elin g under hypo xic condi tions by inducing P4HA1, P4HA2, 729 and PLOD2 expression in fibrobl asts. J Bi ol Chem 288: 10819-10829. 730 G LASSCOCK , E., J . W. Y OO , T . T. C HEN , T . L . K LASSEN a n d J . L . N OEBELS , 2010 Kv 1.1 po tassium chan nel 731 deficiency reveals br ain-driven cardi ac dysfunction as a candida te mecha nism for sudden 732 unexplai ned de ath in ep ilepsy. J Neur osci 30: 5167-5175. 733 G RIFFIN , C., N . L ANZETTA , L . E TER a nd K. S INGER , 2016 Sexually dimor phic myeloid inflammatory an d 734 metabolic respons es to die t-induced ob e sity. Am J Physiol Regul In tegr Comp Physiol 311: R211-735 216. 736 G UO , D. , L. Y OUNG , C. P ATEL , Z. J IAO , Y. W U et al. , 2008 Calcium-activated chlorid e current con tribu tes t o 737 action po ten tial al te rnati ons in left ve ntricular hypert rophy ra bbit . Am J Phy siol Hear t Circ 738 Physiol 295: H97-H104. 739 H ARIRI , N., and L. T HIBAULT , 2010 High-fat diet-induced ob esity in animal mode ls. Nutr R es Rev 23: 270-740 299. 741 H ARTZELL , H. C. , K. Y U , Q . X IAO , L. T. C HIEN and Z. Q U , 2009 Anoct amin/TMEM16 fa mily members a re Ca2 +-742 activated Cl- channels . J Physiol 587: 212 7-2139. 743 I ANCU , O . D . , P . D ARAKJIAN , N. A. W ALTER , B . M ALMANGER , D. O BERBECK e t a l . , 2010 Genet ic diversity and 744 striat al gen e ne tworks: focus on the h ete rogene ous stock-collabo rative cross (HS-CC ) mouse. 745 BMC genomics 11: 585. 746 K ADOWAKI , T., M. A. K IDO , J . H ATAKEYAMA , K. O KAMOTO , T . T SUKUBA e t a l . , 2014 Defective a dipose tissu e 747 developmen t associa ted with h epa tom egaly in cathepsin E-deficie nt mice fed a high-fat diet . 748 Biochem Biophys Res Commun 446: 212- 217. 749 K ANTER , R., an d B. C ABALLERO , 2012 Gl obal gender dispari ties in ob esity: a r eview. A dv Nutr 3: 491-498. 750 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 44 K EANE , T . M . , L . G OODSTADT , P. D ANECEK , M. A. W HITE , K. W ONG et al. , 2011 Mous e genomic variati on an d 751 its effect on pheno types and gen e regula tion. N atu re 477: 289-294. 752 K VALOY , K., B . K ULLE , P. R OMUNDSTAD and T. L. H OLMEN , 2013 Se x-specific effects o f weight-affecting gen e 753 variants in a life course p erspec tive--The HUNT Study, No rway. In t J Ob es (Lond) 37: 1221-1229. 754 L IN , C., M. L. T HEODORIDES , A . H. M C D ANIEL , M. G. T ORDOFF , Q . Z HANG et al. , 2013 QTL analysis of dietary 755 obesity in C57BL/6byj X 129P3 /J F2 mice: diet- and sex-d epend ent effects . PLoS One 8: e68776. 756 M OHR , M ., M. K LEMPT , B . R ATHKOLB , M. H . DE A NGELIS , E. W OLF et al. , 2004 Hyperc holest erolemi a in ENU-757 induced mouse mutan ts. J Lipid Res 45: 2132-2137. 758 M ORGAN , A. P., C. P. F U , C. Y . K AO , C. E. W ELSH , J. P. D IDION et al. , 2016 The Mo u se Universal Gen otyping 759 Array: Fr om Substra ins to Subspeci es. G3 (Bethesda) 6: 263-279. 760 M OTT , R ., C. J . T ALBOT , M. G . T URRI , A . C. C OLLINS and J . F LINT , 2000 A me thod for fin e mapping quan tit ativ e 761 trai t loci in outbr ed animal st ocks. Proc N atl Acad Sci U S A 97: 12649-12654. 762 N ODA , T., H. Y AMAMOTO , I . T AKEMASA , D. Y AMADA , M. U EMURA et al. , 2012 PLOD2 in duced under hypo xia is 763 a novel prognostic factor for hep atoc ellu lar carcinoma afte r curative r esectio n. Liver Int 32: 110-764 118. 765 P ITMAN , W . A . , R . K ORSTANJE , G. A. C HURCHILL , E. N ICODEME , J. J. A LBERS et al. , 2002 Quantit ative t rait locu s 766 mapping of genes tha t r egulat e HDL c holest erol in SM/J and NZB/B1NJ inbre d mice. Physiol 767 Genomics 9: 93-102. 768 P OWER , M . L . , a n d J . S CHULKIN , 2008 Se x differences in fat s torag e, fat me tabolis m, and the heal th risks 769 from obesity: possible evolution ary origi ns. Br J Nu tr 99: 931-940. 770 P Y , B., S . B ASMACIOGULLARI , J. B OUCHET , M. Z ARKA , I. C. M OURA et al. , 2009 The p hos pholipid scramblas es 1 771 and 4 a re c ellular rec epto rs for the secr e tory le ukocyte p rot ease inhibit or and in t eract with CD4 772 at the plasma membran e. PLoS One 4: e 5006. 773 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 45 R AKSHIT , S., A . R AKSHIT a n d J . V . P ATIL , 2012 Multipar ent in tercr oss populati ons in analysis of quantita tive 774 trai ts. J Gene t 91: 111-117. 775 R AT G ENOME , S., C. M APPING , A. B AUD , R. H ERMSEN , V. G URYEV et al. , 2013 Combined sequence-b ased and 776 genetic mapping analysis of complex trai ts in outbr ed ra ts. N at G ene t 45: 767-775. 777 S ALDANA -A LVAREZ , Y., M . G . S ALAS -M ARTINEZ , H. G ARCIA -O RTIZ , A . L UCKIE -D UQUE , G. G ARCIA -C ARDENAS e t a l . , 778 2016 Gender-Dep enden t Associatio n of FTO Polymorphisms with Body Mass Ind ex in Mexicans . 779 PLoS One 11: e0145984. 780 S ONG , Y. B., Y. R. A N , S. J . K IM , H. W. P ARK , J. W. J UNG e t a l . , 2012 Lipid metabol ic effect of Korean red 781 ginseng extr act in mice fed on a high-fat diet. J Sci Food Agric 92: 388-396. 782 S TACKLIES , W., H . R EDESTIG , M . S CHOLZ , D. W ALTHER a n d J . S ELBIG , 2007 pcaMe thods--a bioconduc to r 783 package providing PCA methods for inco mplete da ta . Bioinforma tics 23: 1164-11 67. 784 S TOEHR , J . P., J . E. B YERS , S. M . C LEE , H. L AN , I. V. B ORONENKOV e t a l . , 2004 Identifica tion of major 785 quanti tative trai t loci contr olling body weight variat ion in ob/ob mice. Diabe tes 53: 245-249. 786 S TYLIANOU , I. M. , K. L. S VENSON , S. K. V AN O RMAN , Y. L ANGLE , J. S. M ILLAR e t a l . , 20 09 Novel ENU-induced 787 point muta tion in scavenger recep tor cl ass B, member 1, results in liver specific loss of S CARB1 788 prote in. PLoS One 4: e6521. 789 S U , Z., R. K ORSTANJE , S . W . T SAIH a n d B . P AIGEN , 2008 Candida te genes fo r ob esity reve aled from a 790 C57BL/6J x 129S1/SvImJ intercross. In t J Obes (Lond) 32: 1180-1189. 791 S UHRE , K., a nd C. G IEGER , 2012 G ene tic va riation in met abolic phenotyp es: stu dy designs and applica tions . 792 Nat R ev Gene t 13: 759-769. 793 S VENSON , K. L., D. M. G ATTI , W. V ALDAR , C. E. W ELSH , R. C HENG et al. , 2012 High-resol ution gene tic mapping 794 using the Mouse Diversity ou tbr ed popul ation . Gen etics 190: 437-447. 795 T HREADGILL , D. W., and G . A. C HURCHILL , 20 12 Ten years of the collabo rative cr oss. G 3 2: 153-156. 796 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint 46 V ALDAR , W., L. C. S OLBERG , D. G AUGUIER , S . B URNETT , P. K LENERMAN e t a l . , 2006 Genome-wide gen eti c 797 association of comple x trai ts in he terog e neous stock mice. N atur e gene tics 38: 879-887. 798 VAN DER S LOT , A . J . , A . M . Z UURMOND , A. F. B ARDOEL , C. W IJMENGA , H. E. P RUIJS et al. , 2003 Identificatio n of 799 PLOD2 as telopeptide lysyl hydroxylase , an important en zyme in fibrosis. J Biol Chem 278: 800 40967-40972. 801 V ENKATAREDDY , M. , S. W ANG , Y. Y ANG , S . P ATEL , L. W ICKMAN et al. , 2014 Estima ting podocyte numb er an d 802 density using a single histologic sectio n. J Am Soc Nephrol 25: 1118-1129. 803 W ANG , H., L. H. S TORLIEN and X. F . H UANG , 2002 Effects of dietary fat types on bo dy fatness, leptin , and 804 ARC lepti n r ecept or, NPY, and AgRP m RNA exp ression . Am J Physiol Endoc rin ol Me tab 282: 805 E1352-1359. 806 W ANG , X., R. K ORSTANJE , D. H IGGINS and B. P AIGEN , 2004 Haplotype an alysis in multi ple crosses to id entify a 807 QTL gene. Gen ome Res 14: 1767-1772. 808 W HEAT , J . C., D. S. K RAUSE , T. H. S HIN , X. C HEN , J. W ANG e t a l . , 2014 The core pressor Tle4 is a novel 809 regulat or of murine hema topo iesis and b one developm ent . PLoS One 9: e105557. 810 Y ALCIN , B., J . F LINT a n d R . M OTT , 2005 U sing progenito r str ain informati on to id entify quanti ta tive tr ai t 811 nucleotid es in outb red mice. Gen etics 171: 673-681. 812 Y ATES , A ., W . A KANNI , M . R. A MODE , D. B ARRELL , K. B ILLIS et al. , 2016 Ensembl 2016. Nucleic Acids Res 44: 813 D710-716. 814 Z HENG , W., N. M AST , A . S AADANE a n d I . A . P IKULEVA , 2015 Pathways of chol est erol homeostasis in mous e 815 retin a responsive to die tary and pha rma cologic trea tments . J Lipid Res 56: 81-97. 816 817 818 certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was notthis version posted January 5, 2017. ; https://doi.org/10.1101/098657doi: bioRxiv preprint