Weak Epistasis Generally Stabilizes Phenotypes in a Mouse Intercross

The extent and strength of epistasis is commonly unresolved in genetic studies, and observed epistasis is often difficult to interpret in terms of biological consequences or overall genetic architecture. We investigated the prevalence and consequences of epistasis by analyzing four body composition phenotypes—body weight, body fat percentage, femoral density, and femoral circumference—in a large F2 intercross of B6-lit/lit and C3.B6-lit/lit mice. We used Combined Analysis of Pleiotropy and Epistasis (CAPE) to examine interactions for the four phenotypes simultaneously, which revealed an extensive directed network of genetic loci interacting with each other, circulating IGF1, and sex to influence these phenotypes. The majority of epistatic interactions had small effects relative to additive effects of individual loci, and tended to stabilize phenotypes towards the mean of the population rather than extremes. Interactive effects of two alleles inherited from one parental strain commonly resulted in phenotypes closer to the population mean than the additive effects from the two loci, and often much closer to the mean than either single-locus model. Alternatively, combinations of alleles inherited from different parent strains contribute to more extreme phenotypes not observed in either parental strain. This class of phenotype-stabilizing interactions has effects that are close to additive and are thus difficult to detect except in very large intercrosses. Nevertheless, we found these interactions to be useful in generating hypotheses for functional relationships between genetic loci. Our findings suggest that while epistasis is often weak and unlikely to account for a large proportion of heritable variance, even small-effect genetic interactions can facilitate hypotheses of underlying biology in well-powered studies.


Author Summary
The role of statistical epistasis in the genetic architecture of complex traits has been of great interest to the genetics community since Fisher introduced the concept in 1918. However, assessing epistasis in human and model organism populations has been impeded by limited statistical power. To mitigate this limitation, we analyzed bone and body composition traits in an unusually large mouse intercross population of over 2000 mice, paired with a recently-developed computational approach that leverages

Introduction
The relevance of epistasis in genetic architecture is yet unresolved. In genetic screens of model systems, the evidence for genetic interaction is abundant [1][2][3][4] and has been proven biologically relevant [5][6][7]. However, the situation is less clear in human populations as epistasis is difficult to detect with confidence due to multiple testing across a high number of variants, underpowered samples, evolutionary history, imperfect model selection, or unmeasured confounding variables or noise [8,9]. While most studies detect only additive variance, recent studies have demonstrated a role of epistasis in the genetics of gene expression [10,11] and occasionally link genetic interactions to disease [12]. Thus the extent to which genetic interactions contribute to unexplained variance or provide biological insight in population-based studies is unclear. One strategy to address this problem is controlled experiments in mammalian model systems, in which genotypes are artificially determined and environmental variation is minimized [13][14][15][16].
In this work, we used a multi-trait strategy to investigate the role of epistasis in regulating complex traits in a large mouse intercross. Bone mineral density (BMD) is a complex trait regulated by the interaction of many genetic and environmental factors [17][18][19], and is the best known surrogate measure of fracture risk in patients with osteoporosis [20][21][22]. Human and rodent studies have implicated many candidate quantitative trait loci (QTL) in influencing BMD [17][18][19][23][24][25][26][27]. Many of these loci are pleiotropic and have been found to influence body weight [28], body fat [29][30][31], and bone size [32] in addition to bone density. Epistatic interactions are also common among loci affecting BMD [26,33,34]. A deep understanding of the genetic regulation of BMD, as well as possible intervention points for therapeutics, requires addressing this complex genetic architecture. Assessing genetic interactions and their relation to multiple phenotypes provides an overall picture of the genetic network regulating BMD and related phenotypes. We used a recently developed method, Combined Analysis of Pleiotropy and Epistasis (CAPE) [35], to integrate information across multiple phenotypes to infer directed genetic interactions between loci.
We integrated genetic interactions influencing BMD, femoral circumference, body weight, and body fat percentage in a mouse intercross population [36][37][38][39]. In particular, we were interested in investigating the genetic architecture of these phenotypes in a population with reduced levels of circulating insulin-like growth factor I (IGF1). IGF1 is a major factor involved in bone development and mineralization [40][41][42][43]. Analysis in a population with reduced IGF1 can reveal aspects of bone density that vary at severely reduced levels of IGF1 [44], thereby unmasking more subtle genetic loci involved in this phenotype [45,46]. To avoid major effects of the IGF1/growth hormone (GH) axis we used mice homozygous for the "little" or lit mutation [47], a null mutation in the gene coding for growth hormone releasing hormone receptor (Ghrhr). GH levels, and consequently circulating IGF1 levels, in mice homozygous for the lit mutation are reduced to about 10% of wild type levels [37,47,48]. These mice also exhibit reduced growth, increased fat mass, and decreased bone mass relative to heterozygotes and wild type mice [49]. Thus lit homozygotes offer the opportunity to study the genetics related to both bone growth and body fat composition in a population in which one of the major hormonal axes regulating these phenotypes is greatly reduced.
A population of 2054 F 2 male and female mice derived from a cross between B6-lit/lit and C3.B6-lit/lit parental strains [36][37][38][39] were analyzed. Compared to B6 mice, C3H mice have 20-30% higher circulating IGF1 levels, higher volumetric bone density, higher rates of bone formation, lower rates of bone resorption, and greater breaking strength of bones [50][51][52]. These strain differences persist in the lit/lit homozygotes [52]. We investigated the genetic interactions influencing body weight, percent body fat, femoral circumference and femoral density in the near absence of one of the major contributors to these phenotypes.

Ethics Statement
All animal procedures followed Association for Assessment and Accreditation of Laboratory Animal Care guidelines and were approved by Institutional Animal Care and Use Committee (The Jackson Laboratory, Protocol #99111).

Mice
Inbred mouse strains used in this study were obtained from our research colonies at The Jackson Laboratory, Bar Harbor, Maine. Mice were produced and housed as described in [76]. Briefly, the mice were housed in same-sex groups of 2-5 per cage in a 14:10 light:dark cycle. The mice had free access to acidified water (pH 2.5 with HCl to retard bacterial growth) and irradiated NIH 31 diet (Purina Mills International, Brentwood, MO).

Construction of Congenic Strain and F 2 Intercross
To investigate heritable factors that control BMD in a model where circulating IGF1 levels are reduced, we used a spontaneous mouse mutation, lit, with a non-functional growth hormone releasing hormone receptor (GHRHR). We generated a congenic strain by transferring the lit mutation from the low-BMD C57BL/6J (B6) strain on which it arose to the high-BMD C3H/ HeJ (C3H) strain by backcrossing for eighteen generations. In both C57BL/6J-Ghrhr lit/lit /J (B6lit/lit) and C3H.B6-Ghrhr lit/lit /J (C3.B6-lit/lit) mice, circulating GH is undetectable, serum IGF1 is low, and femoral volumetric BMD by pQCT, femur length, and body mass are reduced compared to heterozygous lit/+ mice [37,47,48]. Although C3.B6-lit/lit mice are of the same body weight and femur length as B6-lit/lit mice, C3.B6-lit/lit mice have higher BMD. Crosses between B6-lit/lit and C3.B6-lit/lit F 1 mice produced the 1008 male and 1062 female F 2 GH/ IGF1 deficient mice analyzed here.

Genetic Analyses
Mice were genotyped at 100 markers using PCR of oligonucleotide primer pairs (MIT markers, www-genome.wi.mit.edu/cgi-bin/mouse/index) from Research Genetics (Birmingham, AL) as described in [76,77]. The pairs amplified strain-specific sequence length polymorphisms, allowing identification of parental strain of origin. Genotypes at each locus were identified as B6/B6, B6/C3H, or C3H/C3H.

Phenotype Measurements
Body weight and femur length. Anesthetized mice were weighed using a routinely calibrated Ohaus electronic scale. Femur length was measured using a digital caliper (Stoelting, Wood Dale, Ill).
Percent body fat by PIXImus. Data were collected on anesthetized mice using the PIXImus small animal DEXA system (LUNAR, Madison, WI), software version 1.43.036.008 as described in [78]. The machine was calibrated daily with a phantom of known density. Measurements of BMD showed low variability: less than 1% for whole body measurements, and about 1.5% for specialized regions. Body fat (BF) and percent body fat (%BF) were derived from measurements of total body weight (TBW), total lean mass (TLM), and total body mineral content (TBM) as follows: BF = TBW-(LBM + WBM); %BF = BF/TBW.
Femoral BMD and periosteal circumference by pQCT. The XCT 960M was used to measure total BMD at 2-mm intervals as described in [77] and [79]. Periosteal circumference was determined at midpoint of the total femur length. Precision of these measurements was determined to be 1.2% through repeated measurement of a single femur. Hydroxyapatite standards (0.050-1.000 mg/mm 3 ) were used for calibration. The correlation between measured density and actual standard density was r = 0.997.

Combined Analysis of Pleiotropy and Epistasis
CAPE is a strategy that detects epistasis and interprets it in terms of directed enhancing and suppressing influences between genetic loci [81]. It has been released as an R package suitable for mouse intercross studies [35]. The method uses regression on pairs of loci to detect interaction effects from each locus pair on each phenotype. It then combines model parameters across phenotypes to infer directed, QTL-to-QTL influences that replace the interaction effects on each phenotype. The result is a pair of directed coefficients modeling how the two loci influence each other's activity, rather than how each pair independently affects each phenotype.
We began the analysis by using R/qtl [82] to impute psuedomarkers in between each measured marker, increasing the number of markers from 100 to 194. We normalized all traits (body weight, body fat percentage, femoral density, and femoral circumference) using rank Z normalization. We then decomposed the normalized traits using singular value decomposition (SVD) to obtain orthogonal eigentraits (ETs) that combined common signals across all phenotypes. All markers were used in the pair-wise interaction scans. However, we filtered marker pairs tested by linkage disequilibrium (LD) to avoid false positive interactions. We excluded all pairs with genotype Pearson correlation r 0.5. This reduced the number of pairs tested from 19,110 (all pairs from 194 markers and 2 covariates) to 18,576 pairs. Linear regression was then performed on all filtered pairs of markers 1 and 2: where U represents ETs, and is an error term. The index i runs from 1 to the number of individuals, and j runs from 1 to 3 (the number of ETs.) x i is the probability of the presence of the C3H allele for individual i at locus j.
For each pair of markers, the regression coefficients are collected for all ETs and reparametrized to two new terms (δ 1 and δ 2 ). These terms represent the additional activity of each variant when the other is present. For example, δ 1 is the additional effect that variant 1 has on each phenotype when variant 2 is present. It should be noted that the δ terms describe the interaction coefficient between the marker pair independent of phenotype. The δ are determined from the pairwise regression parameters via pseudoinversion as follows: To convert the δ terms to directed influence variables, the following transformation was applied: The sign of the m 12 and m 21 terms indicates how each marker influences the other in terms of enhancement (positive) or suppression (negative). A negative value for m 21 indicates that variant 1 reduces the activity of variant 2 on all phenotypes, whereas a positive value increases the activity. To estimate variances of the new model parameters, error estimates were propagated via second-order Taylor expansion [81,83].
Permutation tests were used to calculate p-values for all model parameters. The pair of markers being tested was permuted together relative to the covariates [84]. By combining permutations across all marker pairs, we saved computation time while generating a single null distribution indistinguishable from the null distribution generated by repeatedly permuting each single pair. This null distribution was composed of 164,679 total permutations representing 3 permutations for each of 54,893 locus pairs.
Main effect significance was also determined through permutation testing. The maximum main effect of each locus across all pairwise contexts was selected as the main effect for that locus. All p-values were corrected for multiple testing using the Holm step-down procedure [85] and all variant-to-ET main effects were translated to variant-to-phenotype effects through multiplication by the singular value matrices from the original SVD.

Grouping Linked Markers
To define distinct QTL regions for the interaction network, adjacent markers were combined into linkage blocks. We calculated the correlation matrix for all markers on each chromosome. We used this similarity matrix as an adjacency matrix to construct a weighted network depicting the similarity between all pairs of markers on a single chromosome. Using the fastgreedy community detection algorithm [86] in R/igraph [87], we then calculated the community membership of the vertices in the network. We assigned adjacent markers in the same community to a single QTL region. This process ensured a robust grouping of markers based on their genotypic similarity. A block was considered to have a significant effect on a phenotype if one or more of the resident markers had a significant effect. After this point, we refer to linkage blocks with significant associations as QTL or loci. See S3 Table for markers included in each block.

Identification of Candidate Genes
Some of the groups of linked markers interacted significantly with IGF1, which is a sufficiently specific interaction to allow a candidate gene search. For each QTL that interacted significantly with IGF1 we generated a list of potential candidate genes by finding all genes in the QTL with annotations for bone density. To determine whether any of these genes interacted with IGF1, we used the Integrative Multi-species Prediction tool [88] to generate a network between the query genes and IGF1. We included the maximum of 50 additional genes and adjusted the confidence of the interactions until IGF1 was included in a subnetwork of greater than two genes. From this subnetwork, we extracted genes in the original QTL region, but which were not necessarily included the original query relating to bone density. Candidate polymorphisms in C3H alleles of these genes were identified using the Sanger SNP database [54]. Furthermore, we used previously published expression data from hepatic tissue of chow-fed B6 and C3H mice [56] to identify genes that were differentially expressed between the strains and pairs of genes that had correlated expression. Because most observed regulation of gene expression is in cis, it is reasonable to assume that the genetic variation that influences gene expression will influence expression levels similarly wherever the gene is expressed. Differential expression was determined with Student's t-test, and significance of correlation was the significance of the Pearson correlation coefficient [89].

Motif Analysis
We examined the enrichment of three-node topological patterns, or network motifs [58], in the set of all interactive genetic models. Although each interaction and main effect was independently derived, we grouped the significant effects into a network to detect general patterns. We combined interactions with main effects to generate three-node motifs, which included two interacting variants and a single phenotype. We counted enhancing and suppressing motifs either with two main effects of the same sign (coherent), or two main effects of the opposite sign (incoherent). To determine whether each type of motif was significantly enriched or depleted, we performed permutations by shuffling the signs of the significant interactions independently of the signs of the main effects. This permutation scheme preserved the topology of the network thereby acknowledging constraints caused by shared edges between motifs. By permuting the main effect edges independently of the interaction edges, we prevented spuriously enriching for enhancing or suppressing interactions simply because there were many negative or positive main effects on a given phenotype. We permuted the edge signs 100,000 times to generate a null distribution for each type of motif. The directions of the interactions were not taken into account for this analysis. We used linkage blocks as interacting units and included all interactions and main effects that were significant at a Holm-corrected p 0.01.

Single-Locus Effects on Phenotypes
We used linear regression to determine the association of each locus with each phenotype (Fig 1,  S4 Table). Each of the phenotypes had multiple associated QTL, and these QTL often overlapped across multiple phenotypes. For example, femoral density and femoral circumference shared a large QTL on Chr 4, and body weight, percent fat and femoral circumference showed overlapping QTL on Chr 17. These overlapping QTL indicate the possibility of common genetic factors underlying multiple phenotypes, such that information can be combined across multiple phenotypes to gain information about individual loci. Unique QTL were also observed, providing nonredundant information to discern genetic factors with phenotypic specificity.

Single-Locus Effects on Eigentraits
We decomposed the normalized, mean-centered phenotypes into eigentraits (ETs) (Fig 2A). The first ET represented an average of femoral circumference, percent body fat and body weight and captured 52% of the overall variance. The second ET represented the contrast between femoral circumference and femoral density and captured 25% of the total variance. The third ET captured 19% of the total variance and represented a contrast between femoral circumference and percent body fat. The fourth ET captured only 4% of the total variance, and described a contrast between body weight and percent fat. Because this ET captured a small amount of the total variance, did not include strong contributions from the bone phenotypes, and may add noise to the analysis, we excluded it from this analysis.
Single-locus associations with each ET detected multiple QTL (Fig 2B). Since ET are linear combinations of traits, each QTL indicates a potentially pleiotropic association with varying effect strengths on each trait. For example, data for body weight, percent fat and femoral circumference had overlapping QTL on Chr 17. These phenotypes also contributed substantially to ET1, and there was a corresponding significant QTL for ET1 on Chr 17 representing the common QTL.

An Extensive Network of Weak Genetic Interactions Influences Bone and Body Composition
CAPE analysis of the first three ETs produced a large network of genetic interactions (Fig 3). The high-confidence network (p 0.0005) consisted of a single connected component including 67 QTL linked by 102 directed interactions. Each interaction was directed from a source locus to a target locus and was either suppressing (negative), or enhancing (positive). In suppressing interactions the presence of the C3H allele at the source locus reduced the phenotypic effect of the C3H allele at the target locus regardless sign of the main effect. In enhancing interactions the presence of the C3H allele at the source locus increased the phenotypic effect of the C3H allele at the target locus regardless of the sign of the main effect. The QTL were distributed across the four phenotypes as follows: body weight had 13 QTL; percent fat had 11 QTL; femoral density had 19 QTL; and femoral circumference had 24 QTL. Among the interactions between QTL, 29 were suppressing and 73 were enhancing.
We also used standard linear regression to assess the effect of each marker pair interaction on each normalized phenotype. We calculated the pairwise interaction coefficients for all marker pairs including the two covariates, sex and circulating IGF1 levels. We calculated empirical p-values from 50,000 permutations and corrected the p-values for multiple testing using the Holm step-down procedure. No marker-marker or marker-covariate interactions were significant using the standard epistasis analysis after correction for multiple testing.
We found several significant sex-QTL interactions, all of which were enhancing (Fig 3C). From the single-locus regressions, potential sex-interacting loci were seen on Chrs 1, 6, 7, 10, and 14 (Fig 1). We tested these conditional associations directly and confirmed genetic interactions with sex on Chrs 1, 7, 10 and 14. For example, there was a larger sex difference in femoral density among animals homozygous for the C3H allele at the Chr 1 locus than among animals homozygous for the B6 allele (Fig 4). Thus the C3H status at this locus enhanced the negative effect of the male sex on femoral density, as well as the positive effect of the male sex on circumference, giving females with this genotype increased bone density and reduced circumference relative to males.
An apparent locus on Chr 6 (Fig 1) was not identified as interacting significantly with sex because a consistent directional model could not be fit across all phenotypes. Several markers on Chr 6 had relatively large interaction coefficients in the linear regression with sex, but none of these interactions were significant after correction for multiple testing.

Residual IGF1 Interacts with Genetic Loci
Like sex, IGF1 had significant main effects on all phenotypes. Higher levels of circulating IGF1 were associated with higher body weight, body fat percentage, femoral density, and femoral circumference. IGF1 was also found to interact significantly with several QTL (Fig 3C). Interestingly, all interactions from IGF1 to loci were suppressing, i.e. IGF1 suppressed the effects of these loci. For example, one Chr 9 QTL had negative effects on both femoral density and femoral circumference when IGF1 levels were low, but at high levels of IGF1 this effect was suppressed, and the QTL had a positive effect on these phenotypes (Fig 5). Conversely, this locus also enhanced the effects of IGF1. Looking at body weight as a function of IGF1, for example, it Phenotypes illustrating sex effects enhanced by a QTL on Chr 1. Males have higher body weight, body fat, and femoral circumference, and lower femoral density relative to females. The Chr 1 QTL clearly enhances these effects on femoral circumference, which is evident in the larger phenotypic effect of sex in mice homozygous for the C3H allele at the Chr 1 locus. There are smaller interaction effects for percent fat and femoral density. Importantly, the directed model for this interaction describes the results across all four phenotypes. can be seen that in the B6 homozygotes, IGF1 levels had a positive effect on body weight. There was a slightly larger effect in the heterozygotes, and the largest positive effect of IGF1 on body weight was in the C3H homozygotes (Fig 5). This general pattern is replicated across all phenotypes.
The molecular specificity of an interaction with IGF1 offered the opportunity to search for candidate genes in QTL interacting with IGF1, even though the regions were large (Methods). The second QTL on Chr 10 yielded promising candidates. Using MouseMine, we found 19 genes in the region with annotations to bone density [53]. We used these genes as a query gene set in the IMP tool. The IMP tool determines the likelihood that pairs of genes in a query set interact. It uses databases of known physical interactions, genetic interactions, and correlated expression to pull in additional genes though which genes in the query set may interact. The result is a network of high-confidence interactions that relate the query genes to each other. QTL interaction with IGF1. The effects of IGF1 are enhanced by QTL 9.2. IGF1 has a small positive effect on each phenotype in animals homozygous for the B6 allele at this locus. This can be seen in the positive slope of the blue line in each panel, which shows that each phenotype increases with increasing levels of IGF1. The effect of IGF1 on all phenotypes except femoral density is enhanced in heterozygotes (green line) and enhanced further in C3H homozygotes (purple line). For these genotypes, the increase in slope indicates that there is a more pronounced increase in phenotype as IGF1 levels increase. Y-axes represent rank normalized values of each phenotype, and x-axes represent rank normalized IGF1 levels scaled to range between 0 and 1. Shaded regions show 95% confidence interval for each slope. IMP found a network of 45 genes that linked IGF1 with four bone density-related genes (Esr1, Nr1h4, Kitl, Ctgf) from the query gene set. The minimum confidence for interactions between gene pairs in this network was 91%.
All four genes contain SNPs predicted to be functionally relevant between B6 and C3H [54], including splice site variants, missense variants and variants in regulatory regions (S5 Table). One of the genes, Kitl was differentially expressed in hepatic tissue between C3H and B6, with C3H mice showing lower expression (Student's t-test, p = 0.012) [55,56]. Kitl expression also showed a strong negative correlation with Igf1 expression (r = −0.67, p = 0.016) [55,56]. These findings allow us to hypothesize that Kitl is a potentially causative gene in QTL 10.2. Another potential candidate in this region, Bicc1, was recently found to be related to bone density in mice [57]. However, it was not differentially expressed in the hepatic tissue of C3H and B6 mice (Student's t-test, p = 0.25) (S5 Table), and thus we consider it a less likely candidate in the context of this study. Other candidate regions did not reveal promising causative candidates when analyzed using the same methods.

The Network Is Enriched for Stabilizing Motifs
In addition to individual interactions, we examined the enrichment of three-node patterns, or network motifs [58] (Methods). Each motif consists of two interacting genetic loci and an affected phenotype. At significance p 0.01, we identified a total of 116 motifs influencing body weight, 84 motifs influencing percent fat, 132 motifs influencing femoral density, and 274 motifs influencing femoral circumference (Fig 6). Motifs are classified as coherent, i.e. the main effects are of the same sign, or incoherent, i.e. the main effects are of opposite signs. In addition, motifs can be either suppressing, meaning that there is a suppressing interaction between the genetic loci, or enhancing, with an enhancing interaction between the loci (Fig 6). These different classes of motifs may speak to the general structure of the underlying biological interactions [59][60][61][62][63]. For example, a motif with coherent main effects and a suppressing interaction between them indicates genetic redundancy and may result from proteins operating in series in the same pathway or physical interaction between gene products [60,61]. Alternatively, a synergistic interaction exists in a coherent motif with an enhancing interaction between the variants. In this case, the variants and the interaction between them all push the phenotype in the same direction. Such an interaction may indicate regulatory interactions between parallel regulatory pathways that affect the same process [60].
These motifs can be divided into those that stabilize phenotypes and those that destabilize phenotypes. For example, incoherent enhancing motifs tended to have a stabilizing effect on phenotype because the main effects drove the phenotype in opposite directions. This was the largest class of motif represented in our network (42% of total motifs detected) and was significantly enriched across all phenotypes except femoral density (Fig 6). In contrast, coherent enhancing motifs tended to be destabilizing. The main effects of the interacting loci both drove the phenotype in the same direction, and the enhancing interaction between these loci drove the phenotype further in the same direction, generating extreme phenotypes. These enhancing coherent motifs made up only 17% of the total motifs detected and were significantly depleted across all phenotypes in our network (Fig 6).
Motifs with suppressing interactions were more evenly distributed. Coherent suppressing motifs, which tended to stabilize phenotypes, made up 33% of the total motifs and were enriched for all phenotypes except percent fat. Incoherent suppressing motifs, which tended to be destabilizing, were the smallest class of motif (8%) and were enriched in body weight and femoral circumference (Fig 6). In summary, the network overall consisted of weak interactions (mean β = 0.25±0.16) compared with the additive effects (mean β = 0.52±0.36), and the majority of the interactions (59%) were enhancing. Further, the two largest motif classes, making up 75% of the total motifs, tended to be stabilizing.

Discussion
Bone density and body composition phenotypes, such as percent body fat, are complex traits influenced by many genetic variants, both shared and distinct. Here we used combined analysis of pleiotropy and epistasis (CAPE) to investigate how genetic loci interact in a large population of mice to influence femoral density, femoral circumference, body weight, and percent body fat. Using this exceptionally well-powered mouse intercross we detected many main-effect and interacting QTL associated with these traits and found an extensive network of genetic loci influencing the four phenotypes. These genetic interactions were not detectable through standard regression-based epistasis analysis. We were also able to infer both the directionality and sign of the interactions, which improved our ability to identify candidate QTL genes and provided a uniquely broad view of the genetic architecture of bone and body composition phenotypes. One of the notable features of the network was an asymmetric distribution of enhancing and suppressing interactions, which was apparent for interactions between QTL, as well as for QTL interactions with sex and circulating IGF1. Sex showed a strong bias in both the directionality and sign of interactions (Fig 3B). Most (8 of 10) interactions were QTL that enhanced the phenotypic effects of the male sex on each phenotype. For example, the presence of the C3H alleles in QTL 1.1 increased the positive effect that the male sex had on weight, percent fat, and femoral circumference, as well as the negative effect the male sex had on femoral density. The remaining interactions (2 of 10) showed effects in the opposite direction. For example, being male enhanced the positive main effect that QTL 4.2 had on femoral circumference and density (Fig 3C). Sex and sex hormones are known to influence bone growth and density [64][65][66] and genetic loci that interact with sex to influence bone phenotypes have also been previously identified in rodent models [23,27,33,67,68]. That the QTL were the sources and sex was the target of the majority of these interactions highlights how CAPE determines directionality and interpretation of interactions. Sex was widely pleiotropic, affecting all phenotypes significantly (Fig 3A). An interaction in which sex is a target implies that the QTL influences all of these phenotypes to be identified by CAPE via its modifications on sex. Conversely, sex can target a non-pleiotropic QTL to influence individual phenotypes or a subset of phenotypes. It is possible that the sex-enhancing QTL contain variants in endocrine genes that globally affect sex effects. The two sex-enhanced QTL (4.2 and 13.1), which reciprocally enhance sex effects, are loci for which the uniform enhancement of the sex effects was insufficient to fit all phenotypes simultaneously. These QTL are therefore more likely to be involved in processes that differentially affect the phenotypes, suggesting more specific roles in each phenotype that are responsive to sex difference. These findings imply that interactions between sex and QTL are commonly due to genes that lie "upstream" of processes that underlie sexual dimorphism, rather than "downstream" genes with functions that are altered by sex hormones.
Circulating IGF1, which is reduced to 10% of wild type levels in this lit/lit population, had both main effects and interaction effects. It suppressed the effects of four loci, whereas four loci enhanced the effects of circulating IGF1 (Fig 3C). In contrast to sex, the interactions that IGF1 participates in are relatively balanced between incoming and outgoing interactions. This balance may reflect the fact that, unlike sex, IGF1 is a specific protein that physically interacts with other proteins. We interpret the role of the loci suppressed by IGF1 to be compensatory pathways influencing bone density when IGF1 levels are extremely low. At higher levels of circulating IGF1 the effects of these loci are diminished because the role of the causal genes becomes less relevant. For example, QTL 7.1 has a positive main effect on femoral density, which is suppressed by the presence of density-promoting IGF1. Conversely, the loci that enhance the effects of IGF1 may be targets of IGF1 that act to increase bone density and other phenotypes when IGF1 is present. These QTL, e.g. QTL 10.2 discussed below, can therefore be interpreted to contain genes in pathways that regulate and/or respond to IGF1 signaling. We note that one QTL, Chr 9.2, acts as an enhancer of IGF1 and is suppressed by IGF1 (Fig 3C). This QTL may therefore contain multiple genes involved in both IGF1 pathways and compensatory pathways, or be the result of a gene with an IGF1 signaling role that also serves to trigger alternative pathways in the absence of circulating IGF1.
Although our QTL intervals were too large to decisively identify positional candidate genes, results for main effects and genetic interactions can be combined with prior data to reason about potential candidates. As an example, our model for QTL 10.2 interaction with IGF1 illustrates how hypotheses of causal QTL genes can be generated by requiring consistency in both main effects and interactions. Of the genes in QTL 10.2, Kitl had the best evidence to suggest a role in interacting with IGF1 to influence bone density. In our study we found that the C3H allele at QTL 10.2 had a negative main effect on femoral density. Prior work determined that hepatic Kitl expression is lower in C3H mice than in B6 mice (Student's t-test, p = 0.012) [56] ( Fig 7A) and that low expression is potentially associated with low femoral density [57]. Taken together these data suggest that the C3H allele of Kitl may be a loss of function variant that reduces Kitl transcript levels and consequently femoral density. Furthermore, this hypothesis can account for the directional interaction between QTL 10.2 and IGF1, in which the C3H allele at QTL 10.2 enhanced the positive effects of IGF1 on femoral density. The hepatic expression study found that Kitl and Igf1 are negatively correlated (r = −0.67, p = 0.016) (Fig 7) [56] (Fig 7B). We can therefore hypothesize that a decrease in Kitl expression from the C3H allele corresponds to a rise in IGF1 activity and consequently its positive effects on femoral density. Combining the evidence for an interaction between Kitl and circulating IGF1 with the main effects of Kitl and QTL 10.2 suggests that increased Kitl transcript acts to reduce IGF1 activity in the reference B6 genotype (Fig 7C). When the C3H allele is present Kitl expression decreases, allowing residual Igf1 transcript levels to remain relatively high and thereby enhancing the effect of IGF1 on femoral density in C3H mice relative to B6 mice. While speculative, this capacity to generate specific molecular hypotheses by combining genetic interactions with prior molecular results illustrates the importance of genetic interactions in hypothesis generation, even when they are a minor correction to additive effects. Such analysis is expected to be especially effective in a study with greater genetic mapping resolution and fewer candidate genes per QTL.
Our interaction network was derived in a population homozygous for the lit mutation, a receptor variant that perturbs IGF1 levels. One possible consequence is that the measured phenotypes have been decanalized and some of the QTL we observed may correspond to cryptic genetic variation [45,46,69] that only influence traits in the presence of the lit/lit mutation. This release of cryptic variation may be due to IGF1 effects being reduced to the point that minor genetic effects become observable. Alternatively, low levels of IGF1 induce activity in pathways that are not used in non-lit mice for maintaining growth and bone density, thereby making variation in the genes in these pathways relevant [46]. Potentially cryptic variants could be identified as those which do not replicate in a similar analysis of standard B6 and C3H strains, although resolution of compensatory pathways would likely require a conditional Igf1 knockout.
Overall, the QTL-QTL interaction network in this study had significant enrichment of enhancing interactions (Fig 6). Rather than instances of classic genetic synergy in which two variants combine to amplify a common effect, these enhancing interactions were mostly between variants with incoherent (opposing) main effects. We interpret these interactions as a signature of similar phenotypes between B6 and C3H strains that arise from different combinations of alleles. When alleles from two different strains are recombined in novel ways, unanticipated variation is introduced and extreme phenotypes result (Fig 8). Thus our enhancing interactions indicate a reduction of extreme phenotypes when both loci are homozygous for either parental allele. This moderating effect only occurs if the interaction between incoherent variants is enhancing. For any given phenotype an enhancing interaction from a positive main effect to a negative main effect is equivalent to a suppressing interaction in the reverse direction (Fig 6). However, when this model is applied across multiple phenotypes it is more likely to lead to high phenotype variability in the double homozygotes. This is because the suppressing QTL reduces an opposing main effect thereby favoring its own main effect. In contrast, in the enhancing model the enhancing QTL increases the opposing effect, thus balancing the effect of the overall interaction and bringing the phenotype toward the parental mean.
Although additive genetic variance is usually sufficient to equalize opposing QTL effects, when the positive and negative main effects are of unequal magnitude, genetic interactions are The same regulatory network in a C3H mouse. Relative to the B6 alleles, the C3H alleles have negative, positive, or neutral effects. The C3H alleles combine additively and interactively to balance main effects and yield a phenotype similar to that of the B6 mouse. Taken in pairs, the most common outcomes will be: (i) balanced opposite-effect alleles that additively mimic the reference phenotype; (ii) imbalanced opposite-effect alleles that approach the reference phenotype, thereby revealing an incoherent enhancing interaction; or (iii) same-effect alleles will redundantly affect the phenotype either positively or negatively, thereby exhibiting a coherent suppressing interaction. Synergistic interactions (coherent enhancing or incoherent suppressing) are relatively rare in this model. (C) These network properties are manifest as pair-wise genetic interactions that reduce phenotypic variance when both loci share ancestry, as shown for femoral circumference. Mean phenotypes are plotted by genotype for each incoherent enhancing interaction (N = 141). Genotype refers to the homozygous genotype at each locus, either the source locus or the target locus. Heterozygous animals, which display intermediate phenotypes, are excluded for clarity. Source and target loci are defined by CAPE directional analysis. Interactions with a negative-effect source variant and positive-effect target variant are plotted in blue, with the opposite configuration in brown. Doublehomozygote animals exhibit a reduced range of phenotype values than those with different alleles at the two loci, which have a broader range of phenotypes.
doi:10.1371/journal.pgen.1005805.g008 required to stabilize the phenotypes at the levels seen in the founders. Our network outcomes suggest that genetic interactions in intercross populations will commonly have relatively weak effects, since they are a fractional correction to main effects. This can be contrasted with extreme phenotypes that arise from synergistic (coherent enhancing) interactions in which interaction effects are of greater magnitude than main effects, such as in the classic case of genetic buffering [70]. These interactions are significantly depleted in our network.
Enrichment for motifs that drive phenotypes toward founder values is consistent with previous observations of phenotypic variation in recombinant populations. In populations of mice and Drosophila with introgressed genomic regions, simple additive contributions from all variants would result in phenotypic variance orders of magnitude greater than is observed between founders. That the founders have reduced phenotypic variance between them suggests that most interactions are less than additive [71][72][73][74]. Here, we see an enrichment of motifs leading to a reduction of extreme phenotypes, namely, suppressing coherent motifs, and enhancing incoherent motifs. While these effects are often weaker than main effects and therefore may not substantially improve the heritability accounted for, they nevertheless indicate genes acting together in a common pathway or process [60,61] Our analysis of a large intercross population has revealed a number of features that may be generalized to the genetic architecture of complex traits. First, we have found that a sufficiently powered study paired with a multi-trait analysis method can reveal a large network of genetic interactions between QTL. The systematic patterns in this network, including interactions with sex and a molecular marker, suggest that the interactions are signatures of the pathways and processes involved in the regulation of complex physiological traits. Second, we have found that the most common type of genetic interaction is a fairly subtle signal arising from allelic combinations that drive phenotypes towards median rather than extreme values. These interactions are either minor deviations from additivity or involve alleles with redundant effects, with the former being particularly difficult to detect in all but the largest study populations. These findings are consistent with recent work on a very large meta-analysis of twin studies [75]. Third, we note that the genetic interactions detected here form a connected network involving many interactions between the same subsets of QTL. We speculate that this is because the casual variants reside in groups of co-functional genes that compose specific pathways or processes, and that these pathways vary at multiple points between the B6 and C3H inbred strains. Since pair-wise combinations of C3H alleles have been shown to interactively drive the phenotypes towards the median, we speculate that as more genomic regions from C3H are inherited in a single individual, higher-order combinations of C3H alleles within these pathways will further canalize toward the C3H phenotype rather than cause large phenotypic variation. This may contribute to higher-order epistasis, since many of the variants will have strong effects in isolation that vanish in combinations. In sum, these findings suggest that the most common forms of epistasis may often be difficult to detect, and that the analysis of genetic interactions is nevertheless a powerful means to understand the genetics of the underlying biological pathways and processes.
Supporting Information S1 Table. Data used in this analysis formatted for CAPE. A comma-separated file containing phenotypes and genotypes used in this analysis. The first six columns contain covariates (sex and IGF1) as well as the phenotypes. The remaining columns contain the name, chromosome, position, and genotypes for 100 MIT markers. (CSV) S2 Table. A comma-separated file containing information for all tested pairwise interactions and main effects. The table is sorted by decreasing effect size, and each row contains the names of the source and target markers (or maker and phenotype for main effects), the effect size of the interaction (or main effect), the standard error, the standardized effect size, the empirical p-value as calculated from permutations, and the Holm-corrected p-value. (CSV) S3 Table. A tab-delimited file containing the markers assigned to each linkage block. The first column displays the name of the linkage block, and subsequent columns show the names of the markers in the linkage block. All marker names are prefixed with a chromosome label. Markers designated "locX" are imputed pseudomarkers obtained from R/qtl (Methods). (TXT) S4 Table. Single-Marker Scan Information by Linkage Block. A tab-delimited file containing information about linkage block effects in the single-marker scan. Information includes the number of markers in each linkage block, genomic start and stop positions, standardized effect sizes and significance levels for each eigentrait, as well as a gene count for the block. Only protein coding genes were included in the count. (TXT) S5