The Effects of Mutations are Modified by Genetic Background in Mice

Background: Commonly used genetic models of phenotype assume that alleles have similar effects in all individuals. Mounting evidence suggests that higher-order interactions, which may vary between populations, profoundly affect phenotypic variation. The effect of genetic and environmental perturbations may therefore differ across genetic backgrounds. In this manuscript, we develop a statistical test that determines whether the effect of a mutation on a complex phenotype changes as a function of ancestral background. Results: We apply our test to two large cohorts of mice and observe 49 significant gene by ancestry interaction associations across 14 phenotypes as well as over 1,400 Bonferroni-corrected gene by ancestry interaction associations in gene expression data. We also observe evidence of rapid selection pressure on individual polymorphisms within one of the cohorts. Conclusions: Unlike our prior work in human populations, we observe widespread evidence of ancestryspecific SNP effects, perhaps reflecting the greater divergence present in crosses among laboratory mice.


Background
Genetic association studies in humans and model organisms have identified a number of significant links between individual polymorphisms and phenotypic variability. A fundamental assumption of many of these studies is that an allele will have a similar effect in each member of the population, that is, that epistatic and other higher-order interactions across the genome can largely be ignored [1][2][3][4] . While we have previously observed only modest evidence of ancestry specific genetic effects in humans 5 , model organisms are often further diverged than human populations. For example, we observed radically different phenotypic consequences of null alleles of Tcf7l2 and Cacna1c when expressed on different inbred strain backgrounds 6 . We therefore developed a test, designed for use in model organisms, to identify and examine heterogeneous SNP effects across different genetic backgrounds.
The field of mouse genetics provides numerous opportunities to identify and study the genetic basis of phenotypic variation in mice with obvious extensions to human populations. Many prior studies of epistasis in both mice and humans have used pairwise interactions between loci to attempt to identify interacting loci 4,[7][8][9] , yet this approach is hampered by a substantial increase in the threshold of significance due to the increased numbers of tests (from n to n^2, where n are the number of SNPs tested), which in turn requires larger cohorts to overcome. In this study, we approached epistasis in a different way by examining two distinct mouse cohorts and considering the effects of interactions between individual polymorphisms and global ancestry (Ɵ), which we defined as the percentage of the genome inherited from one of the parental lines of each cohort. We utilized a model we recently developed 5 which explored these interactions in human populations and extended this idea to model organisms by incorporating a linear mixed model capable of accounting for the highly structured relatedness of recombinant inbred (RI) lines 10 . Although our initial model was designed to identify both genetic background and environmental effects, the controlled environment in which mouse cohorts live is intended to minimize sources of environmental variance. This test, which we have called GxƟ, acts as a surrogate measure of the concurrent action of many other SNPs, and allowed us to ask the question of whether or not SNP effect sizes change as a function of overall genetic ancestry. Our observed significant associations validate this hypothesis. In contrast to a pairwise epistatic test where significance means a detected interaction between two specific loci, a natural interpretation of a significant GxƟ interaction is widespread epistasis with the tested genotype G with many loci across the genome, as we showed previously 5 .
We first evaluated the statistical properties of the method on simulated phenotypes from real recombinant inbred line data sets. The GxƟ test successfully distinguishes between effects dependent only on SNP or ancestry and effects which arise from their interaction. We then applied the test to two different populations of mice: recombinant inbred (RI) lines that are a subset of the Hybrid Mouse Diversity Panel (HMDP) 11 , and a 50 th generation intercross between the inbred LG/J and SM/J mouse strains 12 , whose ancestry proportions can be clearly and precisely determined 12 Next, we conducted a replication study demonstrating that our approach is able to replicate findings across similar populations.
Finally, we examined the ancestry distribution at specific sites in the RI lines and observed that for many positions there is a statistically significant depletion or increase of individual ancestries. We interpret this as evidence of selection during the process of RI strain derivation and identify regions of the genome with strong selection up to 4 Mb. These regions are evidence of strong GxƟ interactions which inhibit or promote the transmission of these specific alleles to subsequent generations. These regions are enriched for genes involved in cancer and organogenesis and are enriched (P=1.8E-4) for metallopeptidases, which play key roles in fertility and neo/perinatal lethality [13][14][15] . We also observed several cases in which differential effect sizes were observed for SNP-ancestry interactions. This

Method Overview:
Our objective is to determine whether a SNP has a different effect size a function genetic background in a mouse cross. As described in our previous work<cite>, this is effectively achieved through an interaction test between each SNP and global genetic ancestry Ɵ. Due to differences in relatedness structures of individuals in human populations and mouse crosses, we extended our approach with a linear mixed model to accommodate more complex relatedness structures (see Methods).

Simulated Data:
To examine the properties of our approach in a mouse cross, we applied the method to phenotypes simulated using real genotypes that reflect the underlying relatedness present within the recombinant inbred (RI) lines of the Hybrid Mouse Diversity Panel (HMDP). We Changes to either βG or βGxƟ did not have an effect on the likelihood of identifying an association in the other term, indicating that these two terms are estimated correctly independently of one another. We also explored incorporating a second genetic relationship matrix (GRM) accounting for population structure arising from descent from a single ancestor (the B6 mouse). 17 However, we observed no significant improvement in power when incorporating this second GRM ( Figure 1). Consequently, we used a 1 GRM model in our analyses of real phenotypic data.

Mouse Populations
Next, we applied the method to two large panels of mice to identify GxƟ effects, where a given mutation interacts epistatically with one or many other loci, captured in the model by Ɵ, the global ancestry. Our first cohort is the HMDP, a set of 150+ commercially available inbred strains 18 . Numerous GWAS have been performed in the HMDP, including several using PYLMM, which forms the core of our algorithm [19][20][21] . The largest component of the HMDP is comprised of 122 RI strains (28 AxB, 71 BxD, 12 BxH, 11 CxB). Each RI was constructed from re-derivation of novel inbred lines via brother-sister mating following an F2 cross between the sub-panel's parental lines. In the case of the HMDP, one of the parental strains for each RI strain was the commonly studied strainC57BL/6J (B6). Using B6 as the 'ancestral line', we calculated thetas for each RI strain (Fig 2a). As expected, the average ancestry attributable to B6 was roughly 50% (50.63%) and roughly normally distributed. We removed a single outlier, BXD32/TyJ, whose B6 Ancestry of 25.41% reflects a previously known additional backcross to DBA/2J, resulting in a strain that is 75% DBA and 25% B6. Each study using the HMDP uses a different subset of the entire panel, and we selected a study on heart failure induced by the chronic beta adrenergic agonist isoproterenol 22 for analysis as it used the most RI strains compared to other published HMDP data. We used 123 clinical phenotypes in conjunction with microarray-derived gene expressions measured in the left ventricle and on average tested 67 RI strains per phenotype.
The second cohort consists of 1,063 animals from the F50 -F56 generation of an advanced intercross line (AIL) created by crossing the LG/J and SM/J inbred mouse strains 12 . Unlike the RI strains from the HMDP, AILs are maintained in a manner that minimizes inbreeding. We arbitrarily set LG/J as the 'ancestral' strain of interest and calculated thetas for each of the 1,063 mice in the panel (Fig 2b).
For this study, we focused on a diverse set of 133 phenotypes that had been measured in these mice.
We describe the results of the GxƟ associations in each panel before demonstrating replication of signal in a phenotypic trait as well as in expression data.

Evidence of SNP x Ancestry Interactions in Phenotypic Data
We first applied the GxƟ method to the 123 observed heart-failure related phenotypes from the RI strains. We observed well-calibrated statistics, with λGC equal to 0.978 and λGxƟC equal to 1.045. We observed 44 significant GxƟ loci across 9 phenotypes: E/A Ratio, Free Fatty Acid content in the blood, Cardiac Fibrosis, Fractional Shortening of the heart during contraction, Heart Rate, Internal Diameter of the Left Ventricle, left ventricular mass and left and right atrial weights (Table S1). These GxƟ loci were largely distinct from previously reported GWAS loci in the same phenotypes in this panel of mice 19,23 , yet contained a number of highly relevant genes, as discussed below.
By way of example, we focus on two important phenotypes from the HMDP panel. Cardiac fibrosis is a marker of cardiac dysfunction. Genes identified through the GxƟ screen (Fig 3a) as potential candidates include: Crisp2 (rs6295287 , p=6.74E-7), a secreted biomarker of cardiovascular disease 24 , Top2b (rs31538570, p=2.14E-06) which plays a cardioprotective role in response to stress 25 , Rarb (rs31538570, p=2.14E-06) a known regulator of inflammation with unknown function in the heart 26 and Fibrosin (rs33146511, p=2.54E-6) a major component of the fibrosis pathway 27 . Left ventricular mass increase in response to catecholamine challenge is the primary marker of cardiac hypertrophy in the HMDP. A single GxƟ locus(rs31313229, p=2.91E-06) on chromosome 3 (Fig 3b) contains the gene Lphn2, which has a role in the promotion of cellular adhesion in response to external stimuli. 28 We next applied the GxƟ method to 133 phenotypes measured in the LG/J x SM/J AIL. As with the HMDP data, we observed well-calibrated statistics in this cross, with λGC equal to 0.995 and λGxƟC equal to 1.033. Despite the larger sample size, we observed only 5 significant loci across 4 phenotypes (activity levels in a saline-injected animal on day 3 of a conditioned place preference test; average weight of animals across 5 different time points roughly a week apart; glucose (mg/dL) in blood after a 4 hour fast [Fig 3c]; weight at ~68 days of age) (Table S1). This smaller number of GxƟ interactions could be due to the differences in strain as well as differences in phenotypes.

GxƟ Associations in HMDP Cohort Gene Expression
We next examined gene expression in the hearts of the HMDP cohort, where we had transcriptome microarrays from both a control and treated condition. Each cohort consists of approximately 70 RI lines, with 66 lines overlapping between the two cohorts (full lists of strains in Table   S2). We examined 13,155 expressed and varying (CV > 5%) genes from the left ventricles of the HMDP cohorts using the ~170k SNPs (MAF <= 0.05). We observed 1,486 significant associations with 18 genes at a Bonferroni-corrected P value of 3.2E-10 (135,130 associations with 1,350 genes at GW-significant threshold of 4.2E-6) in the control cohort and 597 significant associations with 39 genes at the same threshold in the treated cohort (32,043 associations with 1,042 genes at 4.2E-6). The reduction of significant associations in the treated cohort is likely due to reduced power due to increased phenotypic

Replication of GxƟ Associations across eQTLs from HMDP Cohort
Although the catecholamine treatment used to induce heart failure in this study does affect gene expression, the majority (12178, 92.6%) of genes were not significantly (Student's T-test) affected by the drug. To demonstrate that the method is able to replicate GxƟ results across cohorts, we examined the reproducibility of expression GxƟ QTLs in the treated and untreated RI lines of the HMDP. Of the 1,486 associations observed in the control data, we observe 305 (21%) with a GxƟ signal (P<.05) in the treated cohort (36 at FDR< 5%), suggesting strong replication between the two cohorts despite differences in genetic background and the reduced power in the treated cohort caused by the effects of the catecholamine drug.

Fitness
Strong GxƟ effects have the potential to affect overall organismal fitness. For example, strong GxƟ effects may lead to a differential frequency of strain-specific genotypes at individual loci, as those loci interact with the rest of the genome to cause changes to fitness and, consequently, retention of those sites. Alternately, it is possible that these variations may occur due to absolute differences in fitness of that SNP without GxƟ effects: however, both of these results are interesting and evidence of selection after admixture.
We conducted a test to search for individual loci with enriched or depleted B6 ancestry in the RI strains. As expected, the average B6 ancestry across all SNPs was 50.68% +/-8.1%, which is indistinguishable from the ancestry by strain average of 50.63% +/-6.7% (Fig 4a). At the level of individual loci, however, we observed significant variation in B6 ancestry across the genome (Fig 4b), with some SNPs displaying very low or very high frequency (Fig 4c). We calculated the statistical likelihood of detecting the observed ancestries at all loci across the genome (λ=1.07, Fig 4d) and identified 614 SNPs from the 170k original SNPs with significantly altered ancestries at a FDR of 5% (Table S3). Nine loci were identified in which 5 or more SNPs were located together (Table 1)

Discussion
We present a test which we call GxƟ, that leverages admixed populations such as inbred mouse strains to identify epistatic interactions between a SNP and an unknown number of other loci summarized into a single genomic ancestry score, Ɵ (available at https://github.com/ChristophRau/GxTheta). One major advantage of this approach is that, unlike other epistasis-focused association approaches, it does not increase the number of tests when compared to a typical GWAS, resulting in similar genome-wide significance thresholds. Our method sought to circumvent the traditional pitfall in Epistasis GWAS by examining each SNP only once for an interaction with a global genomic ancestry. We examined two populations of mice, the recombinant inbred lines of the Hybrid Mouse Diversity Panel 32 and an AIL based on LG/J and SM/J 12 . Despite nearly ten times the number of genetically distinct mice, as well as a larger number of phenotypes and genotypes in the AIL, we observed approximately nine times more significant GxƟ peaks in the HMDP phenotypes (44 vs 5). Several possible reasons exist for this difference. First, outbred mice have significantly lower power due to increased phenotypic variance caused by high rates of heterozygosity at alleles while inbred populations have increased power due to a lack of this variance.
This can also be seen in our regular GWAS results (see Supp) 12,22,40 . Second, the different phenotypes studied in each cohort will have different genetic architecture and experimental noise (e.g. the AIL study includes many behavioral traits while the HMDP does not). Third, the genetic backgrounds of the two cohorts are different, which might also contribute to differences in observed numbers of GxƟ interactions. Although individually, LG and SM are only slightly less genetically diverged from one another as any pair of strains that make up the RI panels of the HMDP (Table S4), the presence of five ancestral lines in the HMDP compared to only two in the F50 cross results in much more genetic diversity in the HMDP compared to the F50 mice. Finally, we observe a higher variance in ancestral background in the HMDP when compared to the AIL. As our method relies on differences in ancestral background to identify sites with different effect sizes in different genetic contexts, the reduced variance in the AIL lines necessarily corresponds to a decrease in the power to detect GxƟ interactions.
Taken together, our method is best suited to datasets with relatively low heterozygosity, clear and numerous differences in genetic background, and with higher variance in the percentage of SNPs attributable to a given ancestry.

Conclusion
The results of our study suggest that heterogenous SNP effects due to differing ancestries is pervasive in mouse populations, especially in diverse populations such as the HMDP. This observation matches prior observation of significant epistatic interactions in inbred strains of mice as well as examples in human studies 5 . Further analyses of these heterogenous-effect SNPs may reveal novel epistatic interactions which drive phenotypic expression, and suggests that careful attention to genetic ancestry should be considered when studying the role of an individual polymorphism on a phenotype.

Mouse Populations and Ancestry Calling
Mouse data were drawn from previously reported studies 12,41 . The LGxSM AIL consists of 1,063 G50-56 mice derived from an original F1 intercross between the LG and SM inbred lines. The Hybrid Mouse Diversity Panel consists of over 150 strains of commercially available inbred mice 32 , of which 122 strains were recombinant inbred lines and suitable for our study. SNPwise ancestries were determined by identifying all SNPs which differed between parental lines (AIL: LG and SM or HMDP:C57BL/6 and A, C3H, DBA/2 or BALB/c). Genotypes from the G50-56 or RI lines were filtered for these SNPs and ancestries calculated using either LG or C57BL/6 as the strain background of interest.

Simulation Framework
We created sets of simulated phenotypes based on the genotypes of the HDMP RI panel, which is an admixed population in which the B6 strain, on average, contributes 50% of each strain's DNA. For each simulated phenotypes, we drew a SNP (MAF > 5%) at random from the HMDP genotypes and created a phenotype based on β, the genetic effect size, , ϕ the effect size of the interaction between global ancestry (Ɵ) and the chosen SNP and three variance terms: 2 , the proportion of variance attributable to genetic effects Ɵ 2 , the proportion of variance attributable to GxƟ effects and 2 , the residual proportional variance attributable to all combined sources of error and variation not considered in this study. Phenotypes were generated both with and without the GxƟ variance term to ascertain the necessity of incorporating a second GRM (K A ) into the algorithm.

GxƟ Model
The equation to determine the effects of a SNP and a SNP x Ancestry term on a phenotype can be written as: Where is the phenotype of person k, μ is the mean phenotypic value, M is the number of markers, β are the weights on the SNPs, X is the mxn array of SNP genotypes, δ is the global weight of the ancestry effect, Ɵ is the ancestries for all N individuals. ϕ are the weights of the GxƟ effect and ε is the combined error term. We want to identify SNPs where ≠ 0 as these are sites where Ancestry is interacting with our genotypes.
We can rewrite the above as: and R01HG006399)). The funding bodies played no role in the study other than funding.

Author Contributions:
CDR designed the project, created the GxƟ software, acquired the HMDP data and analyzed and interpreted it and the AIL data and drafted the manuscript. NMG assisted in the design of the AIL cohort, acquired the AIL data and substantially edited the manuscript. DP assisted in conceptualizing the GxƟ algorithm and substantially revised the manuscript. AAP assisted in the design of the AIL cohort and substantially revised the manuscript. AJL assisted in the design of the HMDP cohort and substantially revised the manuscript. NZ conceptualized and designed the GxƟ algorithm, interpreted the data and helped to draft and substantially revise the manuscript. All authors read and approved the final manuscript. Blue and Purple power curves are the power curves for detecting a significant SNP effect, while the orange and green curves are the power curves for detecting a significant GxƟ effect. Two phenotypic models, one incorporating 1 GRM (1K) correcting for relatedness in the SNPs (green, purple) and one incorporating 2 GRMs (2K) correcting for relatedness in both SNPs and Ancestry (red, blue) were used.