The Genetic Architecture of a Congenital Heart Defect Is Related to Its Fitness Cost

Abstract

In newborns, severe congenital heart defects are rarer than mild ones. This epidemiological relationship between heart defect severity and incidence lacks explanation. Here, an analysis of ~10,000 Nkx2-5^+/− mice from two inbred strain crosses illustrates the fundamental role of epistasis. Modifier genes raise or lower the risk of specific defects via pairwise (G×G_Nkx) and higher-order (G×G×G_Nkx) interactions with Nkx2-5. Their effect sizes correlate with the severity of a defect. The risk loci for mild, atrial septal defects exert predominantly small G×G_Nkx effects, while the loci for severe, atrioventricular septal defects exert large G×G_Nkx and G×G×G_Nkx effects. The loci for moderately severe ventricular septal defects have intermediate effects. Interestingly, G×G×G_Nkx effects are three times more likely to suppress risk when the genotypes at the first two loci are from the same rather than different parental inbred strains. This suggests the genetic coadaptation of interacting G×G×G_Nkx loci, a phenomenon that Dobzhansky first described in Drosophila. Thus, epistasis plays dual roles in the pathogenesis of congenital heart disease and the robustness of cardiac development. The empirical results suggest a relationship between the fitness cost and genetic architecture of a disease phenotype and a means for phenotypic robustness to have evolved.

1. Introduction

Congenital heart disease (CHD) is recognized as the most common birth defect, but cardiac malformation actually encompasses anatomically distinct phenotypes in which the severity and incidence of a defect are inversely related. For example, secundum atrial septal defects (ASD), which are compatible with survival into middle age, are three times more common than atrioventricular septal defects (AVSD), which can cause death during infancy (Figure 1a) [1,2]. No simple or single hypothesis explains the epidemiological relationship, such as a genotype–phenotype correlation or a lower frequency of the causes of severe defects. Furthermore, the fetoplacental circulation enables severe defects to be compatible with survival to birth. The inverse relationship does suggest that each defect has a liability threshold related to its fitness cost: the higher the threshold, the rarer a defect is. In other words, the developmental pathways that lead to a severe defect are more robust to perturbation than the ones to a mild defect.

Questions regarding how the fitness cost of a defect relates to the robustness of the underlying pathways are practical and fundamental. An appreciation of the genetic architecture of a heart defect—i.e., the genes involved, alleles, interactions, and their effects—could inform genetic risk models and strategies for prevention. The genetic architecture encompasses more than the proximate cause of disease. As a rule, the same deleterious mutation in different individuals can manifest as a mild, moderate, or severe defect or no defect at all. The relationship holds in general for causes of syndromic and isolated CHD. For example, despite the well-known associations of trisomy 21 and 22q11.2 deletion with atrioventricular septal defects and tetralogy of Fallot, respectively, their most common cardiac phenotypes are structurally normal heart or an atrial or ventricular septal defect [3,4]. Likewise, the NKX2-5 mutation, which causes isolated CHD, is most commonly associated with atrial septal defects, although severe defects such as tetralogy of Fallot and double outlet right ventricle were reported [5,6].

Modifier genes affect the phenotypic outcome [7]. Genetic heterogeneity and rare alleles, however, hamper the detection of epistasis in humans and the elucidation of complex genetic principles. Inbred strain crosses in mutant mouse models can circumvent these limitations. The number of alleles at a locus is limited to two, and the frequencies of both alleles can be kept equally common.

Inbred strain crosses of Nkx2-5^+/− mice, a model of non-syndromic human CHD, demonstrate that the incidence of specific defects can vary between genetic backgrounds [8]. We previously began two series of inbred strain crosses to establish general conclusions that would not be unique to one cross. First-generation (F1) hybrids of the inbred strains C57BL/6N and FVB/N or A/J were found to have a lower incidence of congenital heart defects than Nkx2-5^+/− mice in the C57BL/6N background or the second-generation (F2) offspring of backcrosses to the parental strains or intercrosses (F1×F1). The results suggest a role in CHD genetics for heterosis, a well-described phenomenon in agricultural genetics in which heterozygosity across the genome increases the average level of robustness in a population [8]. Furthermore, these and other crosses demonstrate that genetic interactions modulate the robustness of specific cardiac developmental pathways to Nkx2-5 mutation and hence the risks of anatomically distinct malformations [9,10].

Epistasis contributes quantifiably to fitness traits in simple organisms, such as yeast [11], Drosophila [12], and C. elegans [13]. On the other hand, quantitative evidence for a contribution of epistasis in higher organisms is scant, and theoretical analyses arrive at disparate conclusions regarding the significance of epistasis in complex traits and disease [14,15]. Why the results conflict is hotly debated [16]. One reason pertinent to congenital heart disease is that the contribution of epistasis could vary with the fitness cost of a phenotype [15]. Natural selection indisputably eliminates deleterious mutations, but selection may also increase the robustness of pathways to perturbation by shaping genetic interactions and networks [17].

We reasoned that a comparative analysis of the genetic architecture of four CHD phenotypes caused by the same mutation could illuminate the subject. The four defects—ASD, muscular and membranous ventricular septal defects (VSD), and AVSD—have mild to severe fitness costs, as defined by the likelihood of survival to reproductive age in humans [2]. We determined the number and effects of quantitative trait loci (QTLs) on the risk of each heart defect in Nkx2-5^+/− mice from two different inbred strain crosses. Each QTL represents a pairwise interaction between a modifier gene and Nkx2-5 (G×G_Nkx). In addition, we examined the non-linear effects of higher-order interactions between two QTLs and Nkx2-5 (G×G×G_Nkx). Empirical data are sorely required regarding the nature of genetic risk in CHD. More fundamentally, a quantitative analysis of epistasis in a set of phenotypes with varying fitness costs could inform a longstanding debate regarding the origins of mutational robustness [18].

2. Materials and Methods

2.1. Mouse Strains and Crosses

Inbred C57BL/6N (B) and FVB/N (F) mice were purchased from Charles River (Wilmington, Massachusetts) and A/J (A) from the Jackson Laboratory (Bar Harbor, Maine). Mice were housed under standard conditions in the same room with access to water and chow ad libitum. The Nkx2-5^+/− mutant line was maintained in the C57BL/6N background [8,19]. To produce F1 hybrids, Nkx2-5^+/− C57BL/6N males were crossed to wild-type A/J or FVB/N females. The Nkx2-5^+/− F1 was then crossed to produce A×B and F×B F2 progeny. We also produced F10 and F14 advanced intercross generations of F×B by random mating from the F2 generation onward. No matings between siblings or first cousins were allowed from the F3 and F4 generations onward [20]. Breeder mice were selected at random; they presumably had no significant congenital heart defect by virtue of their healthy appearance and reproductive fitness. The decision to develop the advanced intercross in F×B rather than A×B was arbitrary. An advanced intercross of A×B was not pursued because of resource constraints.

2.2. Collection and Phenotyping of Hearts

We collected and phenotyped hearts as previously described [8]. Briefly, pups were collected within hours of birth. Thoraxes were fixed in 10% neutral buffered formalin. After Nkx2-5 genotyping, all Nkx2-5^+/− and a subset of wild-type hearts were dissected and embedded in paraffin. Hearts were completely sectioned in the frontal plane at 6 µm thickness. Every section was collected, stained with hematoxylin and eosin, and evaluated by two individuals. Defects were diagnosed by the morphology of the atrial and ventricular septae, semilunar and atrioventricular valves, chambers, and the anatomic relationships between them. Secundum ASDs were distinguished from a patent foramen ovale based on the size of the latter in wild-type hearts.

2.3. Single Nucleotide Polymorphism Genotyping

Genomic DNA was isolated by phenol–chloroform extraction. SNP genotyping was performed on affected and normal Nkx2-5^+/− mice. In the F×B advanced intercross, normal controls were randomly selected from unaffected siblings.

The A×B and F×B F2 intercrosses were genotyped using custom panels of ~120 SNPs on the Sequenom MassARRAY system [10]. The average distance between markers was 19–24 Mb. The F×B advanced intercross generations were genotyped on the Mouse MD Linkage Panel (Illumina, San Diego, CA). A subset of 761 SNPs on the panel was selected for subsequent analyses. We eliminated SNPs that had >10% missing genotypes or low median quality scores (GC < 0.35) or were not in dbSNP (build 141). Fewer SNPs were genotyped in the F2 than the advanced intercross generations because there is less recombination. Combining SNPs from the F×B F2 and advanced intercrosses yielded 887 markers spaced an average of 3 Mb apart. SNP locations were assigned according to GRCm38/mm10. The sex chromosomes were not analyzed. Table S1 lists the SNPs and their minor allele frequencies in each cross.

2.4. Genetic Imputation and Association Analyses

QTLs that modify the risk of a heart defect were mapped via case-control association tests. Each type of heart defect was analyzed separately. The same normal controls were used in each analysis.

We utilized R/qtl (version 1.40-8) to map QTLs under a binary trait model in the A×B F2 intercross [21]. Permutation testing (N = 5000) was performed to determine significance thresholds. Significant and suggestive loci had LOD scores >3.3 (α = 0.05) and >2.7 (α = 0.2), respectively.

We implemented a univariate linear mixed model in GEMMA (version 0.94.1) for a combined analysis of the F×B F2 and advanced intercross generations [22]. A likelihood ratio test compared the alternative hypothesis that an SNP has a non-zero effect against the null hypothesis of zero effect. Each binary phenotype was modeled as a function of SNP genotypes and a random or polygenic effect that accounts for relatedness (Equation (1)):

$$\mathrm{y}=\mathsf{\mu }+\mathrm{x}\mathsf{\beta }+\mathrm{u}+ϵ$$

(1)

where y is the vector of binary phenotypes (0 = control, 1 = heart defect), μ is the mean, x is a vector of SNP marker genotypes, β is a vector of marker effects, u is a vector of random effects, and ϵ is a vector of errors. Genome-wide significance thresholds were Bonferroni corrected. The effective number of tests was 542 after accounting for linkage between SNPs [23,24]. Significant and suggestive loci had nominal P values < 9.2 × 10⁻⁴ (0.05/542) and <1.8 × 10⁻³ (1/542). The combined analysis of the F2 and advanced intercross generations enhances the resolution of mapping but poses two special considerations [25].

First, some SNPs were not genotyped in either the F×B F2 or advanced intercross generations and were imputed. In the F2, 758 SNPs were imputed using R/qtl [21]. In the advanced intercross, 126 SNPs were imputed using QTLRel (version 0.2–15). The F10 and F14 were imputed separately because of differences in recombination patterns between the generations [26]. We estimated that >80% of missing genotypes in the F2 were imputed with >70% confidence; only 2% were imputed with less than 50% confidence. For the advanced intercross generations, 99.9% of missing genotypes were imputed with >70% confidence.

Second, to minimize type I error, we accounted for unequal relatedness between individuals in the advanced intercross generations. The univariate linear mixed model, Equation (1), incorporates a genetic relatedness matrix (GRM) that captures the phenotypic variance due to the genome-wide similarity between individuals. For each case-control analysis, we used GEMMA to estimate centered GRMs using the SNP genotypes from the combined F×B generations. Because F2 individuals are essentially full siblings, their pairwise relatedness coefficients were set to 0.5 and 1 on the diagonal in the GRM. We then adjusted the entire matrix (F2 + advanced intercross generations) to the nearest positive definite matrix using the nearPD function in the “Matrix” R package (version 1.2-6). The adjustment adds a small amount of noise to the GRM to remove negative eigenvalues, which is necessary for the stability of the mixed model. A Mantel test for the similarity between the pre- and post-adjusted GRMs showed >95% similarity in all cases.

2.5. Estimation of Heritability

The proportion of the phenotypic variance explained (PVE) by all genotyped SNPs was estimated using restricted maximum likelihood (REML) in the linear mixed model in GEMMA. This method uses the GRM to estimate the genome-wide PVE and its standard error for each type of heart defect in the A×B and F×B intercrosses [27]. PVE estimates were corrected for ascertainment bias and transformed to the liability scale (PVE_l) according to Equation (2):

$${\mathrm{PVE}}_l=\mathrm{PVE}\times \frac{{\mathrm{P}}^2{\left(1-\mathrm{P}\right)}^2}{\mathrm{A}\left(1-\mathrm{A}\right){\mathrm{z}}^2}$$

(2)

where P is the population disease incidence, A is the incidence in the ascertained cohort, and z is the height of the standard normal distribution at the population liability threshold [28].

2.6. Analysis of G×G_Nkx Effects

The proportion of the phenotypic variance explained by an SNP, PVE_SNP, was calculated according to Equation (3):

$${\mathrm{PVE}}_{SNP}=\frac{\mathrm{MAF}\times \left(1-\mathrm{MAF}\right)\times {\mathsf{\beta }}^2}{\mathrm{PV}}$$

(3)

where MAF is the minor allele frequency, PV is the phenotypic variance, β is the effect of the SNP as estimated in GEMMA for both the A×B and F×B intercrosses. PVE_SNP estimates were corrected and transformed to the liability scale in the same manner as the genome-wide PVE estimates.

We estimated the odds ratios for the G×G_Nkx effects using logistic regression models. The odds ratios for each A×B locus were obtained using the base glm function in R. The odds ratios for each F×B locus were obtained using the logistic mixed model implemented in GMMAT [29]. The mixed model controlled for relatedness via the GRMs above.

2.7. Analysis of G×G×G_Nkx Effects

The contribution of higher-order epistatic effects to the risk of a heart defect was calculated according to the physiological epistasis model [30]. The G×G×G_Nkx effect is determined by the difference between the observed and the expected incidences of a defect at a two-locus genotype under the null model of independently acting loci. The expected incidence is calculated from the unweighted regression of the two-locus genotype incidences onto the two single-locus incidences. To compare G×G×G_Nkx effects between defects that have different incidences, the effects were divided by the unweighted average incidence of the defect across all nine genotype combinations between two loci. Effects that raise or lower risk are deemed positive or negative, respectively. To assess statistical significance, we arcsine-transformed the effects to obtain unbiased variances [31]. Using these variances, we conducted two-tailed T-tests to identify significant, non-zero effects. The significance threshold was Bonferroni-corrected for the number of two-locus genotypes in a pairwise analysis, i.e., P = 0.05/9x, where x is the number of pairwise combinations of loci for a defect. To compare effects between defects that have different incidences, the significant epistatic effects were divided by the unweighted average incidence of each defect [30]. Because only one G×G_Nkx modifier locus was identified for ASD and muscular VSD in the F × B population, we assessed G×G×G_Nkx effects in combination with the next most significant SNP (ASD: rs13478997, Chr. 6, P = 1.94 × 10⁻³; muscular VSD: gnf13.115.241, Chr. 13, P = 2.39 × 10⁻³.)

2.8. Statistical Analyses

GEMMA was run on the Ubuntu operating system (version 14.04 LTS). All other statistical analyses were performed in the R Statistical Computing environment (version 3.3.1). All t- and z-tests for comparisons of means and proportions, respectively, were two-sided. We used the R package “ppcor” [32] to calculate Kendall’s (non-parametric) partial rank correlation τ with a two-sided alternative hypothesis. Partial correlations permitted us to calculate non-parametric correlation while controlling for the intercross and the sample size [33,34]. For the partial correlation calculation of G×G×G_Nkx effects versus defect severity (Figure 3e), we controlled for variation in the observed incidence at two-locus genotypes to avoid bias by an underlying correlation.

3. Results

3.1. Nkx2-5^+/− Mice Mimic the Epidemiological Relationship between the Severity and Incidence of a Heart Defect

We phenotyped ~10,000 Nkx2-5^+/− newborns from intercrosses between the inbred strains A/J and C57BL/6N (A×B, N = 2999) and FVB/N and C57BL/6N (F×B; N = 6958). An ASD is the most common defect associated with human NKX2-5 mutation [6,35]. An ASD is also most common in the mouse crosses, followed by membranous and muscular ventricular septal defects (VSD) and AVSD (Figure 1a). The two anatomic types of VSD have different developmental bases but similar pathophysiological consequences. We note this because epidemiological studies typically combine membranous and muscular VSDs, making VSD appear more common than ASD [1]. Both intercrosses recapitulate the inverse relationship between the severity and incidence of a heart defect (Figure 1b). The attrition of more severely affected embryos does not explain the distribution of defects because Nkx2-5^+/− mice are born at the expected Mendelian ratio. In fact, we obtained pups with very severe heart defects, such as double outlet right ventricle and tricuspid atresia, but too few for quantitative genetic analysis. About 70% of Nkx2-5^+/− mice from either intercross have a normal heart.

3.2. The Effect of G×G_Nkx Interactions between a QTL That Modifies the Risk of a Defect and Nkx2-5 Correlates with the Severity of the Defect

The total, quantitative effect of genetic modifiers on risk could be a function of the number of genes, their effect sizes, or both. To distinguish these possibilities, we mapped the modifier QTLs in the F2 intercross of A×B and in the combined F2 and advanced intercross generations of F×B (Figure S1). Mapping in an F2 intercross has greater power to detect loci, whereas the combined analysis of an F2 and later generations has greater mapping resolution [20]. There were no rare alleles in either cross. The minor allele frequency of SNPs was 0.481 ± 0.012 in the offspring of the A×B cross and 0.446 ± 0.036 in the F×B cross (mean ± s.d.). We mapped QTLs for ASD, membranous and muscular VSD, and AVSD (Figure 2a, Table S2). Each set of QTLs represents G×G_Nkx interactions in developmental pathways leading to a defect.

QTLs that overlap between defects may share a gene involved in the development of multiple cardiac structures (Table S2). Overlapping A×B loci include ones on Chromosome 3 for ASD and AVSD, Chromosome 4 for ASD, membranous VSD and AVSD, and Chromosome 8 for membranous and muscular VSD. Overlapping F×B loci include one on Chromosome 6 for membranous and muscular VSD. An AVSD locus on chromosome 5, in which B carries the risk allele, overlapped between the A×B and F×B crosses.

The number of detected loci is similar between the four heart defects (Figure S2a). G×G_Nkx effect sizes, however, vary with the severity of a defect. G×G_Nkx interactions have a small effect on ASD risk and increasingly larger effects for VSDs and AVSD. The correlation holds whether the effect is calculated as the risk associated with a susceptibility allele or the phenotypic variance explained by a locus (Figure 2b,c). To account for undetected loci, we estimated the phenotypic variability, i.e., heritability, explained by all genotyped SNPs. The heritability is similar between defects (Figure S2b). Therefore, the correlation between G×G_Nkx effect size and severity is not due to a larger total genetic contribution to the risk of severe defects.

3.3. The Effect of G×G×G_Nkx Interactions between Two QTLs and Nkx2-5 Correlates with the Severity of a Defect

While each QTL was mapped according to its individual effect on risk, two loci may interact. Representative examples compare the observed to the expected incidences at two-locus genotypes under the null model of non-interacting loci [30]. Small differences between the observed and expected incidences of ASD and VSD indicate that each locus exerts mainly independent or additive effects (Figure 3a,b; Figure S3 shows the observed incidences between the nine genotypes at all locus pairs). In contrast, the incidence of AVSD at a two-locus genotype can deviate substantially from the null expectation, indicating a non-linear effect. Some differences are very large relative to the population or the unweighted average incidence of a defect across the nine genotypes (Figure 3c,d). The unweighted average incidence permits comparisons between defects of higher-order epistatic effects independently of genotype frequencies [30]. The population and unweighted average incidences are correlated in the inbred strain crosses because there are only two common alleles per locus. In natural populations, multiple and rare alleles obscure this correlation and make epistasis difficult to assess.

Similar to G×G_Nkx interactions, G×G×G_Nkx effect sizes increase with the severity of a defect (Figure 3e). We calculated the G×G×G_Nkx contribution to the observed incidence of a defect at each of the nine genotypes between every pair of significant and suggestive loci (or the next most significant locus in the two cases where there is just one suggestive locus). Statistically significant G×G×G_Nkx effect sizes vary between defects. AVSD effects are larger than ASD and VSD effects. The relationship is not secondary to an underlying correlation, such as between the observed incidence at a two-locus genotype and the G×G×G_Nkx effect size or the severity of a defect.

3.4. Genetic Coadaptation between Interacting G×G×G_Nkx Loci Suppresses the Risk of Cardiac Malformation

Patterns of G×G×G_Nkx effects suggest that positive selection promoted the coadaptation of interacting genes, as opposed to the elimination of intrinsically deleterious alleles during the inbreeding of mouse strains. Among significant G×G×G_Nkx effects, two-locus genotypes that are homozygous at both loci for alleles from the same inbred strain or “syn-homozygous”, e.g., BB at both loci, are more likely to lower risk. Conversely, anti-homozygous G×G×G_Nkx effects, e.g., homozygous FF at one locus and BB at the other, are more likely to raise the risk. Heterozygosity at either locus is equally likely to raise or lower risk (Figure 4a). Syn- and anti-homozygous G×G×G_Nkx interactions can have major effects on the risk of a severe defect. For example, one anti-homozygous, two-locus genotype accounts for 50% of all the AVSDs in the F×B intercross (Figure 3d). The genotypes at either locus are not intrinsically deleterious because Nkx2-5 wild-type mice, including littermates of affected mutants, do not develop AVSDs. Conversely, both syn-homozygous genotypes at two loci in the A×B intercross effectively suppress AVSD risk (Figure 3c). In general, no one- or two-locus genotype skews ASD or VSD risk as extremely (Figure S3). Protective and deleterious effects of syn- and anti-homozygosity, respectively, are associated with congruent effects at the other syn- and anti-homozygous genotypes in a pair of interacting loci (Figure 4b–e). Internally consistent G×G×G_Nkx effects suggest that selection for coadaptation was a recurrent event in the A/J, FVB/N, and C57BL/6N inbred strains.

4. Discussion

Normal cardiac development is crucial for survival to reproductive age, but CHD phenotypes exact varying costs on fitness. One has to wonder whether the fact that severe defects are rarer than mild ones has an evolutionary basis [36]. In humans, highly penetrant, mostly de novo genetic abnormalities account for about one-third of all CHD cases [7], but they are not necessarily associated with the most severe defects. In fact, ASDs accounted for more than 10% of CHD cases in two large studies that examined the role of de novo mutation by whole-exome sequencing [37,38]. We propose instead that the inverse relationship is a consequence of the genetic architectures of individual phenotypes. The present results offer novel perspectives on the nature of genetic risk in CHD and the origins of mutational robustness in cardiac development.

Genetic heterogeneity and rare alleles obscure the biological significance of epistasis in human diseases. In contrast, the power to detect genetic interactions is greater in inbred strain crosses. Through such crosses, we show that the effects of pairwise and higher-order epistasis correlate with the severity of a heart defect in Nkx2-5^+/− mice. Other features of the genetic architecture, such as the number of QTLs, do not correlate. The QTLs for mild ASDs exert mainly small, additive effects on risk, while the QTLs for severe AVSDs exert large, additive, and non-linear effects. The QTLs for moderately severe ventricular septal defects have intermediate effects. Strikingly non-linear G×G×G_Nkx effects can either suppress the risk of or account for a large fraction of AVSDs in the Nkx2-5^+/− population.

A complex genetic model plausibly explains the inverse relationship between the severity and incidence of a defect in humans. Whereas a modestly deleterious mutation may be sufficient to cause a mild defect, additional genetic interactions may be necessary to push a mutation carrier past the liability threshold for a severe defect, as observed in the Nkx2-5^+/− mouse model. Severe defects would thus be rarer because the multi-locus genotypes that give rise to deleterious interactions are improbable. Consistent with this model, a small study of non-syndromic AVSD found that among 34 persons who carried a mutation of one of six AVSD genes (NIPBL, CHD7, CEP152, BMPR1A, ZFPM2, MDM4), eight persons carried two or more mutations of the six genes [39]. Larger human studies are necessary to validate an oligogenic model in which non-linear effects may have an outsized role in severe CHD.

The model suggests a corollary: epistasis should matter less for strongly deleterious mutations that place individuals near the liability threshold for a defect. To exert the same phenotypic effect as a weak mutation and its contributing genetic interactions, a strong mutation must either obstruct one critical pathway or perturb several simultaneously. The dramatic effect of de novo mutations of histone-modifying genes, which regulate the expression of multiple developmental genes, is consistent with the latter possibility [40].

The protective and deleterious effects of syn- and anti-homozygosity in G×G×G_Nkx interactions are consistent with genetic coadaptation, a phenomenon first described by Dobzhansky in 1948. In crosses of Drosophila isolated from geographically distant locales, he noted “gene complexes” that conferred greater fitness when their alleles were from the same, rather than different populations [41,42]. Increased fitness emerged from the interaction between genotypes rather than their independent effects. Genetic coadaptation is probably widespread but under-recognized. Known examples relate to growth under different conditions in yeast [11], male fertility in Drosophila [12], body composition phenotypes in mice [43], and vitamin D receptor and skin color gene coadaptation in humans [44]. An elegant example examined the genetic regulation of fasting blood glucose levels, a diabetes-related trait. In a set of mouse crosses between chromosomal substitution strains, one or two different chromosomes from the A/J strain were introduced into the C57BL/6J background. Five interactions were discovered between pairs of different A/J chromosomes. Compared with the control C57BL/6J strain, mice that carried either one of the A/J chromosomes in an interacting pair (e.g., chromosome 5 or 6) had elevated glucose levels. Under an additive model, mice that carry both A/J chromosomes (e.g., chromosomes 5 and 6) should have higher glucose levels but actually had normal levels [45]. A/J and C57BL/6J appear to have arrived at different genetic solutions to maintain normal blood glucose levels.

Finally, the results suggest a rare example of selection for genetic robustness [18]. Whether robustness to mutation arises by selection or as an inherent consequence of an adapted trait has been a longstanding debate [17,46]. Depending upon the experimental design, investigators lacked either the ancestral population to assess the pre-robust state or the years required to detect selection for robustness in an evolving population [17]. Our application of inbred strains circumvents these challenges. Inbreeding strongly selects for fitness under genome-wide homozygosity. Unknown pressures appear to have selected strain-specific genotypes at interacting loci that enhance the robustness of cardiac developmental pathways to subsequent perturbation. The nkx2-5 mutation was probably not the original selection pressure because there are many potential causes of CHD. Nevertheless, outbreeding in the intercrosses produced individual variation in the robustness to the Nkx2-5 mutation. Modifier locus genotypes permit us to observe the analog of pre- and post-selected states in the same population. If humans similarly vary in the robustness of cardiac developmental pathways, the genetic basis of CHD is a function not only of deleterious variants but also of the degree of allelic compatibility between genes in the individual.