Genome wide association study of behavioral, physiological and gene expression traits in a multigenerational mouse intercross

Genome wide association analyses ( GWAS ) in model organisms have numerous advantages compared 2 to human GWAS, including the ability to use populations with well-defined genetic diversity, the ability to 3 collect tissue for gene expression analysis and the ability to perform experimental manipulations. We examined behavioral, physiological, and gene expression traits in 1,063 male and female mice from a 5 50-generation intercross between two inbred strains (LG/J and SM/J). We used genotyping by 6 sequencing in conjunction with whole genome sequence data from the two founder strains to obtain 7 genotypes at 4.3 million SNPs. As expected, all alleles were common (mean MAF=0.35) and linkage 8 disequilibrium degraded rapidly, providing excellent power and sub-megabase mapping precision. We 9 identified 126 genome-wide significant loci for 50 traits and integrated this information with 7,081 cis - 10 eQTLs and 1,476 trans -eQTLs identified in hippocampus, striatum and prefrontal cortex. We replicated 11 several loci that were identified using an earlier generation of this intercross, including an association 12 between locomotor activity and a locus containing a single gene, Csmd1 . We also showed that Csmd1 13 mutant mice recapitulated the locomotor phenotype. Our results demonstrate the utility of this population, 14 identify numerous novel associations, and provide examples of replication in an independent cohort, 15 which is customary in human genetics, and replication by experimental manipulation, which is a unique 16 advantage of model organisms. 17 demonstrated that Csmd1 homozygous mutant mice of heterozygous show a 54% reduction in Csmd1 Residual expression of Csmd1 in homozygous mutant , We analyzed following a saline in 31 and mutant

1 Introduction Genome-wide association studies (GWAS) have revolutionized psychiatric genetics; however, they have also presented numerous challenges.Some of these challenges can be addressed by using model organisms.For example, human GWAS are confounded by environmental variables, such as childhood trauma, which can reduce power to detect genetic associations.In model organisms, environmental variables can be carefully controlled.Furthermore, it has become clear that phenotypic variation in humans is due to numerous common and rare variants of small effect.In model organisms, genetic diversity can be controlled such that all variants are common.In addition, allelic effect sizes in model organisms are dramatically larger than in humans (Flint and Mackay 2009;Parker and Palmer 2011).
Furthermore, because the majority of associated loci are in noncoding regions, expression quantitative trait loci (eQTLs) are useful for elucidating underlying molecular mechanisms (GTEx Consortium et al. 2017;Albert and Kruglyak 2015).However, it remains challenging to obtain large, high quality samples of human tissue, particularly from the brain.In contrast, tissue for gene expression studies can be collected from model organisms under optimal conditions.Finally, the genomes of model organisms can be edited to assess the functional consequences of specific mutations.
Model organism GWAS often employ multigenerational intercrosses because they promote recombination of ancestral haplotypes.We used an advanced intercross line (AIL) of mice, which is the simplest possible multigenerational intercross.AILs, originally proposed by Darvasi and Soller (Darvasi and Soller 1995), are produced by intercrossing two inbred strains beyond the F2 generation.Because the two inbred strains contribute equally to an AIL, all variants are common, and alleles that are identical by state are necessarily identical by descent (IBD), which greatly simplifies phasing and imputation.We performed a GWAS using the world's most advanced mouse AIL, which was created over 50 generations ago by crossing the LG/J (LG) and SM/J (SM) inbred strains (Ehrich et al. 2005a).We investigated over 100 traits using mice from generations 50-56 (G50-56), including locomotor activity, response to methamphetamine, prepulse inhibition (PPI), body weight, and various muscle and bone phenotypes.We also sequenced mRNA from three brain regions and used those data to map eQTLs and identify quantitative trait genes (QTGs) at each locus.Finally, we explored replication of previous associations identified in LG x SM G34 (Cheng et al. 2010;Samocha et al. 2010;Parker et al. 2014;Lionikas et al. 2010;Parker et al. 2011) and used mutant mice to test one of our strongest candidate QTGs.

Results
We used genotyping by sequencing (GBS) to genotype 1,063 of the 1,123 mice that were phenotyped (60 were not successfully genotyped for technical reasons described in the Supplemental Methods).
After quality control, GBS yielded 38,238 autosomal SNPs.In the 24 AIL mice that were also genotyped on the Giga Mouse Universal Genotyping Array (GigaMUGA; Morgan et al. 2015), only 24,934 markers were polymorphic in LG and SM (Supplemental Fig. S1).
LG and SM have been resequenced (Nikolskiy et al. 2015), which allowed us to impute AIL genotypes at ~4.3 million single nucleotide polymorphisms (SNPs; Fig. 1A).Consistent with the expectation for an AIL, the average minor allele frequency (MAF) was high (Fig. 1B).Linkage disequilibrium (LD) decay, which is critical to mapping resolution, has improved since LG x SM G34 (Fig. 1C; Cheng et al. 2010).

LOCO-LMM effectively reduces the type II error rate
Linear mixed models (LMMs) are commonly used to perform GWAS in AILs and other populations that include close relatives (Gonzales and Palmer 2014).SNP data are used to obtain a genetic relationship matrix (GRM); however, this can lead to an inflation of the type II error rate due to proximal contamination (Cheng et al. 2013;Listgarten et al. 2012).We previously proposed using a leaveone-chromosome-out LMM (LOCO-LMM) to address this issue (Cheng et al. 2013).To demonstrate the appropriateness of a LOCO-LMM, we performed a GWAS for albinism, which is a recessive Mendelian trait, using three approaches: a simple linear model, an LMM and a LOCO-LMM (Fig. 2).GWAS using a LOCO-LMM for albinism yielded an association on chromosome 7 (Fig. 2A); accurately identifying the albino locus (Tyr).As expected, p-values from a genome-wide scan using a linear model, which does not account for relatedness, appeared highly inflated (Fig. 2B).This inflation was greatly reduced by fitting a standard LMM, which included SNPs from chromosome 7 in both the fixed and random effects (Fig. 2C).
The LOCO-LMM, which does not include SNPs from the chromosome being tested in the GRM, showed an intermediate level of inflation (Fig. 2D).Was the inflation observed in Fig. 2B-D due to true signal, or uncontrolled population structure?To address this question, we repeated these analyses after excluding SNPs on chromosome 7 from the fixed effect (Fig. 2E-G).Even in the absence of the causal locus, the simple linear model showed substantial inflation, which can only be explained by population structure (Fig. 2E).The standard LMM appeared overly conservative, which we attributed to proximal contamination (Fig. 2F).The LOCO-LMM showed no inflation, consistent with the absence of Tyr and linked SNPs in the fixed effect (Fig. 2G).These results demonstrate the appropriateness of a LOCO-LMM.

Genetic architecture of complex traits in the LG x SM AIL
We used an LD-pruned set of 523,028 autosomal SNPs genotyped in 1,063 mice from LG x SM G50-56 to perform GWAS for 120 behavioral and physiological traits using a LOCO-LMM (Fig. 3A).We used permutation to define a significance threshold of p=8.06x10 -6 (α=0.05).There were 52 loci associated with 33 behavioral traits and 74 loci associated with 17 physiological traits (Fig. 3A, Supplemental Table S1; Supplemental Fig. S2).
To estimate the heritability attributable to SNPs ('SNP heritability'), we calculated the proportion of trait variance explained by the additive effects of 523,028 SNPs.In general, heritability estimates were larger for physiological traits than for behavioral traits (Fig. 3B, Supplemental Table S2), which is consistent with findings in other rodent GWAS (Parker et al. 2016;Nicod et al. 2016;Rat Genome Sequencing and Mapping Consortium et al. 2013).Mean heritability was 0.355 (se=0.045)for physiological traits and 0.168 (se=0.038)for behavioral traits (conditioned place preference, locomotor sensitization, and habituation to startle were not found to have a genetic component and were excluded from the mean).In general, traits with higher heritabilities yielded more associations (Supplemental Fig. S2).However, there was no significant relationship between heritability and effect size at individual loci (Supplemental Fig. S3), suggesting that high heritability does not reliably predict the presence of large-effect alleles.eQTLs For a subset of phenotyped and genotyped mice, we used RNA-sequencing (RNA-seq) to measure gene expression in the hippocampus (HIP), prefrontal cortex (PFC) and striatum (STR) (α=0.05;Fig. 4, Supplemental Fig. S4).We identified 2,902 cis-eQTLs in HIP, 2,125 cis-eQTLs in PFC and 2,054 cis-eQTLs in STR; 1,087 cis-eQTLs were significant in all three tissues (FDR<0.05;Supplemental Table Figure 3. Manhattan plot and heritability for 120 traits measured in the LG x SM AIL.We identified 126 loci for behavioral and physiological traits using 1,063 mice from G50-56 of the LG x SM AIL.A Manhattan plot of GWAS results is shown in (A).Associations for related traits are grouped by color.For clarity, related traits that mapped to the same locus (Supplemental Table S1) are highlighted only once.The dashed line indicates a permutationderived significance threshold of -log10(p)=5.09(p=8.06x10-6 ; α=0.05).(B) For a representative subset of traits, SNP heritability estimates (percent trait variance explained by 523,028 GWAS SNPs) for a subset of traits are shown.Precise estimates of heritability with standard error are provided for all traits in Supplemental Table S2.
73.72 Mb) that was associated with the expression of 85 genes distributed throughout the genome (Fig. 4; Supplemental Table S4).This locus was present in HIP, but not in PFC or STR.

Integration of eQTLs and behavioral GWAS results
Based on results from human GWAS (Albert and Kruglyak 2015;GTEx Consortium et al. 2017), we hypothesized that most loci associated with behavior were due to gene expression differences.For example, four loci associated with locomotor behavior mapped to the same region on chromosome 17 (Supplemental Table S1; Supplemental Fig. S2).The narrowest of these (D1 side changes, 15-20 min; p=3.60x10 -6 ) identified a locus that contains a single gene, Crim1 (cysteine rich transmembrane BMP regulator 1), which had a significant cis-eQTL in HIP.It would be tempting to conclude that Crim1 is the best candidate to explain the associations with locomotor behavior; however, two nearby genes, Qpct (glutaminyl-peptide cyclotransferase) and Vit (vitrin), though physically located outside of the locus, also had cis-eQTLs within the locomotor-associated region (Supplemental Table S3).We therefore consider all three genes valid candidates to explain the association with locomotor behavior.
One of the most significant loci we identified was an association with the startle response, also on chromosome 17 (p=5.28x10-10 ; Fig. 3; Supplemental Fig. S2).This result replicated a previous association with startle from a prior study using G34 mice (Samocha et al. 2010).We performed a phenome-wide association analysis (PheWAS) which showed that this region pleiotropically affected multiple other traits, including locomotor activity following saline and methamphetamine administration (Supplemental Fig. S6).This region was also implicated in conditioned fear and anxiety in our prior studies of G34 mice (Parker et al. 2014), demonstrating that it has extensive pleiotropic effects on behavior.Because the association with startle identifies a relatively large haplotype that includes over 25 eGenes, our data cannot clarify whether the pleiotropic effects are due to one or several genes in this interval.
We also identified a 0.49-Mb locus on chromosome 8 that was associated with locomotor activity (Fig. 5A; Supplemental Table S1); this region was nominally associated with PPI and multiple other activity traits (Fig. 5B).The region identified in the present study (Fig. 5A-C) replicates a finding from our previous study using G34 mice (Cheng et al. 2010; Fig. 5D).In both cases, the SM allele conferred increased activity (Fig. 5C,D) and the implicated locus contained only one gene: Csmd1 (CUB and sushi multiple domains 1; Fig. 5B; Supplemental Table S2); furthermore, the only cis-eQTL that mapped to this region was for Csmd1 (Supplemental Fig. S7).We obtained mice in which the first exon of Csmd1 was deleted to test the hypothesis that Csmd1 is the QTG for this locus.Csmd1 mutant mice exhibited increased activity compared to heterozygous and wild-type mice (Fig. 5E), similar to the SM allele.This result is consistent with the hypothesis that Csmd1 is the causal gene.
We identified seven overlapping loci for locomotor activity on chromosome 4 (Supplemental Table S1; Supplemental Fig. S2).The strongest locus (D5 activity, 0-30 min; p=6.75x10 -9 ) spanned 2.31 Mb and completely encompassed the narrowest locus, which spanned 0.74 Mb (D5 activity, 25-30 min; p=4.66x10 -8 ); therefore, we focused on the smaller region.Oprd1 (opioid receptor delta-1) was a cis-eGene in all three brain regions; the SM allele conferred an increase in locomotor activity and was associated with decreased expression of Oprd1.Oprd1 knockout mice have been reported to display increased activity relative to wild-type mice (Filliol et al. 2000), suggesting that differential expression of Oprd1 could explain the locomotor effect at this locus.However, the presence of other genes and eGenes within the region make it difficult to determine whether Oprd1 is the only QTG.
This locus contained three cis-eGenes and seven trans-eQTLs (Supplemental Fig. S8).One of the trans-eGenes targeted by the locus (Capn5; calpain 5) was most strongly associated with rs108610974 and may be the QTG (Supplemental Table S4).These results illustrate how the knowledge of both cisand trans-eQTLs informed our search for QTGs.

Pleiotropic effects on physiological traits
Because LG and SM were created by selective breeding for large (LG) and small (SM) body size, this AIL is expected to segregate numerous body size alleles (Lawson and Cheverud 2010;Parker et al. 2011).We measured body weight at ten timepoints throughout development and identified 46 associations.There was extensive pleiotropy among body weight, muscle mass, and bone length (Supplemental Table S1; Supplemental Fig. S9-S12).Accounting for pleiotropic genetic architecture, eight major loci arose that influenced body weight at multiple timepoints (Fig. 3A; Supplemental Table S1).
For example, eight body weight timepoints mapped to a region on chromosome 2, where the LG allele was associated with smaller body mass (Supplemental Table S1; Supplemental Fig. S2,S9).The narrowest region spanned 0.08 Mb and did not contain any genes.However, the 0.08-Mb interval contained a cis-eQTL SNP for Nr4a2 (nuclear receptor subfamily 4, group A, member 2) in PFC.Mice lacking Nr4a2 in midbrain dopamine neurons exhibit a 40% reduction in body weight (Kadkhodaei et al. 2009).Consistent with this, the LG allele was associated with decreased expression of Nr4a2.
All body weight timepoints were associated with a locus on chromosome 7 (Supplemental Table S1; Supplemental Fig. S2).We also identified associations for tibialis anterior (TA), gastrocnemius, plantaris weight that partially overlapped this region (Supplemental Fig. S10).Although the most significant SNP associated with muscle weight was ~5 Mb downstream of the top body weight SNP, the LG allele was associated with greater weight at both loci (Supplemental Table S1).For eight out of ten body weight timepoints, the most significant association fell within Tpp1 (tripeptidyl peptidase 1), which was a cis-eGene in all tissues and a trans-eGene targeted by the master HIP eQTL on chromosome 12 (Fig. 4).To our knowledge, Tpp1 has not been shown to affect body size in mice or humans; however, four other cis-eGenes in the region have been associated with human body mass index (Rpl27a, Stk33, Trim66, and Tub) (Berndt et al. 2013;Locke et al. 2015).Dysfunction of Tub (tubby bipartite transcription factor) causes late-onset obesity in mice, perhaps due to Tub's role in insulin signaling (Stretton et al. 2009).In addition, several trans-eGenes map to this interval, including Gnb1 (G protein subunit beta 1), which forms a complex with Tub (Baehr and Frederick 2009).Another trans-eGene associated with this interval, Crebbp (CREB binding protein), has been associated with juvenile obesity in human GWAS (Comuzzie et al. 2012).

Multiple strong associations identified for muscle mass
We examined five hind limb muscle traits, identifying 22 loci (Supplemental Table S1; Supplemental Fig. S2).No loci were identified for soleus weight.The strongest association we identified in this study was for extensor digitorum longus (EDL) weight (p=2.03x10 -1 ; Supplemental Table S1; Supplemental Fig. S11).An association with gastrocnemius weight provided additional support for the region (p=2.56x10 - ; Supplemental Fig. S11); in both cases, the SM allele was associated with increased muscle mass.Each locus spans less than 0.5 Mb and is flanked by regions of low polymorphism between LG and SM (Supplemental Fig. S11, Supplemental Table S1).Two cis-eGenes within the region, Trappc13 (trafficking protein particle complex 13) and Nln (neurolysin), are differentially expressed in LG and SM soleus (Carroll et al. 2017) and TA (Lionikas et al. 2012), with LG exhibiting increased expression of both genes.While there is no known relationship between Trappc13 and muscle, Nln has been shown to play a role in mouse skeletal muscle (Cavalcanti et al. 2014).
The LG allele was associated with greater EDL, plantaris, and TA weight at another locus on chromosome 4 (Supplemental Table S1; Supplemental Fig. S12).The loci for EDL and plantaris spanned ~0.5 Mb, defining a region that contained six genes (Supplemental Table S1).The top SNPs for EDL (rs239008301; p=7.88x10 -13 ) and plantaris (rs246489756; p=2.25x10 -6 ) were located in an intron of Astn2 (astrotactin 2), which is differentially expressed in LG and SM soleus (Carroll et al. 2017).SM, which exhibits lower expression of Astn2 in soleus relative to LG (Carroll et al. 2017), has a 16 bp insertion in an enhancer region 6.6 kb upstream of Astn2 (ENSMUSR00000192783; Nikolskiy et al. 2015).Two other genes in this region have been associated with muscle or bone phenotypes traits in the mouse: Tlr4 (toll-like receptor 4), which harbors one synonymous coding mutation on the SM background (rs13489095) and Trim32 (tripartite motif-containing 32), which contains no coding polymorphisms between the LG and SM strains.

Discussion
Crosses among well-characterized inbred strains are a mainstay of model organism genetics.In the present study, we used 1,063 male and female mice from LG x SM G50-56 to identify 126 loci for a variety of traits selected for their relevance to human psychiatric and metabolic diseases (Lawson and Cheverud 2010;Swerdlow et al. 2016; de Wit and Phillips 2012) (Fig. 3; Supplemental Table S1; Supplemental Fig. S2).Whereas our previous work established AILs as an effective tool for finemapping loci identified in F2 crosses (Gonzales and Palmer 2014;Carroll et al. 2017;Parker et al. 2014;Cheng et al. 2010;Samocha et al. 2010;Parker et al. 2011;Lionikas et al. 2010), this study demonstrates that AILs are also a powerful fine mapping population in their own right.We show that several QTGs we identified are corroborated by extant human and mouse genetic data.We also replicated a number of our earlier findings.
Classical crosses like F2 and recombinant inbred strains provide poor mapping resolution because the ancestral chromosomes persist as extremely long haplotypes (Parker and Palmer 2011).To address this limitation, we and others have used AILs (Gonzales and Palmer 2014).Because both inbred strains contribute equally to an AIL, there are numerous common variants (Fig. 1A,B), and each successive generation further degrades LD between adjacent SNPs (Fig. 1C).In addition to AILs, a number of other outbred populations have been used in rodent GWAS, including CFW mice (Nicod et al. 2016;Parker et al. 2016), Diversity Outbred (DO) mice (Gatti et al. 2014;Logan et al. 2013), and N/NIH heterogeneous stock rats (HS) (Rat Genome Sequencing and Mapping Consortium et al. 2013;Tsaih et al. 2014).CFW mice are obtained from a commercial vendor, which avoids the expenses of maintaining a colony.In addition, non-siblings can be obtained, which reduces the complicating effects that can occur when close relatives are used in GWAS.However, the CFW founder strains are unavailable, and many alleles exist at low frequencies among CFWs, limiting power and introducing genetic noise (Nicod et al. 2016;Parker et al. 2016).Commercially available DO mice are more expensive than CFWs, but like AILs, the founder strains have been fully sequenced, which allows imputation of SNPs and founder haplotypes.However, three of the eight inbred strains used to produce the DO are so-called wild-derived strains; making DO mice very difficult to handle, which complicates many behavioral procedures (Logan et al. 2013).Furthermore, the causal alleles in the DO are often from one of the wild derived strains, because 8 strains contributed equally to the DO, this means that the causal allele frequencies are often in the range of 0.125.Finally, N/NIH HS rats, which are conceptually very similar to DO mice, have also been used as a fine mapping population (Rat Genome Sequencing and Mapping Consortium et al. 2013;Tsaih et al. 2014).Among these options, AILs stand out for their simplicity, balanced allele frequency and ease of handling.
A major goal of this study was to identify the genes that are responsible for the loci implicated in behavioral and physiological traits.The mapping resolution of the LG x SM AIL was critically important for this goal.However, no matter how precise the resolution, proximity of a gene to the associated SNP is never sufficient to establish causality (Albert and Kruglyak 2015).Therefore, we used RNA-seq to quantify mRNA abundance in three brain regions that are strongly implicated in the behavioral traits: HIP, PFC and STR.We used these data to identify 7,081 cis-eQTLs and 1,372 trans-eQTLs (Fig. 4, Supplemental Fig. S4,S5; Supplemental Tables S3,S4).In a few cases, loci contained only a single eQTL; however, in most cases multiple cis-eQTLs and trans-eQTLs mapped to the implicated loci.Thus, we frequently incorporated functional information, including data about tissue specific expression, coding SNP, mutant mice, and human genetic studies to parse among the implicated genes.
We have previously shown that GBS is a cost-effective strategy for genotyping CFW mice (Fitzpatrick et al. 2013).The advantages of GBS were even greater for this AIL because imputation allowed us to easily obtain 4.3M SNPs while using only half the sequencing depth (Fig. 1A).Even before imputation, GBS yielded nearly 50% more informative SNPs compared to the best available SNP genotyping chip (Morgan et al. 2015) at about half the cost (Supplemental Fig. S1).
One of the goals of this study was to perform GWAS for conditioned place preference (CPP), which is a well-validated measure of the reinforcing effects of drugs (Tzschentke 2007).Unfortunately, the heritability of CPP in this study was not significantly different from zero (Fig. 3B).This result was partially consistent with our prior study in which we used a higher dose of methamphetamine (2 vs. the 1 mg/kg used in the present study) (Bryant et al. 2012).The low heritability of CPP in this AIL likely reflects a lack of relevant genetic variation in this specific population since both panels of inbred strains and genetically engineered mutant alleles show differences in CPP (Tzschentke 2007;Philip et al. 2010;Martinelli et al. 2016), demonstrating the existence of heritable variance in other populations.It is possible that even lower doses of methamphetamine, which might fall on the ascending portion of the dose-response function, would have resulted in higher heritability.Similarly, responses to other drugs or different CPP methodology may have exhibited greater heritability.
We used PheWAS to identify pleiotropic effects of several loci identified in this study.In many cases, pleiotropy involved highly correlated traits such as body weight on different days or behavior at different time points within a single day (Supplemental Fig. S9-S13; Supplemental Table S1).We also observed more surprising examples of pleiotropy, including pleiotropy between locomotor activity and gastrocnemius mass on chromosome 4 (Supplemental Fig. S14), and pleiotropy between locomotor activity and the startle response on chromosome 12 (Supplemental Fig. S15).We also observed extensive pleiotropy on chromosome 17; this locus influenced saline-and methamphetamine-induced locomotor activity and startle response (Supplemental Fig. S6).Moreover, this same region had been previously implicated in anxiety-like behavior (Parker et al. 2014), contextual and conditioned fear (Parker et al. 2014), and startle response (Samocha et al. 2010) in prior studies of LG x SM G34, suggesting that the locus has a broad impact on many behavioral traits.These results support the idea that pleiotropy is a pervasive feature in this AIL and provide further evidence of the replicability of the loci identified by this and prior GWAS.
Discoveries from human GWAS are often considered preliminary until they are replicated in an independent cohort.In model organisms, replication using an independent cohort is rarely employed because it is possible to directly manipulate the implicated gene.We replicated one behavioral locus identified in this study using the criteria of both human and model organism genetics.We had previously identified an association with locomotor activity on chromosome 8 using G34 of this AIL (Cheng et al. 2010), which was replicated in the present study (Fig. 5).In both cases, the SM allele was associated with lower activity (Fig. 5C-D).We also identified a locus for PPI (76 dB) in this region (Fig. 5A; Supplemental Table S1, Supplemental Fig. S2).The loci identified in both G34 and in G50-56 were small and contained just one gene: Csmd1 (Fig. 5B).In the present study we also identified a cis-eQTL for Csmd1 in HIP (Supplemental Fig. S7).Finally, we obtained Csmd1 mutant mice (Distler et al. 2012) and found that they also showed altered locomotor activity (Fig. 5E).Thus, we have demonstrated replication both by performing an independent GWAS and by performing an experimental manipulation that recapitulates the phenotype.

Genetic background
The LG and SM inbred strains were independently selected for high and low body weight at 60 days (Beck et al. 2000).The LG x SM AIL was derived from an F1 intercross of SM females and LG males initiated by Dr. James Cheverud at Washington University in St. Louis (Ehrich et al. 2005a).Subsequent AIL generations were maintained using at least 65 breeder pairs selected by pseudo-random mating (Ehrich et al. 2005b).In 2006, we established an independent AIL colony using 140 G33 mice obtained from Dr. Cheverud (Jmc: LG,SM-G33).Since 2009, we have selected breeders using an R script (Supplemental Methods) that leverages pairwise kinship coefficients estimated from the AIL pedigree to select the most unrelated pairs while also attempting to minimize mean kinship among individuals in the incipient generation (the full pedigree is included in Supplemental File S1).We maintained ~100 breeder pairs in G49-55 to produce the mice for this study.In each generation, we used one male and one female from each nuclear family for phenotyping, and reserved up to three of their siblings for breeding the next generation.
PPI and startle: PPI is the reduction of the acoustic startle response when a loud noise is immediately preceded by a low decibel (dB) prepulse (Graham 1975).PPI and startle are measured across multiple trials that occur over four consecutive blocks of time (Samocha et al. 2010).The primary startle trait is the mean startle amplitude across all pulse-alone trials in blocks 1-4.Habituation to startle is the difference between the mean startle response at the start of the test (block 1) and the end of the test (block 4).PPI, which we measured at three prepulse intensities (3, 6, and 12 dB above 70 dB background noise), is the mean startle response during pulse-alone trials in blocks 2-3 normalized by the mean startle response during prepulse trials in blocks 2-3.

Physiological traits:
We measured body weight (g) on each testing day and at the time of death.
One week after PPI, we measured blood glucose levels (mg/dL) after a four-hour fast.One week after glucose testing, we killed the mice, and measured tail length (cm from base to tip of the tail).We stored spleens in a 1.5 mL solution of 0.9% saline at -80C until DNA extraction.We removed the left hind limb of each mouse just below the pelvis; hind limbs were stored at -80C.Frozen hind limbs were phenotyped by Dr. Arimantas Lionikas at the University of Aberdeen.Phenotyped muscles include two dorsiflexors, TA and EDL, and three plantar flexors: gastrocnemius, plantaris and soleus.We isolated individual muscles under a dissection microscope and weighed them to 0.1 mg precision on a Pioneer balance (Ohaus, Parsippany, NJ, USA).After removing soft tissue from the length of tibia, we measured its length to 0.01 mm precision with a Z22855 digital caliper (OWIM GmbH & Co., Neckarsulm, GER).

Brain tissue:
We collected HIP, PFC and STR for RNA-seq from the brain of one mouse per cage.This allowed us to dissect each brain within five min of removing a cage from the colony room (rapid tissue collection was intended to limit stress-induced changes in gene expression).We preselected brain donors to prevent biased sampling of docile (easily caught) mice and to avoid sampling full siblings, which would reduce our power to detect eQTLs.Intact brains were extracted and submerged in chilled RNALater (Ambion, Carlsbad, CA, USA) for one min before dissection.Individual tissues were stored separately in chilled 0.5-mL tubes of RNALater.All brain tissue was dissected by the same experimenter and subsequently stored at -80°C until extraction.

GBS variant calling and imputation
GBS is a reduced-representation genotyping method (Elshire et al. 2011;Grabowski et al. 2014) that we have adapted for use in mice and rats (Parker et al. 2016;Fitzpatrick et al. 2013).We extracted DNA from spleen using a standard salting-out protocol and prepared GBS libraries by digesting DNA with the restriction enzyme PstI, as described previously (Parker et al. 2016).We sequenced 24 uniquely barcoded samples per lane of an Illumina HiSeq 2500 using single-end, 100 bp reads.We aligned 1,110 GBS libraries to the mm10 reference genome before using GATK (DePristo et al. 2011) to realign reads around known indels in LG and SM (Nikolskiy et al. 2015) (see Supplemental Methods and Supplemental File S3 for details and example commands).We obtained an average of 3.2M reads per sample.We discarded 32 samples with <1M reads aligned to the main chromosome contigs (1-19, X, Y) or with a primary alignment rate <77% (i.e. three s.d.below the mean of 97.4%; Supplemental Fig. S17).
We used ANGSD (Korneliussen et al. 2014) to obtain genotype likelihoods for the remaining 1,078 mice and used Beagle (Browning and Browning 2007) for variant calling, which we performed in two stages.We used first-pass variant calls as input for IBDLD (Han and Abney 2011;Abney 2008), which we used to estimate kinship coefficients for the mice in our sample.Because our sample contained opposite-sex siblings, we were able to identify and resolve sample mix-ups by comparing genetic kinship estimates to kinship estimated from the LG x SM pedigree (described in the Supplemental Methods).In addition, we re-genotyped 24 mice on the GigaMUGA (Morgan et al. 2015) to evaluate GBS variant calls (Supplemental Table S5 lists concordance rates at various stages of our pipeline; see Supplemental Methods for details).
After identifying and correcting sample mix-ups, we discarded 15 samples whose identities could not be resolved (Supplemental Methods).Next, we used Beagle (Browning andBrowning 2016, 2007), in conjunction with LG and SM haplotypes obtained from whole-genome sequencing data (Nikolskiy et al. 2015) to impute 4.3M additional SNPs into the final sample of 1,063 mice.We removed SNPs with low MAFs (<0.1),SNPs with Hardy-Weinberg Equilibrium (HWE) violations (p ≤ 7.62x10 -6 , determined from gene dropping simulations as described in the Supplemental Methods), and SNPs with low imputation a second genome-wide scan for each trait, this time dropping the fixed effect of dosage and including the complete GRM estimated from SNPs on all 19 autosomes.We repeated the procedure using dosage at the most significant SNP as a covariate for each trait and interpreted the difference between the two estimates as the effect size of that locus.

RNA-sequencing and quality control
We extracted mRNA from HIP, PFC and STR as described in Parker et al. (2016) and prepared cDNA libraries from 741 samples with RNA integrity scores ≥ 8.0 (265 HIP; 240 PFC; 236 STR) (Schroeder et al. 2006) as measured on a Bioanalyzer (Agilent, Wilmington, DE, USA).We used Quant-iT kits to quantify RNA (Ribogreen) and cDNA (Picogreen; Fisher Scientific, Pittsburgh, PA, USA).
Barcoded sequencing libraries were prepared with the TruSeq RNA Kit (Illumina, San Diego, USA), pooled in sets of 24, and sequenced on two lanes of an Illumina HiSeq 2500 using 100 bp, single-end reads.
Because mapping quality tends to be higher for reads that closely match the reference genome (Degner et al. 2009), read mapping in an AIL may be biased toward the reference strain (C57BL/6J) (Wang and Clark 2014).We addressed this concern by aligning RNA-seq reads to custom genomes created from LG and SM using whole-genome sequence data (Nikolskiy et al. 2015).We used default parameters in HISAT (Kim et al. 2015) for alignment and GenomicAlignments (Lawrence et al. 2013) for assembly, assigning each read to a gene as defined by Ensembl (Mus_musculus.GRCm38.85)(Aken et al. 2016).We required that each read overlap one unique disjoint region of the gene.If a read contained a region overlapping multiple genes, genes were split into disjoint intervals, and any shared regions between them were hidden.If the read overlapped one of the remaining intervals, it was assigned to the gene that the interval originated from; otherwise, it was discarded.Next, we reassigned the mapping position and CIGAR strings for each read to match mm10 genome coordinates and combined the LG and SM alignment files for each sample by choosing the best mapping.Only uniquely mapped reads were included in the final alignment files.We then used DESeq (Anders and Huber 2010) to obtain normalized read counts for each gene in HIP, PFC and STR.We excluded genes detected in <95% of samples within each tissue.

Figure 1 .
Figure 1.SNPs, minor allele frequencies (MAFs) and linkage disequilibrium (LD) decay in the LG x SM AIL.Imputation provided ~4.3 million SNPs.Filtering for LD (r 2 ≥ 0.95), MAF < 0.1, and HWE (p ≤ 7.62x10 -6 ) resulted in 523,028 SNPs for GWAS.(A) SNP distribution and density of GWAS SNPs are plotted in 500 kb windows for each chromosome.As shown in Supplemental Fig.1, regions with low SNP density correspond to regions predicted to be nearly IBD in LG and SM (Nikolskiy et al. 2015).(B) MAF distributions are shown for ~4.3 million imputed SNPs (gold; unfiltered) and for the 523,028 SNPs used for GWAS (orange; filtered).Mean MAF is the same in both SNP sets.(C) Comparison of LD decay in G50-56 (dark purple) and G34 (light purple) of the LG x SM AIL.Each curve was plotted using the 95 th percentile of r 2 values for SNPs spaced up to 5 Mb apart.

Figure 2 .
Figure 2. GWAS for albinism verifies that the LOCO-LMM effectively controls type I and type II error.We conducted a GWAS for albinism, a Mendelian trait caused by the Tyr locus on mouse chromosome 7, using three models: a linear model, an LMM, and a LOCO-LMM.We also repeated each scan after excluding SNPs on chromosome 7.A Manhattan plot of results from the LOCO-LMM is shown in (A).Quantile-quantile plots of expected vs. observed p-values are shown for (B) a simple linear model that does not account for relatedness; (C) a standard LMM that includes all GWAS SNPs in the genetic relatedness matrix (GRM; i.e. the random effect); and (D) a LOCO-LMM whose GRM excludes SNPs located on the chromosome being tested.Plots (E-G) show results after excluding chromosome 7 from the GWAS.

Figure 5 .
Figure 5. Replication of an association between Csmd1 and locomotor activity.(A) Regional plot drawn from all 4.3M SNPs showing the association between rs33436747 and D5 activity levels.The location of Csmd1, 1.5-LOD interval (gold bar), areas of elevated recombination from Brunschwig et al. (green plus symbols), regions predicted to be nearly IBD between LG and SM by Nikolskiy et al. (grey bars), and SNP MAFs (grey heatmap) are indicated.Points are colored by LD (r 2 ) with rs33436747.The dashed line indicates a significance threshold of -log10(p)=5.09(α=0.05).(B) PheWAS plot of associations between rs33436747 and other behavioral traits measured in G50-56 mice.(C) Bar plot of quantile-normalized residuals of locomotor activity at the Csmd1 locus are plotted for G50-56 mice.(D) Bar plot of quantile-normalized residuals of locomotor activity at the Csmd1 locus for G34 mice from Cheng et al. rs33436747 was not genotyped in G34; therefore, we plotted activity by genotype at the nearest SNP (rs33014260; 6,764 bp upstream of rs33436747).(E) Bar plot of locomotor activity data (distance traveled in 0-30 min) for Csmd1 mutant mice.In panels c-e the number of mice in each genotype class is shown below the corresponding bar.ANOVA and two-sided t-test (95% confidence level) p-values are shown for each comparison.