Detecting major genetic loci controlling phenotypic variability in experimental crosses

Traditional methods for detecting genes that affect complex diseases in humans or animal models, milk production in livestock, or other traits of interest, have asked whether variation in genotype produces a change in that trait’s average value. But focusing on differences in the mean ignores differences in variability about that mean. The robustness, or uniformity, of an individual’s character is not only of great practical importance in medical genetics and food production but is also of scientific and evolutionary interest (e.g., blood pressure in animal models of heart disease, litter size in pigs, flowering time in plants). We describe a method for detecting major genes controlling the phenotypic variance, referring to these as vQTL. Our method uses a double generalized linear model with linear predictors based on probabilities of line origin. We evaluate our method on simulated F₂ and collaborative cross data, and on a real F₂ intercross, demonstrating its accuracy and robustness to the presence of ordinary mean-controlling QTL. We also illustrate the connection between vQTL and QTL involved in epistasis, explaining how these concepts overlap. Our method can be applied to a wide range of commonly used experimental crosses and may be extended to genetic association more generally.QUANTITATIVE trait locus (QTL) analysis has traditionally focused on detection of major genes controlling the expected mean of a phenotype. But there is substantial evidence that not only the mean but also the variance, that is, the stochastic variability of the phenotype about its average value, may itself be under genetic control. The identification of such variance-controlling loci, which we call vQTL, can be helpful in a variety of contexts, including selection of livestock for uniformity, evaluating predictability of response to medical treatment, identification of key biomolecular stabilizers, and assessment of population resilience in ecology and evolution.One way of interpreting an increase in variability is as a decrease in stability. Waddington (1942) described the concept of canalization, whereby natural selection favors the relative constancy of some attributes, for example, well-formed organs and limbs, and thereby leads to the evolution of heritable architectures that buffer the impact of environmental or background genetic variation that would otherwise cause development to go astray. These architectures create virtual “canals” down which developmental programs flow. For a canalized phenotype, which modern usage expands to include nondevelopmental traits, the “zone of canalization” is the range of underlying liability over which potentially disruptive variation may be absorbed without serious consequence to the expressed trait value (Lynch and Walsh 1998). A well-studied example of a stabilizing architecture is that provided by heat-shock protein 90 (Hsp90), which buffers genetic and stochastic variation in the development of plants and flies (Rutherford and Lindquist 1998; Queitsch et al. 2002; Sangster et al. 2008).But in absorbing variation, such stabilizing architectures also hide it from view, and a sensitizing change in the stabilizer that shifts liability outside the zone of canalization can have a dramatic effect on the phenotype. Such shifts release the combined effects of previously “cryptic” genetic variation: now decanalized, the phenotype is more sensitive to internal (including genetic) and external environment, and as a result varies more greatly between individuals (Dworkin 2005; Hornstein and Shomron 2006). In this vein, decanalization has been proposed to explain why the genetic architectures of some diseases in human populations seem more amenable than others to genetic dissection through genome-wide association (Gibson and Goldstein 2007). Specifically, whereas some disease phenotypes in homogeneous populations may be heavily canalized and thereby harder to dissect, others may have been decanalized by modern living conditions (e.g., inflammatory diseases) or modern admixture, while yet others are simply too recent in evolutionary history for buffering networks to have evolved (e.g., response to HIV).Increased variability can also be adaptive. In natural populations disruptive selection favors diversity, with increased “capacitance” (Rice 2008) or “bet-hedging” (Beaumont et al. 2009) spreading risk over a variable fitness landscape. Feinberg and Irizarry (2010) recently proposed a heritable and selectable mechanism for this based on stochastic epigenetic variation. In controlled populations, variability can be increased through directional selection. For example, in a Drosophila selection experiment Clayton and Robertson (1957) reported increased bristle number variance, which is consistent with the idea that genotypes associated with higher environmental variance have a greater chance of being selected under directional selection (Hill and Zhang 2004). Moreover, genetic differences have been observed for phenotypic variability in body weight for chickens (Rowe et al. 2006) and snails (Ros et al. 2004) and litter size in rabbits (Ibanez-Escriche et al. 2008), sheep (Sancristobal-Gaudy et al. 1998), and pigs (Sorensen and Waagepetersen 2003).In natural populations with stabilizing selection we should expect to find alleles minimizing variance for fitness traits (Lande 1980; Houle 1992), whereas directional selection during domestication will favor alleles that increase variance. One may therefore expect to find vQTL in experimental crosses between wild and domestic animals (see Andersson 2001). Nonetheless, genetic buffering that leads to phenotypic robustness need not require an evolutionary explanation to be observed, nor to be useful in medicine and agriculture. Plainly, detecting vQTL and inferring how they arose are separate questions; here we concentrate on the first.

Healing of Euchromatic Chromosome Breaks by Efficient de novo Telomere Addition in Drosophila melanogaster

Many arthropod species are infected with maternally inherited endosymbionts that induce a shift in the sex ratio of their hosts by feminizing or killing males (cytoplasmic sex-ratio distorters, or SRDs). These endosymbionts can have profound impacts on evolutionary processes of their hosts. Here, I derive analytical expressions for the coalescent effective size N_e of populations that are infected with SRDs. Irrespective of the type of SRD, N_e for mitochondrial genes is given by the number of infected females. For nuclear genes, the effective population size generally decreases with increasing prevalence of the SRD and can be considerably lower than the actual size of the population. For example, with male-killing bacteria that have near perfect maternal transmission, N_e is reduced by a factor that is given to a good approximation by the proportion of uninfected individuals in the population. The formulae derived here also yield the effective size of populations infected with mutualistic endosymbionts or maternally inherited bacteria that induce cytoplasmic incompatibility, although in these cases, the reduction in N_e is expected to be less severe than for cytoplasmic SRDs.SIMPLE null models are essential in science. In population genetics, this role is filled by the Wright–Fisher model and its retrospective counterpart, the Kingman coalescent. Both of these models have proven to be immensely useful in spite of the fact that natural populations usually violate the assumptions made in these models. The reason for this is that often, the Wright–Fisher model can be rescaled so that it behaves in many important respects like a more complex population model. This rescaling is achieved through the concept of the effective population size, N_e. Roughly speaking, a complex population model is said to have a certain N_e if the haploid Wright–Fisher model with population size N_e experiences the same amount of random genetic drift as the complex model. Reflecting the different ways in which drift can be measured, N_e can be defined in different ways, e.g., as the inbreeding, the variance, or the coalescent effective population size. Different definitions often produce the same value for N_e, but may also yield drastically different numbers (Kimura and Crow 1963).The coalescent effective population size is defined through the factor by which time needs to be rescaled in a complex population model to produce the standard coalescent with time scale given by the population size N (Nordborg and Krone 2002). It has been argued that this is the most useful definition for N_e because “the coalescent essentially embodies all of the information that can be found in sampled genetic data” (Sjödin et al. 2005). More recently, Wakeley and Sargsyan (2009) have proposed two extensions of the coalescent effective population size in which they advocate including a mutation parameter in the definition and also allowing for a nonlinear relationship between N_e and N.One frequently encountered feature in natural populations that complicates population genetics is infection with maternally inherited endosymbionts. In particular, many arthropod species harbor a great number of phylogenetically diverse microorganisms that influence their hosts'' biology in different ways (Bourtzis and Miller 2003; Bourtzis and Miller 2006; Bourtzis and Miller 2009). Because of their maternal transmission, many of these microorganisms—for example, the bacteria Wolbachia pipientis and Cardinium hertigii—have evolved intricate manipulations of their hosts'' reproductive system that allows them to spread in a host population through exploitation of male hosts (reproductive parasitism, reviewed in Engelstädter and Hurst 2009a). Most manipulations involve a shift in the sex ratio of their hosts (both primary and at the population level), and the inducing endosymbionts are consequently referred to as cytoplasmic sex-ratio distorters (SRDs). In some species, genetic males develop into females if they are infected (“feminization”: Martin et al. 1973; Rigaud 1997; Bouchon et al. 1998; Hiroki et al. 2002). More commonly, infected males are killed by the endosymbionts early in their development (male killing: reviewed in Hurst et al. 2003). The adaptive advantage of this strategy is seen in an early fitness boost in the surviving females in a brood, for example, through reduced sibling competition or cannibalism of the dead brothers (Hurst 1991; Hurst and Majerus 1993; Jaenike et al. 2003). Some examples for species infected with male-killing or feminizing SRDs are given in Huigens and Stouthamer 2003).

TABLE 1

Empirical examples for cytoplasmic SRDs with parameter estimates

Efficient Transposition of Tol2 in the Mouse Germline

Insertional mutagenesis screens play an integral part in the annotating of functional data for all sequenced genes in the postgenomic era. Chemical mutagenesis screens are highly efficient but identifying the causative gene can be a laborious task. Other mutagenesis platforms, such as transposable elements, have been successfully applied for insertional mutagenesis screens in both the mouse and rat. However, relatively low transposition efficiency has hampered their use as a high-throughput forward genetic mutagenesis screen. Here we report the first evidence of germline activity in the mouse using a naturally active DNA transposon derived from the medaka fish called Tol2, as an alternative system for high-throughput forward genetic mutagenesis screening tool.THE Tol2 (transposable element of Oryzias latipes, number 2) element belongs to the hAT family (hobo of Drosophilia, Activator of maize and Tam3 of snapdragon) of transposons and was the first known autonomously active vertebrate type II transposable element (Koga et al. 1996; Kawakami et al. 1998). Unlike other DNA-type transposons like Sleeping Beauty (SB) (Ivics et al. 1997) or piggyBac (PB) (Fraser et al. 1996), Tol2 does not exhibit any known strong site specificity for integration nor does it exhibit any significant overexpression inhibition activity (Kawakami and Noda 2004; Balciunas et al. 2006) as seen in SB (Geurts et al. 2003). Recently, Tol2 was shown to effectively carry large DNA cargo of up to 10 kb in human and mouse cells without affecting its transposition efficiency (Balciunas et al. 2006). To date, Tol2 has also been demonstrated to transpose efficiently in zebrafish, frog, chicken, mouse cells, and human cells (Kawakami et al. 2000, 2004; Koga et al. 2003; Kawakami and Noda 2004; Balciunas et al. 2006; Hamlet et al. 2006; Sato et al. 2007).Germline mutagenesis using the SB transposon system has been demonstrated in both the mouse (Dupuy et al. 2001; Horie et al. 2001) and rat (Kitada et al. 2007; Lu et al. 2007). In addition, PB germline mutagenesis in mice has also been demonstrated (Ding et al. 2005; Wu et al. 2007). However, the relatively low germline transposition efficiency of both transposon systems reported so far has hampered their use in a high-throughput forward genetic mutagenesis screen (Keng et al. 2005; Kitada et al. 2007).

TABLE 1

Germline transposition frequency in various transposon systems

Retrograde Intraflagellar Transport Mutants Identify Complex A Proteins With Multiple Genetic Interactions in Chlamydomonas reinhardtii

The intraflagellar transport machinery is required for the assembly of cilia. It has been investigated by biochemical, genetic, and computational methods that have identified at least 21 proteins that assemble into two subcomplexes. It has been hypothesized that complex A is required for retrograde transport. Temperature-sensitive mutations in FLA15 and FLA17 show defects in retrograde intraflagellar transport (IFT) in Chlamydomonas. We show that IFT144 and IFT139, two complex A proteins, are encoded by FLA15 and FLA17, respectively. The fla15 allele is a missense mutation in a conserved cysteine and the fla17 allele is an in-frame deletion of three exons. The flagellar assembly defect of each mutant is rescued by the respective transgenes. In fla15 and fla17 mutants, bulges form in the distal one-third of the flagella at the permissive temperature and this phenotype is also rescued by the transgenes. These bulges contain the complex B component IFT74/72, but not α-tubulin or p28, a component of an inner dynein arm, which suggests specificity with respect to the proteins that accumulate in these bulges. IFT144 and IFT139 are likely to interact with each other and other proteins on the basis of three distinct genetic tests: (1) Double mutants display synthetic flagellar assembly defects at the permissive temperature, (2) heterozygous diploid strains exhibit second-site noncomplemention, and (3) transgenes confer two-copy suppression. Since these tests show different levels of phenotypic sensitivity, we propose they illustrate different gradations of gene interaction between complex A proteins themselves and with a complex B protein (IFT172).CILIA and flagella are microtubule-based organelles that are found on most mammalian cells. They provide motility to cells and participate in many sensory processes. Defects in or loss of cilia/flagella cause a variety of human diseases that include polycystic kidney disease, retinal degeneration, infertility, obesity, respiratory defects, left–right axis determination, and polydactyly (Fliegauf et al. 2007). Mouse mutants demonstrate that cilia are essential for viability, neural tube closure, and bone development (Eggenschwiler and Anderson 2007; Fliegauf et al. 2007). Cilia and flagella are also present in protists, algae, moss, and some fungi.The assembly and maintenance of cilia and flagella require intraflagellar transport (IFT) (Kozminski et al. 1995). IFT involves the movement of 100- to 200-nm-long protein particles from the basal body located in the cell body to the tip of the flagella using the heterotrimeric kinesin-2 (anterograde movement) (Kozminski et al. 1995) and movement back to the cell body (retrograde movement) using the cytoplasmic dynein complex (Pazour et al. 1999; Porter et al. 1999). IFT particles change their direction of movement as well as their size, speed, and frequency at the ends of the flagella as they switch from anterograde to retrograde movement (Iomini et al. 2001). Biochemical isolation of IFT particles reveals that they are composed of at least 16 proteins and that these particles can be dissociated into two complexes in vitro by changing the salt concentration (Cole et al. 1998; Piperno et al. 1998). Recent genetic and bioinformatics analysis adds at least 7 more proteins to the IFT particle (Follit et al. 2009) (Eggenschwiler and Anderson 2007).

TABLE 1

Proteins and gene names for the intraflagellar transport particles in Chlamydomonas, C. elegans, and mouse

Variation in Genomic Recombination Rates Among Heterogeneous Stock Mice

We used a large panel of pedigreed, genetically admixed house mice to study patterns of recombination rate variation in a leading mammalian model system. We found considerable inter-individual differences in genomic recombination rates and documented a significant heritable component to this variation. These findings point to clear variation in recombination rate among common laboratory strains, a result that carries important implications for genetic analysis in the house mouse.THE rate of recombination—the amount of crossing over per unit DNA—is a key parameter governing the fidelity of meiosis. Recombination rates that are too high or too low frequently give rise to aneuploid gametes or prematurely arrest the meiotic cell cycle (Hassold and Hunt 2001). As a consequence, recombination rates should experience strong selective pressures to lie within the range defined by the demands of meiosis (Coop and Przeworski 2007). Nonetheless, classical genetic studies in Drosophila (Chinnici 1971; Kidwell 1972; Brooks and Marks 1986), crickets (Shaw 1972), flour beetles (Dewees 1975), and lima beans (Allard 1963) have shown that considerable inter-individual variation for recombination rate is present within populations. Recent studies examining the transmission of haplotypes in human pedigrees have corroborated these findings (Broman et al. 1998; Kong et al. 2002; Coop et al. 2008).Here, we use a large panel of heterogeneous stock (HS) mice to study variation in genomic recombination rates in a genetic model system. These mice are genetically admixed, derived from an initial generation of pseudorandom mating among eight common inbred laboratory strains (DBA/2J, C3H/HeJ, AKR/J, A/J, BALB/cJ, CBA/J, C57BL/6J, and LP/J), followed by >50 generations of pseudorandom mating in subsequent hybrid cohorts (Mott et al. 2000; Demarest et al. 2001). The familial relationships among animals in recent generations were tracked to organize the mice into pedigrees. In total, this HS panel includes ∼2300 animals comprising 85 families, 8 of which span multiple generations. The remainder consists of nuclear families (sibships) that range from 1 to 34 sibs, with an average of 9.6 sibs (Valdar et al. 2006) (Mott et al. 2000; Demarest et al. 2001; Shifman et al. 2006).

TABLE 1

Heterogeneous stock mouse pedigrees

On the Classification of Epistatic Interactions

Modern genomewide association studies are characterized by the problem of “missing heritability.” Epistasis, or genetic interaction, has been suggested as a possible explanation for the relatively small contribution of single significant associations to the fraction of variance explained. Of particular concern to investigators of genetic interactions is how to best represent and define epistasis. Previous studies have found that the use of different quantitative definitions for genetic interaction can lead to different conclusions when constructing genetic interaction networks and when addressing evolutionary questions. We suggest that instead, multiple representations of epistasis, or epistatic “subtypes,” may be valid within a given system. Selecting among these epistatic subtypes may provide additional insight into the biological and functional relationships among pairs of genes. In this study, we propose maximum-likelihood and model selection methods in a hypothesis-testing framework to choose epistatic subtypes that best represent functional relationships for pairs of genes on the basis of fitness data from both single and double mutants in haploid systems. We gauge the performance of our method with extensive simulations under various interaction scenarios. Our approach performs reasonably well in detecting the most likely epistatic subtype for pairs of genes, as well as in reducing bias when estimating the epistatic parameter (ɛ). We apply our approach to two available data sets from yeast (Saccharomyces cerevisiae) and demonstrate through overlap of our identified epistatic pairs with experimentally verified interactions and functional links that our results are likely of biological significance in understanding interaction mechanisms. We anticipate that our method will improve detection of epistatic interactions and will help to unravel the mysteries of complex biological systems.UNDERSTANDING the nature of genetic interactions is crucial to obtaining a more complete picture of complex biological systems and their evolution. The discovery of genetic interactions has been the goal of many researchers studying a number of model systems, including but not limited to Saccharomyces cerevisiae, Caenorhabditis elegans, and Escherichia coli (You and Yin 2002; Burch et al. 2003; Burch and Chao 2004; Tong et al. 2004; Drees et al. 2005; Sanjuán et al. 2005; Segre et al. 2005; Pan et al. 2006; Zhong and Sternberg 2006; Jasnos and Korona 2007; St. Onge et al. 2007; Decourty et al. 2008). Recently, high-throughput experimental approaches, such as epistatic mini-array profiles (E-MAPs) and genetic interaction analysis technology for E. coli (GIANT-coli), have enabled the study of epistasis on a large scale (Schuldiner et al. 2005, 2006; Collins et al. 2006, 2007; Typas et al. 2008). However, it remains unclear whether the computational and statistical methods currently in use to identify these interactions are indeed the most appropriate.The study of genetic interaction, or “epistasis,” has had a long and somewhat convoluted history. Bateson (1909) first used the term epistasis to describe the ability of a gene at one locus to “mask” the mutational influence of a gene at another locus (Cordell 2002). The term “epistacy” was later coined by Fisher (1918) to denote the statistical deviation of multilocus genotype values from an additive linear model for the value of a phenotype (Phillips 1998, 2008).These origins are the basis for the two main current interpretations of epistasis. The first, as introduced by Bateson (1909), is the “biological,” “physiological,” or “compositional” form of epistasis, concerned with the influence of an individual''s genetic background on an allele''s effect on phenotype (Cheverud and Routman 1995; Phillips 1998, 2008; Cordell 2002; Moore and Williams 2005). The second interpretation, attributed to Fisher, is “statistical” epistasis, which in its linear regression framework places the phenomenon of epistasis in the context of a population (Wagner et al. 1998; Wade et al. 2001; Wilke and Adami 2001; Moore and Williams 2005; Phillips 2008). Each of these approaches is equally valid in studying genetic interactions; however, confusion still exists about how to best reconcile the methods and results of the two (Phillips 1998, 2008; Cordell 2002; Moore and Williams 2005; Liberman and Feldman 2006; Aylor and Zeng 2008).Aside from the distinction between the statistical and the physiological definitions of epistasis, inconsistencies exist when studying solely physiological epistasis. For categorical traits, physiological epistasis is clear as a “masking” effect. When noncategorical or numerical traits are measured, epistasis is defined as the deviation of the phenotype of the multiple mutant from that expected under independence of the underlying genes.The “expectation” of the phenotype under independence, that is, in the absence of epistasis, is not defined consistently between studies. For clarity, consider epistasis between pairs of genes and, without loss of generality, consider fitness as the phenotype. The first commonly used definition of independence, originating from additivity, defines the effect of two independent mutations to be equal to the sum of the individual mutational effects. A second, motivated by the use of fitness as a phenotype, defines the effect of the two mutations as the product of the individual effects (Elena and Lenski 1997; Desai et al. 2007; Phillips 2008). A third definition of independence has been referred to as “minimum,” where alleles at two loci are independent if the double mutant has the same fitness as the less-fit single mutant. Mani et al. (2008) claim that this has been used when identifying pairwise epistasis by searching for synthetic lethal double mutants (Tong et al. 2001, 2004; Pan et al. 2004, 2006; Davierwala et al. 2005). A fourth is the “Log” definition presented by Mani et al. (2008) and Sanjuan and Elena (2006). The less-frequently used “scaled ɛ” (Segre et al. 2005) measure of epistasis takes the multiplicative definition of independence with a scaling factor.These different definitions of independence are partly due to distinct measurement “scales.” For some traits, a multiplicative definition of independence may be necessary to identify epistasis between two genes, whereas for other traits, additivity may be appropriate (Falconer and Mackay 1995; Wade et al. 2001; Mani et al. 2008; Phillips 2008). An interaction found under one independence definition may not necessarily be found under another, leading to different biological conclusions (Mani et al. 2008).Mani et al. (2008) suggest that there may be an “ideal” definition of independence for all gene pairs for identifying functional relationships. However, it is plausible that different representations of independence for two genes may reflect different biological properties of the relationship (Kupper and Hogan 1978; Rothman et al. 1980). “Two categories of general interest [the additive and multiplicative definitions, respectively] are those in which etiologic factors act interchangeably in the same step in a multistep process, or alternatively act at different steps in the process” (Rothman et al. 1980, p. 468). In some cases, the discovery of epistasis may merely be an artifact of using an incorrect null model (Kupper and Hogan 1978). It may be necessary to represent “independence” differently, resulting in different statistical measures of interactions, for different pairs of genes depending on their functions.Previous studies have suggested that different pairs of loci may have different modes of interaction and have attempted to subclassify genetic interactions into regulatory hierarchies and mutually exclusive “interaction subtypes” to elucidate underlying biological properties (Avery and Wasserman 1992; Drees et al. 2005; St. Onge et al. 2007). We suggest that epistatic relationships can be divided into several subtypes, or forms, corresponding to the aforementioned definitions of independence. As a particular gene pair may deviate from independence according to several criteria, we do not claim that these subtypes are necessarily mutually exclusive. We attempt to select the most likely epistatic subtype that is the best statistical representation of the relationship between two genes. To further subclassify interactions, epistasis among deleterious mutations can take one of two commonly used forms: positive (equivalently alleviating, antagonistic, or buffering) epistasis, where the phenotype of the double mutant is less severe than expected under independence, and negative (equivalently aggravating, synergistic, or synthetic), where the phenotype is more severe than expected (Segre et al. 2005; Collins et al. 2006; Desai et al. 2007; Mani et al. 2008).Another objective of such distinctions is to reduce the bias of the estimator of the epistatic parameter (ɛ), which measures the extent and direction of epistasis for a given gene pair. Mani et al. (2008), assuming that the overall distribution of ɛ should be centered around 0, find that inaccurately choosing a definition of independence can result in increased bias when estimating ɛ. For example, using the minimum definition results in the most severe bias when single mutants have moderate fitness effects, and the additive definition results in the largest positive bias when at least one gene has an extreme fitness defect (Mani et al. 2008). Therefore, it is important to select an optimal estimator for ɛ for each pair of genes from among the subtypes of epistatic interactions.Epistasis may be important to consider in genomic association studies, as a gene with a weak main effect may be identified only through its interaction with another gene or other genes (Frankel and Schork 1996; Culverhouse et al. 2002; Moore 2003; Cordell 2009; Moore and Williams 2009). Epistasis has also been studied extensively in the context of the evolution of sex and recombination. The mutational deterministic hypothesis proposes that the evolution of sex and recombination would be favored by negative epistatic interactions (Feldman et al. 1980; Kondrashov 1994); many other studies have also studied the importance of the form of epistasis (Elena and Lenski 1997; Otto and Feldman 1997; Burch and Chao 2004; Keightley and Otto 2006; Desai et al. 2007; MacCarthy and Bergman 2007). Indeed, according to Mani et al. (2008, p. 3466), “the choice of definition [of epistasis] alters conclusions relevant to the adaptive value of sex and recombination.”Given fitness data from single and double mutants in haploid organisms, we implement a likelihood method to determine the subtype that is the best statistical representation of the epistatic interaction for pairs of genes. We use maximum-likelihood estimation and the Bayesian information criteria (BIC) (Schwarz 1978) with a likelihood-ratio test to select the most appropriate null or epistatic model for each putative interaction. We conduct extensive simulations to gauge the performance of our method and demonstrate that it performs reasonably well under various interaction scenarios. We apply our method to two data sets with fitness measurements obtained from yeast (Jasnos and Korona 2007; St. Onge et al. 2007), whose authors assume only multiplicative epistasis for all interactions. By examining functional links and experimentally validated interactions among epistatic pairs, we demonstrate that our results are biologically meaningful. Studying a random selection of genes, we find that minimum epistasis is more prevalent than both additive and multiplicative epistasis and that the overall distribution of ɛ is not significantly different from zero (as Jasnos and Korona 2007 suggest). For genes in a particular pathway, we advise selecting among fewer epistatic subtypes. We believe that our method of epistatic subtype classification will aid in understanding genetic interactions and their properties.

St. Onge et al. (2007) data set:

St. Onge et al. (2007) examined 26 nonessential genes known to confer resistance to MMS, constructed double-deletion strains for 323 double-mutant strains (all but two of the total possible pairs), and assumed the multiplicative form of epistasis for all interactions (see Methods: Analysis of experimental data). Following these authors, we focus on single- and double-mutant fitnesses measured in the presence of MMS. (For results in the absence of MMS, see File S1 and File S1_2.)Using the resampling method described in Analysis of experimental data and File S1, 222 gene pairs pass the cutoff of having epistasis inferred in at least 900 of 1000 replicates. This does not include 5 synthetic lethal gene pairs. Hypothesis testing and a multiple-testing procedure (for 222 simultaneous hypotheses) are necessary to determine the final epistatic pairs.To select one among the three multiple-testing procedures, we follow St. Onge et al. (2007) and examine gene pairs that share specific functional links (see Analysis of experimental data). The Bonferroni method is likely too conservative, yielding only 25 significantly epistatic pairs with only one functional link among them; alternatively, the pFDR procedure appears to be too lenient in rejecting independence for all 222 pairs. Therefore, we use the FDR procedure (although the number of functional links is not significant) and detect 193 epistatic pairs, of which 5 (2.6%) are synthetic lethals, 19 (9.8%) have additive epistasis, 33 (17.1%) have multiplicative epistasis, and 136 (70.5%) have minimum epistasis (File S1_1). We find 29 gene pairs with positive (alleviating) epistasis and 159 gene pairs with negative (aggravating) epistasis.

TABLE 2

Summary of gene pairs with the indicated epistatic subtypes, inferred using the FDR procedure with the BIC method that considers all three epistatic subtypes and their corresponding null models

Surprising Fitness Consequences of GC-Biased Gene Conversion: I. Mutation Load and Inbreeding Depression

GC-biased gene conversion (gBGC) is a recombination-associated process mimicking selection in favor of G and C alleles. It is increasingly recognized as a widespread force in shaping the genomic nucleotide landscape. In recombination hotspots, gBGC can lead to bursts of fixation of GC nucleotides and to accelerated nucleotide substitution rates. It was recently shown that these episodes of strong gBGC could give spurious signatures of adaptation and/or relaxed selection. There is also evidence that gBGC could drive the fixation of deleterious amino acid mutations in some primate genes. This raises the question of the potential fitness effects of gBGC. While gBGC has been metaphorically termed the “Achilles'' heel” of our genome, we do not know whether interference between gBGC and selection merely has practical consequences for the analysis of sequence data or whether it has broader fundamental implications for individuals and populations. I developed a population genetics model to predict the consequences of gBGC on the mutation load and inbreeding depression. I also used estimates available for humans to quantitatively evaluate the fitness impact of gBGC. Surprising features emerged from this model: (i) Contrary to classical mutation load models, gBGC generates a fixation load independent of population size and could contribute to a significant part of the load; (ii) gBGC can maintain recessive deleterious mutations for a long time at intermediate frequency, in a similar way to overdominance, and these mutations generate high inbreeding depression, even if they are slightly deleterious; (iii) since mating systems affect both the selection efficacy and gBGC intensity, gBGC challenges classical predictions concerning the interaction between mating systems and deleterious mutations, and gBGC could constitute an additional cost of outcrossing; and (iv) if mutations are biased toward A and T alleles, very low gBGC levels can reduce the load. A robust prediction is that the gBGC level minimizing the load depends only on the mutational bias and population size. These surprising results suggest that gBGC may have nonnegligible fitness consequences and could play a significant role in the evolution of genetic systems. They also shed light on the evolution of gBGC itself.GC-BIASED gene conversion (gBGC) is increasingly recognized as a widespread force in shaping genome evolution. In different species, gene conversion occurring during double-strand break recombination repair is thought to be biased toward G and C alleles. In heterozygotes, GC alleles undergo a kind of molecular meiotic drive that mimics selection (reviewed in Marais 2003). This process can rapidly increase the GC content, especially around recombination hotspots (Spencer et al. 2006), and, more broadly, can affect genome-wide nucleotide landscapes (Duret and Galtier 2009a). For instance, it is thought to play a role in shaping isochore structure evolution in mammals (Galtier et al. 2001; Meunier and Duret 2004; Duret et al. 2006) and birds (Webster et al. 2006). Direct experimental evidence of gBGC mainly comes from studies in yeast (Birdsell 2002; Mancera et al. 2008; but see Marsolier-Kergoat and Yeramian 2009) and humans (Brown and Jiricny 1987). However, associations between recombination and the nucleotide landscape and frequency spectra biased toward GC alleles provide indirect evidence in very diverse organisms (OrganismsDirect evidenceIndirect evidenceAchille''s heel evidenceReferencesYeastMeiotic segregation biasMancera et al. (2008)Mitotic and mitotic heteromismatch correction biasCorrelation between GC and recombinationBirdsell (2002)MammalsMitotic heteromismatch correction biasBrown and Jiricny (1987)Correlation between GC*/GC and recombinationDuret and Arndt (2008); Meunier and Duret (2004)Biased frequency spectrum toward GC allelesGaltier et al. (2001); Spencer et al. (2006)GC bias associated with high d_N/d_S near recombination hotspotBerglund et al. (2009; Galtier et al. (2009)BirdsCorrelation between GC and recombinationInternational Chicken Genome Sequencing Consortium (2004)TurtlesCorrelation between GC and chromosome sizeKuraku et al. (2006)DrosophilaCorrelation between GC and recombinationMarais et al. (2003)Biased frequency spectrum toward GC allelesGaltier et al. (2006)NematodesCorrelation between GC and recombinationMarais et al. (2001)GrassesCorrelation between GC and outcrossing/selfingGlémin et al. (2006)Correlation between GC* and recombination and outcrossing/selfingOutcrossing increases d_N/d_S for genes with high GC*Haudry et al. (2008)Green algaeCorrelation between GC and recombinationJancek et al. (2008)ParameciumCorrelation between GC and chromosome sizeDuret et al. (2008)Open in a separate windowThe impact of gBGC on noncoding sequences and synonymous sites has been studied in depth, especially because of confounding effects with selection on codon usage (Marais et al. 2001). More recently, Galtier and Duret (2007) pointed out that gBGC may also interfere with selection when affecting functional sequences. They argued that gBGC could leave spurious signatures of adaptive selection and proposed to extend the null hypothesis of molecular evolution. Indeed, gBGC can lead to a ratio of nonsynonymous (d_N) over synonymous (d_S) substitutions above one (Berglund et al. 2009; Galtier et al. 2009), i.e., a typical signature of positive selection (Nielsen 2005). This hypothesis has been widely debated for human-accelerated regions (HARs). These regions are extremely conserved across mammals but show evidence of accelerated evolution along the human lineage, which has been interpreted as evidence of positive selection (Pollard et al. 2006a,b; Prabhakar et al. 2006, 2008). On the contrary, other authors argued that patterns observed in HARs, such as the AT → GC substitution bias, the absence of a selective sweep signature, or the propensity to occur within or close to recombination hotspots, are more likely explained by gBGC rather than positive selection (Galtier and Duret 2007; Berglund et al. 2009; Duret and Galtier 2009b; but see also Pollard et al. 2006a who also suggested that gBGC might play a role in HARs evolution). It is thus crucial to take gBGC into account when interpreting genomic data.Moreover, Galtier and Duret (2007) initially suggested that gBGC hotspots could contribute to the fixation of slightly deleterious AT → GC mutations and could represent the Achilles'' heel of our genome. This hypothesis was reinforced later in primates, with evidence of gBGC-driven fixation of deleterious mutations in proteins (Galtier et al. 2009). A similar result was also found in some grass species, whose genomes are also supposed to be affected by gBGC (Glémin et al. 2006). Haudry et al. (2008) compared two outcrossing and two selfing grass species and showed that GC-biased genes exhibit higher d_N/d_S ratio in outcrossing than in selfing lineages. The reverse pattern would be expected under pure selective models because of the reduced selection efficacy in selfers (Charlesworth 1992; Glémin 2007). This pattern is in agreement with a genomic Achilles'' heel associated with outcrossing, while gBGC is inefficient in selfing species because they are mainly homozygous.Twenty years ago, Bengtsson (1990) already pointed out that biased conversion can generally affect the mutation load. The mutation load is the reduction in the mean fitness of a population due to mutation accumulation, which could lead to population extinction if it is too high (Lynch et al. 1995). At this time, Bengtsson concluded that “it is impossible to know if biased conversion plays a major role in determining the magnitude of the mutation load in organisms such as ourselves, but the possibility must be considered and further investigated (Bengtsson 1990, p. 186).” Now, one can propose gBGC could be such a widespread biased conversion process. It thus appears timely to thoroughly investigate the fitness consequences of gBGC through its potential effects on the dynamics of deleterious mutations. The fitness consequences of gBGC were also pointed out as a major future issue to be addressed by Duret and Galtier (2009a). In addition to the load, deleterious mutations have many other evolutionary consequences (for review see Charlesworth and Charlesworth 1998). They are thought to be the main determinant of inbreeding depression, i.e., the reduction in fitness of inbred individuals compared to outbred ones. They also play a key role in the evolution of genetic systems (sexual reproduction and recombination, inbreeding avoidance mechanisms, ploidy cycles), of senescence, or in the degeneration of nonrecombining regions, such as Y chromosomes. So far, we know little, if anything, about how gBGC might affect these processes.In his seminal work, Bengtsson (1990) did not address several important points. First, he did not include genetic drift in his model. Nearly neutral mutations, for which drift and selection are of similar intensities, are the most damaging ones because they can drift to fixation, unlike strongly deleterious mutations that are maintained at low frequency (Crow 1993; Lande 1994, 1998). While gBGC intensities are rather weak (Birdsell 2002; Spencer et al. 2006), they could markedly affect the fate of nearly neutral mutations (see also Galtier et al. 2009). Second, Bengtsson did not study the effect of gene conversion on inbreeding depression, while he showed that recessive mutations, mostly involved in inbreeding depression, are the most affected by gene conversion. Third, he did not envisage systematic GC bias with its opposite effects on A/T and G/C deleterious alleles. Fourth, while he noted that selfing affects both the efficacy of selection and that of conversion, he did not fully investigate the effect of mating systems. On one hand, selfing is efficient in purging strongly deleterious mutations causing inbreeding depression. However, since selfing is expected to increase drift, weakly deleterious mutations can fix in selfing species, contributing to the so-called “drift load” (Charlesworth 1992; Glémin 2007). Self-fertilizing populations are thus expected to exhibit low inbreeding depression and high drift load. On the other hand, gBGC, and thus its cost, vanishes as the selfing rate and homozygosity increase (Marais et al. 2004). gBGC could thus challenge classical views on mating systems and it was even speculated that gBGC could affect their evolution (Haudry et al. 2008).Here I present a population genetics model that includes mutation, selection, drift, and gBGC, which extends previous studies (Gutz and Leslie 1976; Lamb and Helmi 1982; Nagylaki 1983a,b; Bengtsson 1990). I specifically examine how gBGC can affect inbreeding depression and the mutation load. I also focus on the effect of mating system, which is especially interesting with regard to the interaction between biased conversion and selection. Finally, I discuss how these results could give insight into how gBGC evolved.

Impacts of gBGC on inbreeding depression:

Inbreeding depression is defined as the reduction in fitness of selfed (and more generally inbred) individuals compared to outcrossed individuals,(15)where and are the mean fitness of outcrosses and selfcrosses, respectively (Charlesworth and Charlesworth 1987; Charlesworth and Willis 2009). The approximation is very good in most conditions, because under weak (s ≪ 1) and strong selection (x ≪ 1) (see Glémin et al. 2003). Similar to the load, considering both sites for which either S or W alleles are deleterious, in proportion q and 1 – q, respectively, we get(16)

gBGC and the genetic basis of inbreeding depression in panmictic populations:

In infinite panmictic populations without gBGC, inbreeding depression depends only on mutation rates and dominance levels. Partially recessive mutations () contribute only to inbreeding depression, and the more recessive they are, the higher the inbreeding depression (Charlesworth and Charlesworth 1987). In finite populations, deterministic results hold for strongly deleterious mutations (s ≫ 1/N_e), which contribute mostly to inbreeding depression. Contrary to the load, weakly deleterious mutations (∼s ≤ 1/N_e) contribute little to inbreeding depression (Figure 4, a and c, and see Bataillon and Kirkpatrick 2000).Open in a separate windowFigure 4.—Inbreeding depression (×10⁶) as a function of s without (a and c) or with (b and d) gBGC (b = 0.0002). (a and b) h = 0.2: thick lines, N = 5000; thin lines, N = 10,000; dashed lines, N = 50,000; dotted lines, N = 100,000. (c and d) N = 10,000: thick lines, h = 0.4; thin lines, h = 0.2; dashed lines, h = 0.1; dotted lines, h = 0.05. u = 10⁻⁶, λ = 2.Like the load, gBGC affects both the magnitude and the structure of inbreeding depression. In infinite populations, and more generally for strongly deleterious alleles (N_es ≫ 1), replacing x by x_eq given by Equations 4 in Equations 15 and 16 leads to(17a)(17b)(17c)The effect of gBGC on inbreeding depression is not monotonic. Like the load, gBGC increases inbreeding depression if b > hs(1 − 2q/(q + λ − qλ)). However, contrary to the load, a strong gBGC decreases inbreeding depression, which tends to 0 as b increases, while the load tends to qs (Equation 10c). An analysis of Equation 17b shows that mutations that maximize inbreeding depression are those that also maximize the load, i.e., S deleterious mutations with s ≈ 2b.In finite populations, inbreeding depression must be integrated over the Φ distribution, which leads to(18)(see also Glémin et al. 2003). While it is not possible to get an analytical expression of (18), numerical computations (see appendix b) show that S deleterious mutations with s ≈ 2b also maximize inbreeding depression in finite populations (Figure 4). More broadly, inbreeding depression is maximal under the overdominant-like selection regime (gray area in Figure 2). Once again, even low to moderate gBGC markedly affects the genetic structure of inbreeding depression. First, mutations of intermediate effects contribute the most to inbreeding depression, i.e., up to one order of magnitude higher than strongly deleterious mutations (compare Figure 4a with 4b). Second, even nearly additive mutations can have a substantial effect (compare Figure 4c with 4d).Since little is known about the distribution of dominance coefficients, especially the dominance of mildly deleterious mutations (of the order of b), it is difficult to quantitatively predict the full impact of gBGC on inbreeding depression. We can conclude that, on average, gBGC should increase inbreeding depression. However, further insight into mutational parameters is crucial to assess the quantitative impact of gBGC.

Joint effect of gBGC and mating system on the load and inbreeding depression:

Selfing, or more generally inbreeding, slightly reduces the segregating load through the purging of recessive mutations (Ohta and Cockerham 1974), but can substantially increase the fixation load because of the effective population size reduction under inbreeding: (see above and Pollak 1987; Nordborg 1997; Glémin 2007). In numerical examples, I assumed that α decreases with F according to the background selection model (Charlesworth et al. 1993; Nordborg et al. 1996), as in Glémin (2007). With gBGC, selfing thus has two opposite effects on the fixation load. Selfing increases the drift load sensu stricto but decreases the fixation load due to gBGC. A surprising consequence is that the load can be higher in outcrossing than in selfing populations (Figure 5). Quantitatively this is also expected, even with a gBGC hotspot affecting just 3% of the genome (Figure 5 and Open in a separate windowFigure 5.—Effective population size (a and b) and the load (×10⁶) (c–f) as a function of F for different gBGC intensities (thick lines, b = 0; thin lines, b = 0.0001; dashed lines, b = 0.0002; dotted lines, b = 0.0005). The effective population size depends on F under the background selection (BS) model (Charlesworth et al. 1993), using Equations 16 and 17 in Glémin (2007): , where U is the genomic deleterious mutation rate, R is the genomic recombination rate, s_d is the mean selection coefficient against strongly deleterious mutations, and h_d is their dominance coefficient. N = 10,000, U = 0.2, h_d = 0.1, and s_d = 0.05. (a, c, and e) R = 5, “weak” BS; (b, d, and f) R = 0.5, “strong” BS. (c and d) Load averaged over half GC and half AT deleterious alleles, with a bias in favor of AT alleles. (e and f) Load averaged over 10% of GC deleterious alleles and 90% of AT deleterious alleles with a bias in favor of AT alleles; see Figure 3. h = 0.5, u = 10⁻⁶, and λ = 2.Generally, the effect of selfing is simpler for inbreeding depression. Purging, N_e reduction, and suppression of gBGC contribute to decreasing inbreeding depression in selfing populations (Figure 6a). However, there are special cases in which maximum inbreeding depression is reached for intermediate selfing rates (Figure 6b). In such cases, in outcrossing populations, gBGC is strong enough to sweep polymorphism out and reduce inbreeding depression (b > s, regime 1 in Figure 2). As the selfing rate increases, gBGC declines, and the selection dynamics become overdominant-like (regime 2, Figure 2), thus maximizing inbreeding depression. For high selfing rates, gBGC vanishes (regime 3 in Figure 2) and deleterious alleles are either purged or fixed if there is substantial drift. This is similar to the effect of selfing on inbreeding depression caused by asymmetrical overdominance, where inbreeding depression also peaks for intermediate selfing rates (Ziehe and Roberds 1989; Charlesworth and Charlesworth 1990). In the present case, the range of parameters leading to this peculiar behavior is narrow because the overdominant-like region depends on the selfing rates and can vanish either for low or for high selfing rates (Figure 2).Open in a separate windowFigure 6.—Inbreeding depression (×10⁶) as a function of F for different gBGC intensities (thick lines, b = 0; thin lines, b = 0.0001; dashed lines, b = 0.0002; dotted lines, b = 0.0005). Inbreeding depression is averaged over half GC and half AT deleterious alleles. The effective population size depends on F as in Figure 5 (same parameters). (a) s = 0.002; (b) s = 0.0005; (c) s = 0.0002. h = 0.2, u = 10⁻⁶, and λ = 2.

Minimum load and the evolution of gBGC and recombination landscapes:

Although gBGC may have deleterious fitness consequences, it is surprising that it evolved in many taxa (Duret and Galtier 2009a). Birdsell (2002) initially suggested that gBGC may have evolved as a response to mutational bias toward AT (λ > 1, here). Indeed, I show that a minimum load is reached for weak gBGC (b ≈ ln(λ)/4N, Equation 14). This result is very general whatever the distribution of fitness effects of mutations (appendix d). However, the range of optimal gBGC is narrow, and gBGC increases the load as far as b > ln(λ)/2N (appendix c). In humans, using N = 10,000 and λ = 2, gBGC levels that minimize the load are ∼1.17 × 10⁻⁵, i.e., one order of magnitude lower than the average bias observed in recombination hotspots (Myers et al. 2005). However, selection on conversion modifiers will not necessarily minimize the load because of gametic disequilibrium generated between modifiers and fitness loci (Bengtsson and Uyenoyama 1990). Selection for limitation of somatic AT-biased mutations could also have selected for GC-biased mismatch repair machinery (Brown and Jiricny 1987). If the bias level that would be selected for somatic reasons is >ln(λ)/2N, a side effect would be the generation of a substantial load at the population level. Finally, it is interesting to note that when synonymous codon positions are under selection for translation accuracy, optimal gBGC levels can be higher than gBGC levels that minimize the protein load, especially when most optimal codons end in G or C ().Conversely, gBGC could also affect the evolution of recombination landscapes, which could evolve to reduce the gBGC load. Surprisingly, for a given recombination/conversion level, the hotspot distribution does not appear to be optimal (Nishant and Rao 2005), one can speculate that the hotspot localization outside genes could be a response to avoid the deleterious effects of gBGC.Up to now, these verbal arguments have not been assessed theoretically (but see Bengtsson and Uyenoyama 1990 for a different kind of conversion bias). Population genetics models are necessary to test these hypotheses concerning the evolution of gBGC and recombination landscapes and to pinpoint the key parameters that might govern their evolution.

gBGC and the evolution of mating systems:

Deleterious mutations also play a crucial role in the evolution of mating systems. They are the main source of inbreeding depression, which balances the automatic advantage of selfing. The drift load is also thought to contribute to the extinction of selfing species. Since they are mainly homozygous, selfing species are mostly free from gBGC and its deleterious impacts. I discuss below how this might affect the evolution of mating systems.

Inbreeding depression and the shift in mating systems:

Inbreeding depression plays a key role in the evolution of mating systems (Charlesworth and Charlesworth 1987; Charlesworth 2006b). Since it balances the automatic advantage of selfing, high inbreeding depression favors outcrossing, while selfing can evolve when it is low. Moreover, selfing helps to purge strongly deleterious mutations, thus decreasing inbreeding depression. This positive feedback reinforces the disruptive selection on the selfing rate and prevents the transition from selfing to outcrossing (Lande and Schemske 1985).Theoretical results suggest that, in most conditions, gBGC would reinforce inbreeding depression in outcrossing populations (Figure 6), which would prevent the evolution of selfing. In reverse, if selfing is initially selected for, recurrent selfing would reduce the load through both purging and avoidance of gBGC. Under this scenario, gBGC would reinforce disruptive selection on mating systems. However, under some conditions (see Figure 6), inbreeding depression peaks at intermediate selfing rates, as observed for asymmetrical overdominance (Ziehe and Roberds 1989; Charlesworth and Charlesworth 1990). In theory, this could prevent the shift toward complete selfing and maintain stable mixed mating systems (Charlesworth and Charlesworth 1990; Uyenoyama and Waller 1991). However, this pattern is observed under restrictive conditions and it is very unlikely on the whole-genome scale. Dominance patterns are crucial for predicting inbreeding depression, especially with gBGC. Contrary to the load, it is thus difficult to evaluate the quantitative impact of gBGC on inbreeding depression. However, increased inbreeding depression in outcrossing species subject to gBGC seems to be the most likely scenario.

gBGC and the long-term evolution of mating systems:

In the long term, the gBGC-induced load also challenges the “dead-end hypothesis,” which posits that, because of the reduction of selection efficacy, self-fertilizing species would accumulate weakly deleterious mutations in the long term, eventually leading to extinction (Takebayashi and Morrell 2001). Because of gBGC, not drift, outcrossing species could also accumulate a load of weakly deleterious mutations (Figure 7), and they could suffer from a higher load than highly self-fertilizing species (Haudry et al. (2008) found that in two outcrossing grass species, but not in two self-fertilizing ones, the d_N/d_S ratio is significantly higher for genes exhibiting GC enrichment. They speculated that substitutions in these genes might contribute to increasing the load in these two outcrossing grass species. Such results are still very sparse. In plants, evidence of strong gBGC is mainly restricted to grasses (but see Wright et al. 2007). It will be necessary to conduct more in-depth studies to assess the phylogenetic distribution of gBGC in plants and other hermaphrodite organisms and to further test the genomic Achilles'' heel hypothesis in relation to mating systems. While theoretically possible, the quantitative effect of gBGC on the evolution of mating systems remains a new, open, and challenging question.

Conclusion:

I showed that the interaction between gBGC and selection might have surprising qualitative consequences on load and inbreeding depression patterns. Given the few quantitative data available on gBGC levels and selection intensities (mainly in humans), it turns out that even weak genome-wide gBGC can have significant fitness impacts. gBGC should be taken into account not only for sequence analyses (Berglund et al. 2009; Galtier et al. 2009), but also for its potential fitness consequences, for instance concerning genetic diseases. Interferences between gBGC and selection also give rise to new questions on the evolution of mating systems. However, most of the challenging conclusions given here have yet to be quantitatively evaluated. Quantification of gBGC and its interaction with selection in various organisms will be crucial in the future.     

Coalescent Simulation of Intracodon Recombination

The coalescent with recombination is a very useful tool in molecular population genetics. Under this framework, genealogies often represent the evolution of the substitution unit, and because of this, the few coalescent algorithms implemented for the simulation of coding sequences force recombination to occur only between codons. However, it is clear that recombination is expected to occur most often within codons. Here we have developed an algorithm that can evolve coding sequences under an ancestral recombination graph that represents the genealogies at each nucleotide site, thereby allowing for intracodon recombination. The algorithm is a modification of Hudson''s coalescent in which, in addition to keeping track of events occurring in the ancestral material that reaches the sample, we need to keep track of events occurring in ancestral material that does not reach the sample but that is produced by intracodon recombination. We are able to show that at typical substitution rates the number of nonsynonymous changes induced by intracodon recombination is small and that intracodon recombination does not generally result in inflated estimates of the overall nonsynonymous/synonymous substitution ratio (ω). On the other hand, recombination can bias the estimation of ω at particular codons, resulting in apparent rate variation among sites and in the spurious identification of positively selected sites. Importantly, in this case, allowing for variable synonymous rates across sites greatly reduces the false-positive rate and recovers statistical power. Finally, coalescent simulations with intracodon recombination could be used to better represent the evolution of nuclear coding genes or fast-evolving pathogens such as HIV-1.We have implemented this algorithm in a computer program called NetRecodon, freely available at http://darwin.uvigo.es.THE coalescent (Kingman 1982; Hudson 1990) provides an efficient sampling of genealogical histories from a theoretical population evolving under a neutral Wright–Fisher model (Ewens 1979; Kingman 1982; Hudson 1990). Coalescent simulations are commonly used in molecular population genetics to understand the behavior and interactions among evolutionary processes under different scenarios (Innan et al. 2005), such as hypothesis testing (DeChaine and Martin 2006), evaluation and comparison of different analytical methods (Carvajal-Rodriguez et al. 2006), or estimation of population genetic parameters (Beaumont et al. 2002). Indeed, to obtain meaningful biological inferences from these simulations, it is very important that the underlying model is as realistic as possible. In this regard, a number of models have been developed during the last decade that consider different evolutionary processes such as recombination (Simonsen and Churchill 1997; Wiuf and Posada 2003), gene conversion (Wiuf and Hein 2000), selection (Hudson and Kaplan 1988, 1995), and gene flow or demographic history (Slatkin 1987; Pybus and Rambaut 2002).Despite these advances, and in the face of a plethora of coalescent simulators (Excoffier et al. 2000; Hudson 2002; Posada and Wiuf 2003; Spencer and Coop 2004; Mailund et al. 2005; Schaffner et al. 2005; Marjoram and Wall 2006; Arenas and Posada 2007; Hellenthal and Stephens 2007; Liang et al. 2007), it was not possible until very recently to simulate recombining protein-coding DNA sequences within this framework (Anisimova et al. 2003; Arenas and Posada 2007). Importantly, to our knowledge, the algorithms described or implemented so far allow recombination only between codons, not within them. The reason for this unrealistic constraint is that standard codon models describe the probabilities of change along a lineage from one codon to another (Yang 2006), whereas recombination can occur between any two nucleotides, potentially resulting in one or more lineages not being shared by all the positions of the codon. In other words, although the unit for substitution in coding sequences is the codon, the unit for recombination in these sequences is still the nucleotide. Here we describe a new algorithm that overcomes this limitation by allowing for the evolution of different positions of the same codon in distinct genealogies. Furthermore, we use this algorithm to evaluate the effect of intracodon recombination on the generation of nonsynonymous (NS) diversity and on the estimation of the ratio of nonsynonymous-to-synonymous substitution rates (ω or d_N/d_S) (Li and Gojobori 1983) and the hypotheses derived from it.     

Adaptive Divergence in Experimental Populations of Pseudomonas fluorescens. IV. Genetic Constraints Guide Evolutionary Trajectories in a Parallel Adaptive Radiation

The capacity for phenotypic evolution is dependent upon complex webs of functional interactions that connect genotype and phenotype. Wrinkly spreader (WS) genotypes arise repeatedly during the course of a model Pseudomonas adaptive radiation. Previous work showed that the evolution of WS variation was explained in part by spontaneous mutations in wspF, a component of the Wsp-signaling module, but also drew attention to the existence of unknown mutational causes. Here, we identify two new mutational pathways (Aws and Mws) that allow realization of the WS phenotype: in common with the Wsp module these pathways contain a di-guanylate cyclase-encoding gene subject to negative regulation. Together, mutations in the Wsp, Aws, and Mws regulatory modules account for the spectrum of WS phenotype-generating mutations found among a collection of 26 spontaneously arising WS genotypes obtained from independent adaptive radiations. Despite a large number of potential mutational pathways, the repeated discovery of mutations in a small number of loci (parallel evolution) prompted the construction of an ancestral genotype devoid of known (Wsp, Aws, and Mws) regulatory modules to see whether the types derived from this genotype could converge upon the WS phenotype via a novel route. Such types—with equivalent fitness effects—did emerge, although they took significantly longer to do so. Together our data provide an explanation for why WS evolution follows a limited number of mutational pathways and show how genetic architecture can bias the molecular variation presented to selection.UNDERSTANDING—and importantly, predicting—phenotypic evolution requires knowledge of the factors that affect the translation of mutation into phenotypic variation—the raw material of adaptive evolution. While much is known about mutation rate (e.g., Drake et al. 1998; Hudson et al. 2002), knowledge of the processes affecting the translation of DNA sequence variation into phenotypic variation is minimal.Advances in knowledge on at least two fronts suggest that progress in understanding the rules governing the generation of phenotypic variation is possible (Stern and Orgogozo 2009). The first stems from increased awareness of the genetic architecture underlying specific adaptive phenotypes and recognition of the fact that the capacity for evolutionary change is likely to be constrained by this architecture (Schlichting and Murren 2004; Hansen 2006). The second is the growing number of reports of parallel evolution (e.g., Pigeon et al. 1997; ffrench-Constant et al. 1998; Allender et al. 2003; Colosimo et al. 2004; Zhong et al. 2004; Boughman et al. 2005; Shindo et al. 2005; Kronforst et al. 2006; Woods et al. 2006; Zhang 2006; Bantinaki et al. 2007; McGregor et al. 2007; Ostrowski et al. 2008)—that is, the independent evolution of similar or identical features in two or more lineages—which suggests the possibility that evolution may follow a limited number of pathways (Schluter 1996). Indeed, giving substance to this idea are studies that show that mutations underlying parallel phenotypic evolution are nonrandomly distributed and typically clustered in homologous genes (Stern and Orgogozo 2008).While the nonrandom distribution of mutations during parallel genetic evolution may reflect constraints due to genetic architecture, some have argued that the primary cause is strong selection (e.g., Wichman et al. 1999; Woods et al. 2006). A means of disentangling the roles of population processes (selection) from genetic architecture is necessary for progress (Maynard Smith et al. 1985; Brakefield 2006); also necessary is insight into precisely how genetic architecture might bias the production of mutations presented to selection.Despite their relative simplicity, microbial populations offer opportunities to advance knowledge. The wrinkly spreader (WS) morphotype is one of many different niche specialist genotypes that emerge when experimental populations of Pseudomonas fluorescens are propagated in spatially structured microcosms (Rainey and Travisano 1998). Previous studies defined, via gene inactivation, the essential phenotypic and genetic traits that define a single WS genotype known as LSWS (Spiers et al. 2002, 2003) (Figure 1). LSWS differs from the ancestral SM genotype by a single nonsynonymous nucleotide change in wspF. Functionally (see Figure 2), WspF is a methyl esterase and negative regulator of the WspR di-guanylate cyclase (DGC) (Goymer et al. 2006) that is responsible for the biosynthesis of c-di-GMP (Malone et al. 2007), the allosteric activator of cellulose synthesis enzymes (Ross et al. 1987). The net effect of the wspF mutation is to promote physiological changes that lead to the formation of a microbial mat at the air–liquid interface of static broth microcosms (Rainey and Rainey 2003).Open in a separate windowFigure 1.—Outline of experimental strategy for elucidation of WS-generating mutations and their subsequent identity and distribution among a collection of independently evolved, spontaneously arising WS genotypes. The strategy involves, first, the genetic analysis of a specific WS genotype (e.g., LSWS) to identify the causal mutation, and second, a survey of DNA sequence variation at specific loci known to harbor causal mutations among a collection of spontaneously arising WS genotypes. For example, suppressor analysis of LSWS using a transposon to inactivate genes necessary for expression of the wrinkly morphology delivered a large number of candidate genes (top left) (Spiers et al. 2002). Genetic and functional analysis of these candidate genes (e.g., Goymer et al. 2006) led eventually to the identity of the spontaneous mutation (in wspF) responsible for the evolution of LSWS from the ancestral SM genotype (Bantinaki et al. 2007). Subsequent analysis of the wspF sequence among 26 independent WS genotypes (bottom) showed that 50% harbored spontaneous mutations (of different kinds; see Open in a separate windowFigure 2.—Network diagram of DGC-encoding pathways underpinning the evolution of the WS phenotype and their regulation. Overproduction of c-di-GMP results in overproduction of cellulose and other adhesive factors that determine the WS phenotype. The ancestral SBW25 genome contains 39 putative DGCs, each in principle capable of synthesizing the production of c-di-GMP, and yet WS genotypes arise most commonly as a consequence of mutations in just three DGC-containing pathways: Wsp, Aws, and Mws. In each instance, the causal mutations are most commonly in the negative regulatory component: wspF, awsX, and the phosphodiesterase domain of mwsR (see text).To determine whether spontaneous mutations in wspF are a common cause of the WS phenotype, the nucleotide sequence of this gene was obtained from a collection of 26 spontaneously arising WS genotypes (WS_A-Z) taken from 26 independent adaptive radiations, each founded by the same ancestral SM genotype (Figure 1): 13 contained mutations in wspF (Bantinaki et al. 2007). The existence of additional mutational pathways to WS provided the initial motivation for this study.

TABLE 1

Mutational causes of WS

The Conserved miR-51 microRNA Family Is Redundantly Required for Embryonic Development and Pharynx Attachment in Caenorhabditis elegans

A major question about cytokinesis concerns the role of the septin proteins, which localize to the division site in all animal and fungal cells but are essential for cytokinesis only in some cell types. For example, in Schizosaccharomyces pombe, four septins localize to the division site, but deletion of the four genes produces only a modest delay in cell separation. To ask if the S. pombe septins function redundantly in cytokinesis, we conducted a synthetic-lethal screen in a septin-deficient strain and identified seven mutations. One mutation affects Cdc4, a myosin light chain that is an essential component of the cytokinetic actomyosin ring. Five others cause frequent cell lysis during cell separation and map to two loci. These mutations and their dosage suppressors define a signaling pathway (including Rho1 and a novel arrestin) for repairing cell-wall damage. The seventh mutation affects the poorly understood RNA-binding protein Scw1 and severely delays cell separation when combined either with a septin mutation or with a mutation affecting the septin-interacting, anillin-like protein Mid2, suggesting that Scw1 functions in a pathway parallel to that of the septins. Taken together, our results suggest that the S. pombe septins participate redundantly in one or more pathways that cooperate with the actomyosin ring during cytokinesis and that a septin defect causes septum defects that can be repaired effectively only when the cell-integrity pathway is intact.THE fission yeast Schizosaccharomyces pombe provides an outstanding model system for studies of cytokinesis (McCollum and Gould 2001; Balasubramanian et al. 2004; Pollard and Wu 2010). As in most animal cells, successful cytokinesis in S. pombe requires an actomyosin ring (AMR). The AMR begins to assemble at the G₂/M transition and involves the type II myosin heavy chains Myo2 and Myp2 and the light chains Cdc4 and Rlc1 (Wu et al. 2003). Myo2 and Cdc4 are essential for cytokinesis under all known conditions, Rlc1 is important at all temperatures but essential only at low temperatures, and Myp2 is essential only under stress conditions. As the AMR constricts, a septum of cell wall is formed between the daughter cells. The primary septum is sandwiched by secondary septa and subsequently digested to allow cell separation (Humbel et al. 2001; Sipiczki 2007). Because of the internal turgor pressure of the cells, the proper assembly and structural integrity of the septal layers are essential for cell survival.Septum formation involves the β-glucan synthases Bgs1/Cps1/Drc1, Bgs3, and Bgs4 (Ishiguro et al. 1997; Le Goff et al. 1999; Liu et al. 1999, 2002; Martín et al. 2003; Cortés et al. 2005) and the α-glucan synthase Ags1/Mok1 (Hochstenbach et al. 1998; Katayama et al. 1999). These synthases are regulated by the Rho GTPases Rho1 and Rho2 and the protein kinase C isoforms Pck1 and Pck2 (Arellano et al. 1996, 1997, 1999; Nakano et al. 1997; Hirata et al. 1998; Calonge et al. 2000; Sayers et al. 2000; Ma et al. 2006; Barba et al. 2008; García et al. 2009b). The Rho GTPases themselves appear to be regulated by both GTPase-activating proteins (GAPs) and guanine-nucleotide-exchange factors (GEFs) (Nakano et al. 2001; Calonge et al. 2003; Iwaki et al. 2003; Tajadura et al. 2004; Morrell-Falvey et al. 2005; Mutoh et al. 2005; García et al. 2006, 2009a,b). In addition, septum formation and AMR function appear to be interdependent. In the absence of a normal AMR, cells form aberrant septa and/or deposit septal materials at random locations, whereas a mutant defective in septum formation (bgs1) is also defective in AMR constriction (Gould and Simanis 1997; Le Goff et al. 1999; Liu et al. 1999, 2000). Both AMR constriction and septum formation also depend on the septation initiation network involving the small GTPase Spg1 (McCollum and Gould 2001; Krapp and Simanis 2008). Despite this considerable progress, many questions remain about the mechanisms and regulation of septum formation and its relationships to the function of the AMR.One major question concerns the role(s) of the septins. Proteins of this family are ubiquitous in fungal and animal cells and typically localize to the cell cortex, where they appear to serve as scaffolds and diffusion barriers for other proteins that participate in a wide variety of cellular processes (Longtine et al. 1996; Gladfelter et al. 2001; Hall et al. 2008; Caudron and Barral 2009). Despite the recent progress in elucidating the mechanisms of septin assembly (John et al. 2007; Sirajuddin et al. 2007; Bertin et al. 2008; McMurray and Thorner 2008), the details of septin function remain obscure. However, one prominent role of the septins and associated proteins is in cytokinesis. Septins concentrate at the division site in every cell type that has been examined, and in Saccharomyces cerevisiae (Hartwell 1971; Longtine et al. 1996; Lippincott et al. 2001; Dobbelaere and Barral 2004) and at least some Drosophila (Neufeld and Rubin 1994; Adam et al. 2000) and mammalian (Kinoshita et al. 1997; Surka et al. 2002) cell types, the septins are essential for cytokinesis. In S. cerevisiae, the septins are required for formation of the AMR (Bi et al. 1998; Lippincott and Li 1998). However, this cannot be their only role, because the AMR itself is not essential for cytokinesis in this organism (Bi et al. 1998; Korinek et al. 2000; Schmidt et al. 2002). Moreover, there is no evidence that the septins are necessary for AMR formation or function in any other organism. A further complication is that in some cell types, including most Caenorhabditis elegans cells (Nguyen et al. 2000; Maddox et al. 2007) and some Drosophila cells (Adam et al. 2000; Field et al. 2008), the septins do not appear to be essential for cytokinesis even though they localize to the division site.S. pombe has seven septins, four of which (Spn1, Spn2, Spn3, and Spn4) are expressed in vegetative cells and localize to the division site shortly before AMR constriction and septum formation (Longtine et al. 1996; Berlin et al. 2003; Tasto et al. 2003; Wu et al. 2003; An et al. 2004; Petit et al. 2005; Pan et al. 2007; Onishi et al. 2010). Spn1 and Spn4 appear to be the core members of the septin complex (An et al. 2004; McMurray and Thorner 2008), and mutants lacking either of these proteins do not assemble the others at the division site. Assembly of a normal septin ring also depends on the anillin-like protein Mid2, which colocalizes with the septins (Berlin et al. 2003; Tasto et al. 2003). Surprisingly, mutants lacking the septins are viable and form seemingly complete septa with approximately normal timing. These mutants do, however, display a variable delay in separation of the daughter cells, suggesting that the septins play some role(s) in the proper completion of the septum or in subsequent processes necessary for cell separation (Longtine et al. 1996; An et al. 2004; Martín-Cuadrado et al. 2005).It is possible that the septins localize to the division site and yet are nonessential for division in some cell types because their role is redundant with that of some other protein(s) or pathway(s). To explore this possibility in S. pombe, we screened for mutations that were lethal in combination with a lack of septins. The results suggest that the septins cooperate with the AMR during cytokinesis and that, in the absence of septin function, the septum is not formed properly, so that an intact system for recognizing and repairing cell-wall damage becomes critical for cell survival.     

Immunomodulation by Mesenchymal Stem Cells in Veterinary Species

Mesenchymal stem cells (MSC) are adult-derived multipotent stem cells that have been derived from almost every tissue. They are classically defined as spindle-shaped, plastic-adherent cells capable of adipogenic, chondrogenic, and osteogenic differentiation. This capacity for trilineage differentiation has been the foundation for research into the use of MSC to regenerate damaged tissues. Recent studies have shown that MSC interact with cells of the immune system and modulate their function. Although many of the details underlying the mechanisms by which MSC modulate the immune system have been defined for human and rodent (mouse and rat) MSC, much less is known about MSC from other veterinary species. This knowledge gap is particularly important because the clinical use of MSC in veterinary medicine is increasing and far exceeds the use of MSC in human medicine. It is crucial to determine how MSC modulate the immune system for each animal species as well as for MSC derived from any given tissue source. A comparative approach provides a unique translational opportunity to bring novel cell-based therapies to the veterinary market as well as enhance the utility of animal models for human disorders. The current review covers what is currently known about MSC and their immunomodulatory functions in veterinary species, excluding laboratory rodents.Abbreviations: AT, adipose tissue; BM, Bone marrow; CB, umbilical cord blood; CT, umbilical cord tissue; DC, dendritic cell; IDO, indoleamine 2;3-dioxygenase; MSC, mesenchymal stem cells; PGE₂, prostaglandin E2; VEGF, vascular endothelial growth factorMesenchymal stem cells (MSC, alternatively known as mesenchymal stromal cells) were first reported in the literature in 1968.³⁹ MSC are thought to be of pericyte origin (cells that line the vasculature)^21,22 and typically are isolated from highly vascular tissues. In humans and mice, MSC have been isolated from fat, placental tissues (placenta, Wharton jelly, umbilical cord, umbilical cord blood), hair follicles, tendon, synovial membrane, periodontal ligament, and every major organ (brain, spleen, liver, kidney, lung, bone marrow, muscle, thymus, pancreas, skin).^23,121 For most current clinical applications, MSC are isolated from adipose tissue (AT), bone marrow (BM), umbilical cord blood (CB), and umbilical cord tissue (CT; 11,87,99 Clinical trials in human medicine focus on the use of MSC both for their antiinflammatory properties (graft-versus-host disease, irritable bowel syndrome) and their ability to aid in tissue and bone regeneration in combination with growth factors and bone scaffolds (clinicaltrials.gov).¹³¹ For tissue regeneration, the abilities of MSC to differentiate and to secrete mediators and interact with cells of the immune system likely contribute to tissue healing (Figure 1). The current review will not address the specific use of MSC for orthopedic applications and tissue regeneration, although the topic is covered widely in current literature for both human and veterinary medicine.^57,62,90

Table 1.

Tissues from which MSC have been isolated

Topologies of the Conditional Ancestral Trees and Full-Likelihood-Based Inference in the General Coalescent Tree Framework

The general coalescent tree framework is a family of models for determining ancestries among random samples of DNA sequences at a nonrecombining locus. The ancestral models included in this framework can be derived under various evolutionary scenarios. Here, a computationally tractable full-likelihood-based inference method for neutral polymorphisms is presented, using the general coalescent tree framework and the infinite-sites model for mutations in DNA sequences. First, an exact sampling scheme is developed to determine the topologies of conditional ancestral trees. However, this scheme has some computational limitations and to overcome these limitations a second scheme based on importance sampling is provided. Next, these schemes are combined with Monte Carlo integrations to estimate the likelihood of full polymorphism data, the ages of mutations in the sample, and the time of the most recent common ancestor. In addition, this article shows how to apply this method for estimating the likelihood of neutral polymorphism data in a sample of DNA sequences completely linked to a mutant allele of interest. This method is illustrated using the data in a sample of DNA sequences at the APOE gene locus.THE interest in analyzing polymorphism data in contemporary samples of DNA sequences under various evolutionary scenarios creates a demand to design computationally tractable full-likelihood-based inference methods. For an evolutionary scenario of interest, an ancestral-mutation model can be used to design such a method. The ancestral-mutation model for a sample of DNA sequences at a nonrecombining locus is a combination of two processes: one is an ancestral process that traces the lineages of the sample back in time until the most recent common ancestor, constructing an ancestral tree for the sample. The second is a mutation process that is superimposed on the ancestral tree. The complexities of ancestral-mutation models make the design of such methods challenging. Full data are used instead of summary statistics, which can result in loss of important information in the data (see Felsenstein 1992; Donnelly and Tavaré 1995). In addition, current methods use specific features of the underlying ancestral-mutation models, so they lose flexibility to be applicable to other ancestral-mutation models.More specifically, Griffiths and Tavaré (1994c, 1995) and Kuhner et al. (1995) developed full-likelihood-based inference methods for neutral polymorphisms at a nonrecombining locus. They used the combinations of the standard coalescent (Kingman 1982a,b,c; Hudson 1983; Tajima 1983) with the finite-sites or infinite-sites (Watterson 1975) models as ancestral-mutation models. Stephens and Donnelly (2000) designed an importance sampling method to estimate the full likelihood of the data using the same settings for the ancestral-mutation models. Hobolth et al. (2008) provided another importance sampling scheme restricted to the infinite-sites model. The last two methods are computationally more efficient than the first two methods, but they lose flexibility to be applicable to ancestral models without standard coalescent features with independent coalescence waiting times, such as the coalescent processes with exponential growth (Slatkin and Hudson 1991; Griffiths and Tavaré 1994b).To incorporate the coalescent processes with exponential growth, Kuhner et al. (1998) and Griffiths and Tavaré (1994a, 1999) extended their previous methods. For example, the method of Griffiths and Tavaré (1994a, 1999) allows one to consider ancestral models based on coalescent processes with variable population sizes. Coop and Griffiths (2004) modified this inference method and made it applicable for analyzing full polymorphism data in a sample of DNA sequences from a nonrecombining locus completely linked to a mutant allele of interest, either neutral or under selection. Additionally, ancestral models have been developed for this type of sample, where the mutant allele is either neutral (Griffiths and Tavaré 1998, 2003; Wiuf and Donnelly 1999; Stephens 2000) or under selection (Slatkin and Rannala 1997; Stephens and Donnelly 2003). The ancestral model of Slatkin and Rannala (1997) is part of a family of ancestral models derived by Thompson (1975), Nee et al. (1994), and Rannala (1997), using a linear birth–death process as an evolutionary process in a population. Although all the ancestral models mentioned above differ in their properties and evolutionary scenarios, they are part of the general coalescent tree framework (Griffiths and Tavaré 1998). Therefore, a computationally tractable full-likelihood-based inference method based on this general framework is of great interest.For a sample of n sequences, an ancestral model in the general coalescent tree framework is described as a bifurcating rooted tree with n − 1 internal nodes and n leaves, where the internal nodes are coalescent events that happen one at a time. The tree is a combination of two independent components: the topology and the branch lengths. The topology of the tree is constructed going backward in time by combining two randomly chosen ancestral lineages of the sample at each node; the branch lengths of the tree are defined by the joint distribution function of the coalescence waiting times. Note that any density function for coalescence waiting times can define an ancestral model in the general coalescent tree framework.The n leaves (and the sequences in the sample) are labeled from 1 to n; and the n − 1 internal nodes of the ancestral tree are labeled from 1 to n − 1 (in order of occurrence of the coalescent events backward in time). Thus, the topology of an ancestral tree is a leaf-labeled bifurcating rooted tree with totally ordered interior vertices. These trees are called topological trees.When using the general coalescent tree framework and the infinite-sites model, an evolutionary process that generates polymorphism data in a sample of DNA sequences can be described in the following way. An ancestral tree is constructed, as described above, and mutations are added independently on different branches of the ancestral tree as Poisson processes with equal rates, θ/2, in which θ is the mutation rate at the locus. Then, at the mutation events, the ancestral sequences of the sample are changed according to the infinite-sites model; that is, each mutation occurs at a site of an ancestral sequence at which no previous mutations occurred. Thus, these changes define polymorphism data.Naively, this probabilistic framework can be used to estimate the likelihood of the full observed data in a sample of n sequences. That is, data sets are simulated independently as described above and each simulated data set is compared to the observed data. The proportion of the simulated data sets that match the observed data is an estimate of the likelihood of the observed data. Although this approach provides an estimate for the likelihood of the observed data, this method is computationally infeasible, because the topologies of the ancestral trees of the generated data sets are sampled from the space of all the possible topological trees with n leaves. This space has size n!(n − 1)!/2ⁿ⁻¹ (Edwards 1970), which is huge for moderate values of n. The topologies of the ancestral trees of the generated data sets that match the observed data represent a small portion of that space. Thus, designing a method that samples topologies of the ancestral trees from this subspace can make the method computationally tractable.On the basis of this idea, I use the general coalescent tree framework with the infinite-sites model to develop a computationally tractable full-likelihood-based inference method for polymorphisms in DNA sequences at a nonrecombining locus. First, an exact sampling scheme for topologies of the conditional ancestral trees is developed. This method has some computational limitations, so to overcome these limitations a second scheme based on an importance sampling is provided. These sampling schemes are combined with Monte Carlo integrations to estimate the likelihood of the full data, the ages of the mutations in the sample, and the time of the most recent common ancestor of the sample. I describe an application of this method for neutral polymorphism data in a sample of DNA sequences at a nonrecombining locus that is completely linked to a mutant allele of interest, either neutral or under selection. The method is illustrated using the data in a sample of DNA sequences at the APOE gene locus from Fullerton et al. (2000).     

Polymorphic Genes of Major Effect: Consequences for Variation,Selection and Evolution in Arabidopsis thaliana

The importance of genes of major effect for evolutionary trajectories within and among natural populations has long been the subject of intense debate. For example, if allelic variation at a major-effect locus fundamentally alters the structure of quantitative trait variation, then fixation of a single locus can have rapid and profound effects on the rate or direction of subsequent evolutionary change. Using an Arabidopsis thaliana RIL mapping population, we compare G-matrix structure between lines possessing different alleles at ERECTA, a locus known to affect ecologically relevant variation in plant architecture. We find that the allele present at ERECTA significantly alters G-matrix structure—in particular the genetic correlations between branch number and flowering time traits—and may also modulate the strength of natural selection on these traits. Despite these differences, however, when we extend our analysis to determine how evolution might differ depending on the ERECTA allele, we find that predicted responses to selection are similar. To compare responses to selection between allele classes, we developed a resampling strategy that incorporates uncertainty in estimates of selection that can also be used for statistical comparisons of G matrices.THE structure of the genetic variation that underlies phenotypic traits has important consequences for understanding the evolution of quantitative traits (Fisher 1930; Lande 1979; Bulmer 1980; Kimura 1983; Orr 1998; Agrawal et al. 2001). Despite the infinitesimal model''s allure and theoretical tractability (see Orr and Coyne 1992; Orr 1998, 2005a,b for reviews of its influence), evidence has accumulated from several sources (artificial selection experiments, experimental evolution, and QTL mapping) to suggest that genes of major effect often contribute to quantitative traits. Thus, the frequency and role of genes of major effect in evolutionary quantitative genetics have been a subject of intense debate and investigation for close to 80 years (Fisher 1930; Kimura 1983; Orr 1998, 2005a,b). Beyond the conceptual implications, the prevalence of major-effect loci also affects our ability to determine the genetic basis of adaptations and species differences (e.g., Bradshaw et al. 1995, 1998).Although the existence of genes of major effect is no longer in doubt, we still lack basic empirical data on how segregating variation at such genes affects key components of evolutionary process (but see Carrière and Roff 1995). In other words, How does polymorphism at genes of major effect alter patterns of genetic variation and covariation, natural selection, and the likely response to selection? The lack of data stems, in part, from the methods used to detect genes of major effect: experimental evolution (e.g., Bull et al. 1997; Zeyl 2005) and QTL analysis (see Erickson et al. 2004 for a review) often detect such genes retrospectively after they have become fixed in experimental populations or the species pairs used to generate the mapping population. The consequences of polymorphism at these genes on patterns of variation, covariation, selection, and the response to selection—which can be transient (Agrawal et al. 2001)—are thus often unobserved.A partial exception to the absence of data on the effects of major genes comes from artificial selection experiments, in which a substantial evolutionary response to selection in the phenotype after a plateau is often interpreted as evidence for the fixation of a major-effect locus (Frankham et al. 1968; Yoo 1980a,b; Frankham 1980; Shrimpton and Robertson 1988a,b; Caballero et al. 1991; Keightley 1998; see Mackay 1990 and Hill and Caballero 1992 for reviews). However, many of these experiments report only data on the selected phenotype (e.g., bristle number) or, alternatively, the selected phenotype and some measure of fitness (e.g., Frankham et al. 1968, Yoo 1980b; Caballero et al. 1991; Mackay et al. 1994; Fry et al. 1995; Nuzhdin et al. 1995; Zur Lage et al. 1997), making it difficult to infer how a mutation will affect variation, covariation, selection, and evolutionary responses for a suite of traits that might affect fitness themselves. One approach is to document how variation at individual genes of major effect affects the genetic variance–covariance matrix (“G matrix”; Lande 1979), which represents the additive genetic variance and covariance between traits.Although direct evidence for variation at major-effect genes altering patterns of genetic variation, covariation, and selection is rare, there is abundant evidence for the genetic mechanisms that could produce these dynamics. A gene of major effect could have these consequences due to any of at least three genetic mechanisms: (1) pleiotropy, where a gene of major effect influences several traits, including potentially fitness, simultaneously, (2) physical linkage or linkage disequilibrium (LD), in which a gene of major effect is either physically linked or in LD with other genes that influence other traits under selection, and (3) epistasis, in which the allele present at a major-effect gene alters the phenotypic effect of other loci and potentially phenotypes under selection. Evidence for these three evolutionary genetic mechanisms leading to changes in suites of traits comes from a variety of sources, including mutation accumulation experiments (Clark et al. 1995; Fernandez and Lopez-Fanjul 1996), mutation induction experiments (Keightley and Ohnishi 1998), artificial selection experiments (Long et al. 1995), and transposable element insertions (Rollmann et al. 2006). For pleiotropy in particular, major-effect genes that have consequences on several phenotypic traits are well known from the domestication and livestock breeding literature [e.g., myostatin mutations in Belgian blue cattle and whippets (Arthur 1995; Grobet et al. 1997; Mosher et al. 2007), halothane genes in pigs (Christian and Rothschild 1991; Fujii et al. 1991), and Booroola and Inverdale genes in sheep (Amer et al. 1999; Visscher et al. 2000)]. While these data suggest that variation at major-effect genes could—and probably does—influence variation, covariation, and selection on quantitative traits, data on the magnitude of these consequences remain lacking.Recombinant inbred line (RIL) populations are a promising tool for investigating the influence of major-effect loci. During advancement of the lines from F₂''s to RILs, alternate alleles at major-effect genes (and most of the rest of the genome) will be made homozygous, simplifying comparisons among genotypic classes. Because of the high homozygosity, individuals within RILs are nearly genetically identical, facilitating phenotyping of many genotypes under a range of environments. In addition, because of recombination, alternative alleles are randomized across genetic backgrounds—facilitating robust comparisons between sets of lines differing at a major-effect locus.Here we investigate how polymorphism at an artificially induced mutation, the erecta locus in Arabidopsis thaliana, affects the magnitude of these important evolutionary genetic parameters under ecologically realistic field conditions. We use the Landsberg erecta (Ler) × Columbia (Col) RIL population of A. thaliana to examine how variation at a gene of major effect influences genetic variation, covariation, and selection on quantitative traits in a field setting. The Ler × Col RIL population is particularly suitable, because it segregates for an artificially induced mutation at the erecta locus, which has been shown to influence a wide variety of plant traits. The Ler × Col population thus allows a powerful test of the effects of segregating variation at a gene—chosen a priori—with numerous pleiotropic effects. The ERECTA gene is a leucine-rich receptor-like kinase (LRR-RLK) (Torii et al. 1996) and has been shown to affect plant growth rates (El-Lithy et al. 2004), stomatal patterning and transpiration efficiency (Masle et al. 2005; Shpak et al. 2005), bacterial pathogen resistance (Godiard et al. 2003), inflorescence and floral organ size and shape (Douglas et al. 2002; Shpak et al. 2003, 2004), and leaf polarity (Xu et al. 2003; Qi et al. 2004).Specifically, we sought to answer the following questions: (1) Is variation at erecta significantly associated with changes to the G matrix? (2) Is variation at erecta associated with changes in natural selection on genetically variable traits? And (3) is variation at erecta associated with significantly different projected evolutionary responses to selection?     

Directionality of Epistasis in a Murine Intercross Population

Directional epistasis describes a situation in which epistasis consistently increases or decreases the effect of allele substitutions, thereby affecting the amount of additive genetic variance available for selection in a given direction. This study applies a recent parameterization of directionality of epistasis to empirical data. Data stems from a QTL mapping study on an intercross between inbred mouse (Mus musculus) strains LG/J and SM/J, originally selected for large and small body mass, respectively. Results show a negative average directionality of epistasis for body-composition traits, predicting a reduction in additive allelic effects and in the response to selection for increased size. Focusing on average modification of additive effect of single loci, we find a more complex picture, whereby the effects of some loci are enhanced consistently across backgrounds, while effects of other loci are decreased, potentially contributing to either enhancement or reduction of allelic effects when selection acts at single loci. We demonstrate and discuss how the interpretation of the overall measurement of directionality depends on the complexity of the genotype–phenotype map. The measure of directionality changes with the power of scale in a predictable way; however, its expected effect with respect to the modification of additive genetic effects remains constant.EPISTASIS is present when the effect of a genetic substitution depends on the genotypes at other loci. At the population level, this means that average allelic effects change as allele frequencies at other loci change, and thus that gene effects can evolve. The evolutionary significance of epistasis has been recognized mainly in relation to the allele-frequency changes that are caused by genetic drift (e.g., Goodnight 1987, 1988, 1995; Cheverud and Routman 1996; Barton and Turelli 2004; De Brito et al. 2005; Turelli and Barton 2006), whereas the epistatic effects under directional selection have been treated only recently (e.g., Carter et al. 2005; Weinreich et al. 2005; Carlborg et al. 2006; Hansen et al. 2006; Yukilevich et al. 2008). Technically, the effects of epistatic interaction can be considered as having two aspects: the architecture itself (i.e., the existence of a nonadditive component, the so-called functional aspect, see below), and the effect of allele frequency on genetic variance (i.e., the statistical aspect). An additional consideration is crucial for the response to selection of a given trait, namely that the response is generated by the joint action of many epistatic interactions. Each of the interactions can either enhance or diminish the additive genetic effect in any specific phenotypic dimension. Their composite effect depends on the pattern, i.e., whether the effects accumulate or cancel each other out (Hansen and Wagner 2001a,b; Carter et al. 2005; Hansen et al. 2006). In the following we provide a brief general account of epistasis and then focus on the effect of its composite pattern, the empirical assessment of which is the goal of this study.The traditional population-genetic approach to selection response initially emphasized additive genetic variance and treated any variance unexplained by the additive effects, including variance due to interactions within or between loci, as residual variance (Fisher 1918). Later this model was extended to account for epistasis (Cockerham 1954, Kempthorne 1954). The interaction component of this residual variance is dependent on population allele frequencies at the interacting loci. Starting with nonadditive effects within a single locus (dominance), Falconer (1960) described the effect of allele frequencies on the statistical measure of average allelic effect. Cheverud and Routman (1995) explored the analogous effect at the two-locus level, leading to distinction between allele-frequency-dependent statistical epistasis (contributing to epistatic variance) on population level, and allele-frequency-independent physiological (or functional) epistasis on the individual level, which contributes to all the genetic variance components, i.e.,, additive, dominance, and epistatic. Thus while physiological epistasis describes the genetic architecture of a given phenotype defining the potential for epistatic effects, statistical epistasis describes the realization of individual-level epistatic effects in terms of allele-frequency-dependent genetic variance components. Several authors have worked out tools to estimate gene interaction effects at the individual level independently of the allele frequencies to distinguish between the physiological and statistical epistasis (Cheverud and Routman 1995; Wagner et al. 1998; Hansen and Wagner 2001a; Barton and Turelli 2004; Yang 2004; Zeng et al. 2005; Wang and Zeng 2006). These methods have been recently generalized in a common framework for measuring epistasis and translating between the population and individual levels (NOIA: Alvarez-Castro and Carlborg 2007; Alvarez-Castro et al. 2008; see Le Rouzic 2008 for R software package). An understanding of the pattern of epistatic architecture enables prediction of its effects at any given set of allele frequencies.The original physiological epistasis has often referred to isolated pairwise interactions. The genetic basis of a complex trait is affected by more than two interacting loci; thus the trait''s genetic variance is affected by the combination of these interaction effects. Effects at the individual loci can add up, or cancel out. For example, when the alleles at background loci become fixed (e.g., due to a bottleneck), previously background-dependent genetic effects at loci A and B can become additive effects of the same or of the opposite sign. Depending on the sign and size of their effects, the two allele substitutions at loci A and B therefore add up to increase, decrease, or have no overall effect on additive genetic variance. Hansen and Wagner (2001a) introduced the notion of directionality of epistasis to emphasize the importance of the pattern of epistatic architecture for the system''s evolvability. Directionality measures the consistency of epistatic effects on additive variance for a specific locus and trait across the genome, given a defined reference genotype. It describes whether epistasis tends to enhance or diminish the additive effects of interacting loci on a trait in a specified phenotypic direction. The population-dynamic studies analyzing effects of the epistatic pattern show that the directionality of epistasis can be a major determinant of evolution on time scales beyond a few generations (Hermisson et al. 2003; Carter et al. 2005; Hansen et al. 2006; Yukilevich et al. 2008; Alvarez-Castro et al. 2009; Fierst and Hansen 2010). Averaged across interactions, positive directional epistasis increases additive variance and the response to positive selection relative to that predicted by the additive genetic effects alone, while negative directional epistasis tends to decrease the response to selection in that same phenotypic direction. An absence of epistatic directionality occurs when positive and negative directional epistatic effects cancel out on average or, trivially, from the absence of epistasis. Thus directionality describes the local, population-specific curvature of the genotype–phenotype map (Figure 1).Open in a separate windowFigure 1.—Directionality of epistasis describes the local curvature of the genotype–phenotype map. The genotypic values of the trait of interest in the two parental inbred populations and in an intercross are plotted on y-axis. Given a defined direction on a trait axis and a reference point for measurement of genetic effects, the directionality () describes whether the epistasis increases or decreases the effect of allelic substitutions relative to the value predicted by additive effects alone. The extrapolation of the curvature beyond the local effects requires the knowledge of higher-order epistatic effects; however, the local effects can be calculated for different reference points and different phenotypic directions (see Hansen and Wagner 2001a,b).Aspects of epistatic directionality have been addressed previously with respect to the effect on fitness (reviewed in Phillips et al. 2000), with focus on the synergistic epistasis for deleterious effect enhancement (e.g., Kondrashov 1988; Charlesworth 1990; Hansen and Wagner 2001b). However, the empirical study of overall directionality of epistasis, and its measurement in morphological or physiological traits, is still lacking (Hansen 2006). Implicit indications from empirical studies are ambiguous (see examples in Carter et al. 2005). Here, we present an assessment of directionality of two-way epistasis between QTL in an intercross population of laboratory mice, paying special attention to scale effects.