DYNAMICS OF HYBRID INCOMPATIBILITY IN GENE NETWORKS IN A CONSTANT ENVIRONMENT (original) (raw)

Abstract

After an ancestral population splits into two allopatric populations, different mutations may fix in each. When pairs of mutations are brought together in a hybrid offspring, epistasis may cause reduced fitness. Such pairs are known as Bateson–Dobzhansky–Muller (BDM) incompatibilities. A well-known model of BDM incompatibility due to Orr suggests that the fitness load on hybrids should initially accelerate, and continue to increase as the number of potentially incompatible substitutions increases (the “snowball effect”). In the gene networks model, which violates a key assumption of Orr's model (independence of fixation probabilities), the snowball effect often does not occur. Instead, we describe three possible dynamics in a constant environment: (1) Stabilizing selection can constrain two allopatric populations to remain near-perfectly compatible. (2) Despite constancy of environment, punctuated evolution may obtain; populations may experience rare adaptations asynchronously, permitting incompatibility. (3) Despite stabilizing selection, developmental system drift may permit genetic change, allowing two populations to drift in and out of compatibility. We reinterpret Orr's model in terms of genetic distance. We extend Orr's model to the finite loci case, which can limit incompatibility. Finally, we suggest that neutral evolution of gene regulation in nature, to the point of speciation, is a distinct possibility.

Darwin (1859) compared attempts to cross two varieties of a single species (in many cases producing hybrid offspring of reduced fertility or viability) and attempts to cross two truly distinct species (producing almost no viable, or fertile, hybrid offspring). He came to the general conclusion that, “… there is no fundamental distinction between species and varieties,” in the sense that the two cases differ only in the degree to which similar mechanisms cause infertility or inviability (with rare exceptions). Darwin also recognized, but never solved, the problem of how hybrid infertility or inviability might increase in apparent opposition to selection, “inasmuch as the sterility of hybrids could not possibly be of any advantage to them, and therefore could not have been acquired by the continued preservation of successive profitable degrees of sterility.”Dobzhansky (1936) and Muller (1942) used a genetic model to solve the paradox. In this model, pairs of incompatible alleles can arise by two mutations in two allopatric populations. The important feature is that the alleles are deleterious only in combination, and that this combination is never “tested” by natural selection until the alleles are finally brought together in an unfortunate hybrid individual. Preexisting reproductive isolation between two populations would allow the accumulation of such incompatibilities, whose collective effect could make reproductive isolation permanent even if the original cause of isolation were removed. This model had been previously proposed by Bateson (1909), though apparently neither Dobzhansky nor Muller were aware of this (Orr 1996). Such pairs (or triplets, etc.) of alleles are therefore referred to as Bateson–Dobzhansky–Muller incompatibilities (BDMIs).

Orr (1995) derived a relationship between K, the number of substitutions since the ancestral population split into two allopatric populations, and L, the expected “fitness load” suffered by hybrids in terms of the average deleterious effect of each such potentially incompatible pair of alleles

1

For —the start of the speciation process—Orr suggests that the strength of reproductive isolation should increase nearly as fast as the square of the number of substitutions

2

This initial rapid acceleration of the fitness load on hybrids is known as the “snowball effect” (Orr and Turelli 2001). As speciation continues, (K increases) L asymptotically approaches one, corresponding to complete reproductive isolation. Orr's result does not refer to any particular mode of fixation (e.g., under stabilizing selection, directional selection, or neutral drift) but assumes that fixations of alleles are independent and equiprobable events, and that one locus can experience only one substitution in either population, which eliminates the possibility of parallel genetic evolution.

Here we explore some violations of Orr's assumptions: specifically, we allow the fixations of different alleles to have different probabilities, and those fixation probabilities not to be independent, but instead to depend on what fixations have previously occurred. In addition, we permit multiple substitutions across two populations at the same locus, thereby allowing parallel evolution. We do not seek to falsify Orr's (1995) results, but rather to extend them.

Our experimental model for this exploration is in silico populations of “gene networks.” This model has previously been used to study evolvability and canalization (Wagner 1996; Siegal and Bergman 2002). Gene networks model the expression dynamics of sets of genes that are up- and down-regulated by one another's gene products. The model permits complex epistasis and pleiotropy; hence fixation probabilities are not independent.

In this article, we consider the case of a constant environment, identical in the two allopatric populations. We choose this because the case of directional selection has already been covered well using a similar model: Johnson and Porter (Johnson and Porter 2000, 2007; Porter and Johnson 2002) showed that linear gene regulatory pathways can accumulate BDMIs in the presence of directional selection (a slowly ever-receding phenotypic optimum, which induces the population mean phenotype to “chase” the optimum). That directional selection would accelerate speciation is perhaps intuitive, because it results in substitutions, but theoretical predictions indicate that hybrid incompatibility should also increase due to substitutions occurring by drift alone, even under conditions of stabilizing selection (i.e., Orr's model does not specify the cause of the substitutions). Notably, Johnson and Porter did not observe accumulation of hybrid incompatibility under stabilizing selection. “The interaction of genetic drift and mutation, even at very high rates, did not reduce hybrid fitness at all on the time-scales we considered” (Johnson and Porter 2000). In contrast, using the gene networks model, we demonstrate the growth of hybrid incompatibility (as well as other dynamics), even in a constant environment. We attribute this difference to the multilocus epistasis and abundant antagonistic pleiotropy present in the gene networks model, but not in the linear regulatory pathways of the Johnson and Porter model. Johnson and Porter (2007) did include a single point of pleiotropic influence; however, much larger neutral networks are apparently required for the effects we describe to emerge. The gene networks model also easily generates hybrid incompatibility under directional selection, but that is not the focus of this article.

To summarize our results, we find that the gene networks model does not in general obey the predictions of Orr (1995). The hybrid fitness load L does not show an initial increase proportional to _K_2, and, in most cases, does not asymptotically approach one, as predicted by equations (1 and 2), above. Instead, depending on the several parameters of the model, under constant selection, we observe three dynamics: (1) hybrid incompatibility between two allopatric populations may not increase at all; or (2) it may increase to large values; or (3) a pair of populations may “drift” in and out of compatibility.

The first of these occurs with certain model parameter values for which a constant environment produces a basin of attraction toward an optimal genotypic state. In this case, a constant environment is both stabilizing and prevents the development of hybrid incompatibility. This affirms an observation of Johnson and Porter, but conflicts with Orr's predictions.

The second dynamic, accumulation of hybrid incompatibility to high levels, is observed when punctuated equilibria (Bergman and Feldman 2003) occur over an extended period despite a constant environment. During this process, the mean fitness intermittently increases as the mean phenotype jumps toward, but does not reach, a static optimum. The third dynamic, drift of a pair of populations in and out of compatibility, occurs in a constant environment when selection is truly stabilizing. In this case, the mean phenotype and mean fitness do not change, but “developmental system drift” (True and Haag 2001) occurs, permitting hybrid incompatibility to arise. (This is the case that Johnson and Porter did not observe in their model.)

To make the distinction clear: a constant environment does not necessarily imply stabilizing selection. If only deleterious variation is present, we say that constant selection is stabilizing: it removes the deleterious variation, selecting for the current, relatively adapted phenotype. However, in the same constant environment, if adaptive variation occurs rarely, then periods of stasis will be interrupted by periods of adaptation. Pleiotropy and complex epistasis can permit the second dynamic to occur. Furthermore, antagonistic pleiotropy increases the size of the near-neutral networks in the genetic neighborhood of local fitness optima, permitting the third dynamic to occur by developmental system drift.

We also extend the model of Orr to the case that the number of interacting loci, Λ, is finite, in which case the growth of hybrid incompatibilities may be curtailed.

PARALLEL EVOLUTION

We define parallel evolution as the independent evolution of similar characters in allopatric populations sharing a common ancestral state. Parallel evolution may be genetic or phenotypic; where relevant, the distinction is indicated below.

A central hypothesis of this article is that parallel evolution can prevent the development of hybrid compatibility. Parallel genetic evolution (absence of divergent substitutions in two populations) acts directly to prevent the development of BDMIs. Parallel phenotypic evolution (parallel evolution in the external phenotype) may or may not entail parallel genetic evolution. If it does (see Fig. 4, for example), then the growth of hybrid incompatibility may be curtailed. Parallel phenotypic evolution may be caused by directional selection, which can strongly drive two allopatric populations in the same evolutionary direction. Alternatively, even if two populations are evolving neutrally, variational constraints (the topology of the neutral network, in the neutral case) dictate the evolutionary pathways along which they may drift. Furthermore, even when there is no fitness increase in individual phenotypes, there can still be indirect selection for mutational robustness, such that certain parts of the neutral network will tend to be more populated (van Nimwegen et al. 1999), resulting in parallel genetic evolution, and curtailing hybrid incompatibility.

4

Parallel genetic evolution in two allopatric populations. (Parameters used: _a_= 1, integer target, _c_= 1.0, μ= 0.01, _N_= 100, _Td_= 1000. Truncation selection was applied.)

Methods

THE PHENOTYPE AND THE GENOTYPE

In the gene networks model of Wagner (1996) and Siegal and Bergman (2002), the phenotype of an individual, shown on the right of Figure 1, is an _M_-dimensional vector representing the amounts of M gene products it produces. The heritable genotype of the individual is represented by an M_×_M matrix, shown on the left of Figure 1. Element (i, j) of the matrix represents the (positive, negative, or zero) regulatory action of gene product j upon the expression of gene i.

1

The “developmental process” of the gene networks model. The M_×_M genotype matrix ∥_gij_∥ and an initial vector are provided. The individual's _M_-dimensional phenotype vector (elements _p_1 through pM) is produced by repeated iterations, each consisting of a matrix multiplication, followed by normalization of the elements of the resulting vector. After Siegal and Bergman (2002).

Initially, a random fraction c (the parameter c is fixed for any given experiment) of the matrix elements is drawn from the standard normal distribution; the rest (1 −c) are set to zero, which permits some control of the amount of epistasis present in the network.

SEGREGATION BY ROWS, AND MUTATION

Upon mating, an offspring is generated, row-by-row, by copying each new row from the matrix of one parent or the other with equal probability, akin to segregation of chromosomes (rows), with no recombination within chromosomes. Then mutation is applied at rate μ per genome and entails replacement of a random nonzero matrix element with another number from the standard normal distribution. Thus the per (nonzero) element mutation rate will be μ/(_cM_2).

THE DEVELOPMENTAL PROCESS

At the beginning of the developmental process, each individual's phenotype is set to a fixed initial vector, the elements of which are in the closed interval [0,1]. Each iteration of the developmental process, illustrated in Figure 1, consists of two parts: first, the genotype matrix and the current phenotype vector are multiplied; and second, the entries of the resulting vector are normalized. The following sigmoid normalization function is used to map each element x of the product back to the interval [0, 1]:

3

The vector of normalized elements becomes the phenotype vector for the next iteration of the developmental process. The parameter a (fixed for any given experiment) controls the steepness of the sigmoid function, which is plotted for two values of a in Figure 2.

2

The sigmoid normalization function sig(x, a) from equation (3). Each element of the product of the genotype matrix and the current phenotype vector is normalized back to the interval [0, 1] with this function, producing the corresponding element of the next iteration of the phenotype vector.

A high value of a (e.g., 10) creates a tendency for the elements of the phenotype vector to assume near-zero or near-one values. A value near 0.5 at an a of 10 requires a precise balance of regulatory effects. Low values of a (e.g., 1) permit elements of the phenotype vector to more easily assume intermediate values and make it difficult to attain near-zero or near-one values.

We restrict the phenotypic elements to the interval [0, 1], rather than the interval [−1, 1] used by Wagner (1996) and Siegal and Bergman (2002), because we interpret the value of each phenotypic element as the presence (positive value) or absence (zero value) of the corresponding gene product. Positive genotype matrix entries imply upregulation, and negative entries downregulation. Therefore, a mutation is required to reverse the direction of regulation; i.e., there is no critical point about which a slight change in the value of a phenotypic element will suddenly reverse the direction of regulation.

INTEGER- AND REAL-VALUED PHENOTYPIC TARGETS

An additional model parameter is the “target type.” The elements of the initial phenotype vector, and of the target phenotype vector, may either be “integer-valued” (chosen with equal probability to be zero or one), or “real-valued” (chosen uniformly from the real numbers in the closed interval [0, 1]).The initial phenotype vector in a given experiment will have the same type as the target phenotype vector. A ten-dimensional integer-valued target is plotted as the 10 large orange circles in both panels of the first row of Figure 3; the vector represented is (1, 0, 1, 0, 1, 1, 0, 0, 0, 1). A real-valued target is plotted as the 10 large orange circles in each panel of the second row. (The other features of this figure are discussed in detail below.)

3

Phenotypic plots at generation 0 (left column) and generation 100 (right column) for four runs using four combinations of a values (_a_= 1 or _a_= 10) and target types (integer or real). The 10 large orange circles in each plot represent the 10-dimensional target vector. Overlaid on the orange circles are 1000 smaller colored dots, indicating the 10-dimensional phenotype vectors of the 100 individuals in the population. The phenotype of one individual is plotted as 10 dots, one in each column, all of the same color. This color represents the distance d (see eq. 4) from that individual's phenotype to the target phenotype: dark blue indicates _d_= 1.0; dark red indicates _d_= 0.0; and an intermediate green color indicates _d_= 0.5. It is not generally possible to distinguish the 10 dots corresponding to a specific individual's phenotype vector; rather, the plots give an overall picture of the distribution in the population of the values of each element, in comparison to the target. (Other parameters used: _c_= 0.75, μ= 0.1, _N_= 100.)

Because high values of a tend to produce phenotypes with near-one or near-zero elements, experiments with real-valued targets and a high value of a tend to have difficulty producing a phenotype near the target phenotype, but they easily approximate integer-valued targets. In contrast, if _a_= 1, intermediate phenotype values are easily attained, but extreme matrix values are required to produce phenotypes with near-zero or near-one elements.

VIABILITY, FITNESS, AND MATING

At the beginning of the developmental process, each individual's phenotype is set to the fixed initial vector uinit. If the individual's phenotype converges to any vector in less than 100 developmental iterations (rather than entering a cycle, or not converging), then that individual is declared to be viable; otherwise, it is removed from the population. Convergence is determined as follows. Let d be the distance metric between two phenotype vectors ua and ub defined by:

4

At the end of each iteration, an average phenotypic vector uavg over the last τ= 10 iterations is computed. If the average d between the last τ phenotype vectors and uavg is below the threshold value ɛ= 10−4, then convergence has occurred, and the individual is declared viable.

At the beginning of each experiment, an _M_-dimensional phenotypic “target vector” is randomly chosen; this target defines the best possible phenotype. For all viable offspring, we compute a probability of survival to maturity (“fitness”), f, using the distance between the final, converged phenotype and the target vector:

5

The parameter σ permits control of the severity of selection; we use a value of σ= 0.1 here.

At each generation, we repeat the following operations until there are exactly N surviving individuals: random mating with segregation of chromosomes; mutation; viability selection (i.e., whether the phenotype converges); fitness selection by equation (5).

POPULATION METRICS

We define I, the incompatibility between two allopatric populations, as:

6

where _W_1 and _W_2 are the average fitnesses within the two allopatric populations, and H is the average fitness (counting nonviables as having zero fitness) of the hybrid offspring produced by mating pairs with one parent coming from each population. This is intended to be analogous to Orr's hybrid fitness load L in equations (1 and 2); thus Orr's predictions should apply to I.

Define the “self-incompatibility,”Is, within a single population as:

7

where Woffspring is the average fitness of offspring of the current generation (nonviables have zero fitness), and W is the average fitness of the current generation. Typically in our models, Is is positive immediately after the founding of a population, but it quickly decreases to near-zero.

Gavrilets (2004) noted two stages in the evolutionary dynamics of a simple two-locus BDM model: the first, driven by strong selection, during which most of the population rapidly moves onto a neutral fitness ridge, and a second in which the dynamics are driven by weak forces such as mutation. The first is analogous to the rapid elimination of Is, which we routinely observe. Nei et al. (1983) noted a similar effect in their models.

We measure Ki, the number of substitutions in population i since the population split, as follows. At the moment of population split, each Ki is set to zero. The frequencies of all polymorphisms at each matrix element of the N genotypes in each population are tracked. When the frequency of a particular allele first exceeds 0.95 in population i, Ki is incremented by one. That allele is then marked as fixed unless its frequency subsequently decreases below 0.5. In comparing two of the allopatric populations, i and j, we sum the number of substitutions between them: _K_=Ki+Kj.

Results

All experimental populations described in this article were created from six founder individuals that were randomly generated and subjected to one round of viability and fitness selection. The populations were thereafter maintained at the constant size of _N_= 100 (Figs. 3 and 4) or _N_= 250 (all subsequent Figures involving simulation of gene networks) individuals. For Figures 3–6 and 8, animated versions of the figures are available in the online Supporting information.

5

Incompatibility, I, for the 10 replicates r1, r2, …, r10 as a function of _K_2, for a = 1. Left: integer-valued target; right: real-valued target. (_c_= 1.0, μ= 0.003, _N_= 250, _Td_= 1000).

6

Incompatibility, I, for the 10 replicates r1, r2, …, r10 as a function of _K_2, for _a_= 10. Left: integer-valued target; right: real-valued target. (_c_= 1.0, μ= 0.003, _N_= 250, _Td_= 1000).

8

Incompatibility, I, for the 10 replicates r1, r2, …, r10 as a function of _K_2, for _a_= 10 and real target, _c_= 0.1 (μ= 0.003, _N_= 250, _Td_= 1000).

A VIEW OF THE PHENOTYPE FOR FOUR-PARAMETER SETS

Figure 3 illustrates four experimental runs (four rows) of a single population at two points in time: generation 0 (left column) and generation 100 (right column). The four runs correspond to the four combinations of two values of the parameter a (_a_= 1, or _a_= 10) and two types of phenotypic target (integer-valued or real-valued). The left plot in the first row shows this run at generation 0. Because a is low (_a_= 1), the values of the elements of the phenotype vectors are well distributed through the intermediate values between zero and one. The dots are generally colored blue and blue-green, indicating that no phenotypes very near the target have yet evolved. In contrast, the right side of the first row of the figure shows generation 100 of the same run: the elements of the phenotypes are clustered closer to their target values, and the orange and red colors of the dots reflect the evolved closeness of each phenotype vector to the target vector.

The second row of Figure 3 shows a run with _a_= 1 and a real-valued target. At generation 0, the dots representing the phenotypes are more green than blue. Because the elements of the target vector (orange circles) have intermediate rather than extreme values, the initial phenotype vectors are in general closer to the target from the outset. At generation 100, the phenotype vectors again cluster near the target vector, as indicated by the orange–red colors of the dots.

The third row of Figure 3 shows a run with _a_= 10 and an integer-valued target. The initial distribution (left) of the elements of the phenotype vectors is biased toward the extremes of 0 and 1 due to the sharpness of the sigmoid normalization curve (see Fig. 2), which makes it difficult to obtain intermediate element values. At generation 100, the elements of the phenotypes are clustered about the elements of the target, as reflected by the dark red color of the dots. The high value of a causes a bias toward the production of phenotypes with extreme values, causing the population to evolve very fit phenotypes in this case.

The fourth row of Figure 3 has _a_= 10 and a real-valued target. The initial distribution of phenotypes is similar to that of the third row with elements biased toward extreme values (0 or 1). However, in this case, the target has intermediate values, so the initial phenotypes are far from it, as reflected by the blue shades of the dots. By generation 100, the population is still unable to evolve a very fit phenotype: the colors of the dots are an intermediate green, and, in most cases, the phenotypic elements are not clustered near the target elements.

We will use these four combinations of a and the target type as our primary means of manipulating the evolutionary challenge faced by our populations. Later, we will also discuss the effect of varying the c parameter.

PARALLEL EVOLUTION IN THE GENOTYPE

Each experiment below starts with a single population. At generation Td (“time of divergence”), the population is cloned into two populations, between which there is no subsequent migration. We typically wait to split the population until it has become well adapted (attained near-optimal fitness), except when this is not possible due to very slow adaptation, as described below. This is also sufficient time to eliminate the initial Is that is due to the founding of the population with randomly generated individuals.

To generate a visually striking example of parallel genetic evolution, we increase the selective intensity by applying truncation selection. Specifically, in addition to our normal selective regime, the bottom 25% of individuals by fitness rank were not permitted to mate. Genotypes, rather than phenotypes, are plotted in Figure 4. At generation 0 (first column), the experiment starts with a single population. The 10 × 10 grid of colored squares represents the genotype matrix of one individual randomly selected from the population of 100 individuals. The darkest reds represent high-valued matrix elements (+3 and higher), and the darkest blues represent low-valued matrix elements (−3 and lower). Zero is represented by an intermediate green color. At generation 1000 (second column), the population is cloned, forming two identical populations. The genotype of one randomly selected individual is displayed from the original Population 1 (top), and from the cloned Population 2 (bottom), respectively. There is no further migration between the two populations for the rest of the experiment.

By generation 10,000 (third column), it is visually apparent that parallel genetic evolution has occurred. The target vector for this run was the integer-valued target (i.e., the orange circles) shown in both panels in the first row of Figure 3. The patterns of red and blue rows of both of the genotype matrices in the last column of Figure 4 correspond exactly to the pattern of high and low elements of the target vector, i.e., (1, 0, 1, 0, 1, 1, 0, 0, 0, 1): high values in row i of the genotype matrix upregulate gene product i, and low values do not. The standard fitness model (without truncation selection), used in the rest of our results, also produces parallel evolution, but not to such a visually striking degree as in Figure 4.

THE BIOLOGICAL INTERPRETATION OF THE PARAMETER “_a_”

We interpret the parameter a as specifying the sensitivity of the regulatory response to the intermediate phenotype vectors that are computed during development. Do the genes act additively, at least in the middle of the range of their minimum and maximum combined effect (low a)? Or do they behave more like threshold characters, for which there is a sudden qualitative change when a threshold is passed (high a)?

When _a_= 1, a broad basin of attraction is produced around local phenotypic optima; the population can quickly “hill climb” the fitness gradient, which is nonzero nearly everywhere, toward local optima, and to the global optimum. At the same time, _a_= 1 makes it difficult to produce extreme (near-zero or one) phenotypic values. There is, therefore, a distinct fitness peak with a very small neutral network at the top. In terms of a metaphorical “fitness landscape,”_a_= 1 corresponds to broad-based, sloping hills that are peaked, rather than flat, on top.

When _a_= 10, it is much easier to produce extreme phenotypic element values, because for production of such extremes it is not necessary that multiple regulatory pathways act in the same direction. Therefore much larger neutral networks are created about the optima. At the same time, the basins of attraction about optima become much narrower, because for most of the input range of the sigmoid function, the gradient is near-zero. For the metaphorical “fitness landscape,” there is a low plateau containing mesas: narrow elevated areas with steep sides and a flat top.

The parameter a could also be placed under evolutionary control, an analysis not pursued here.

I VERSUS _K_2

The extent of hybrid incompatibility was measured by our metric I, for the four combinations of a and target type illustrated in Figure 3. We performed 10 replicate runs of eight parameter sets created using: two values of a, the steepness of the sigmoid curve (1, and 10); two target types (integer-valued, and real-valued); and two values of c (0.1, and 1.0). Following Johnson and Porter (2000), we used a population of _N_= 250 and a genomic mutation rate of μ= 0.003. The initial population was split into two allopatric populations of size 250 at time _Td_= 1000.

All 10 replicates of each parameter set were started with the same, fixed random seed. Thus, for each repetition of the same parameter set, each detail of the simulation (the phenotypic targets, the founder individuals, and all other individuals in the population at each generation) was identical, up to the moment of population split. At the moment of split, a new random seed was set, depending on the repetition number. This scheme permits us to eliminate variation between runs due to random events occurring before the population split. We further replicated all of the above (the 10 replicates of eight parameter sets) with 10 different initial random seeds (yielding a total of 800 runs), the choice of initial random seed had no significant effect. One set of 80 runs is shown below, but all sets are available in the Supporting information.

We chose a Td of 1000, because this delay allows all combinations of a and the target type (with the exception of _a_= 10 and real-valued target) to attain high fitness before the population is split; it also allows the self-incompatibility Is to be largely eliminated. All experiments ran for _Td_+ 30,000 generations.

THE _a_= 1 CASE: PARALLEL EVOLUTION LIMITS THE INCREASE OF INCOMPATIBILITY

When _a_= 1, there is a single evolutionary dynamic. Whether the target is real- or integer-valued, the population can quickly “hill climb” to the global optimum where there is little neutral drift. Instead strong stabilizing selection keeps the population near the optimum genotype, and pairs of allopatric populations are prevented from developing hybrid incompatibility.

On the left of Figure 5, we plot I against _K_2 for the 10 runs with _a_= 1, an integer-valued target. Clearly I does not increase as the square of K, as would be predicted by Orr. The range of the I values is very small: by the end of each experiment (corresponding to 31,000 generations, or _K_2∼ 20,000–35,000 in the figure), I is less than 0.001, indicating that the populations are extremely compatible despite many generations of allopatric evolution. Speciation, which would require values of I approaching 1.0, does not occur. We expect strong parallel evolution for this parameter set, since, for _a_= 1, extreme matrix values are required to produce near-zero- and near-one-valued phenotypic elements to match the integer-valued target elements. The average fitness (see online Supporting information) hovers around 0.98–0.99.

On the right of Figure 5, a similar plot is shown for the case of _a_= 1 and real-valued target. Significant parallel evolution is still expected for this parameter set (because even a real-valued target has some extreme values, driving some parallel evolution); and incompatibility remains nearly as low as for the first parameter set: less than one percent by the end of each run. The average fitness (see Supporting information) quickly attains 0.98–0.99 and remains at about that level on average.

THE _a_= 10 CASE: TWO DIFFERENT MECHANISMS PERMIT HYBRID INCOMPATIBILITY

For _a_= 10, two distinct evolutionary dynamics occur: (1) For integer targets, the variational bias toward extreme output values allows the population to quickly find the tops of the metaphorical “mesas” described above; the population then drifts about these large neutral networks. (2) In the case of real targets, variation when _a_= 10 is so biased toward producing extreme phenotypic values that the population has great difficulty finding local optima; a dynamic of punctuated equilibria ensues, as occasional adaptation occurs, followed by long periods of stasis. Near-optimal fitness is never attained during the experiments. Under this dynamic, two allopatric populations may evolve asynchronously (one ahead of the other), or they may optimize the elements of the phenotype in differing order; either type of asynchrony allows the temporary development of hybrid incompatibility, which may later be reduced as the second population “catches up” to the first.

True and Haag (2001) describe a process they name “Developmental System Drift” (DSD), which means change in the developmental system (implying genetic change) even without phenotypic change. True and Haag state that because DSD happens by chance alone, rather than by selection, “drift” is an appropriate term. DSD does occur in the gene networks model.

We suggest, however, that evolution of the developmental system can also be caused by indirect selection for robustness (van Nimwegen et al. 1999). Such evolution can be considered “neutral” in the sense that all (viable) phenotypes produced may be of identical fitness (probability of survival to adulthood) although, depending on the topology of the neutral network, certain genotypes may experience greater mutational load (producing a higher proportion of inviable mutants), resulting in an indirect form of selection.

Furthermore, evolution of the developmental system is also important during episodes of adaptation, not only during movement of the population on a neutral ridge. Consider a population split in two with the two allopatric populations subject to the same directional selection. If we observe identical fitness increase and parallel phenotypic evolution, this does not imply that the populations have also experienced parallel genetic evolution. If the two populations do not differ genetically, even if this difference is cryptic (not expressed in any phenotypic difference), then again there is the potential for hybrid incompatibility to arise. We have observed this effect in gene networks simulations with time-varying phenotypic targets (data not shown). Hence, we prefer to use the term “Developmental System Evolution” (DSE) rather than True and Haag's DSD. We may then speak of “DSE under neutrality” or “DSE under phenotypic constancy” to indicate constancy of fitness, or phenotype, respectively, or “divergent DSE despite parallel phenotypic evolution” to indicate that two populations whose phenotypes have evolved in parallel are nonetheless genetically distinct.

On the left side of Figure 6, I is plotted against _K_2 for the case of _a_= 10 and integer-valued target. Here I does increase to significant values (e.g., 20–40%) during several runs. However, increase in I often seems to be temporary, as I wanders to a peak, and then falls to low values (less than 5% incompatibility).

Fitness (see Supporting information) is already near-optimal (>0.997) before the population splits, and remains extremely high. However, there is a neutral network about the global optimum, and drift about this neutral network by DSE allows the development of hybrid incompatibility, which can then fade away as pairs of populations drift back into compatibility. Hence, this mechanism of generation of incompatibility is potentially reversible; whether it reverses depends on the size and topology of the neutral network at the “top of the mesa.”Gavrilets (2004, p. 174) suggests that hybrid incompatibility should increase with genetic distance. (This relationship has been found in _Drosophila_[Coyne and Orr 1989, 1997], frogs [Sasa et al. 1998], birds [Lijtmaer et al. 2003], and Lepidoptera [Presgraves 2002], although it was not detected in darters [Mendelson 2003] and was found in only one of three genera of angiosperms [Moyle et al. 2004].) If genetic distance were to increase over time, as would be expected for very large neutral networks, then hybrid incompatibility should increase over time by drift. In contrast, we believe that the dynamic of “drifting” in and out of compatibility, as illustrated on the left side of Figure 6, is best explained in terms of random walks on neutral fitness ridges that are relatively confined, and where it cannot be assumed that genetic distance will tend to increase by drift. Rather, genetic distance may wax and wane as two populations randomly move relative to one another about a confined set of neutral ridges; hybrid incompatibility (as a roughly monotonic function of genetic distance) may therefore also wax and wane.

On the right side of Figure 6, I is plotted against _K_2 for _a_= 10 with a real-valued target vector. Here, I does increase, in a volatile manner, to significant levels: 30–80%, depending on the run. Plotting the average fitness, W, against _K_2 (see Supporting information), this volatility is clearly due to punctuated equilibrium effects. The average fitnesses of all populations show clear, persistent plateaus, punctuated when an adaptive variant is discovered by sudden jumps up to a higher plateau, which then persists for a time. In this case, despite a constant environment, the selective pressures are not constant: long periods of stasis due to stabilizing selection persist when adaptive variation is absent, but a new adaptive variant may sweep through its population, potentially causing a cascade of fixations that can contribute to incompatibility.

This incompatibility is also potentially reversible in an interesting way: suppose that two populations optimize the several elements of the phenotype one at a time (slowly, via punctuated equilibria), and in differing orders. Further suppose that after a very long time, they ultimately optimize all phenotypic elements. The intermediate state could include many possible phenotypes (e.g., there are many permutations of five of 10 “correct,” i.e., near-optimal elements, and a correspondingly large number of (potentially incompatible) genotypes that could produce these phenotypes. In contrast, there is only one way to get all 10 elements “correct,” and fewer genetic alternatives to produce this phenotype. Therefore one could expect an initial increase in incompatibility, followed by a decline of incompatibility, despite allopatric conditions. (See also (Orr 2005) regarding the probability of parallel evolution.)

This situation (continuing rare adaptation, even over extended periods) may be common in biological systems in a constant environment. For example, in populations of E. coli that were maintained in a constant environment for more than 20,000 generations, adaptation continued throughout the entire period (Cooper and Lenski 2000). Parallel genetic and phenotypic evolution among the several populations was directly observed (Cooper et al. 2003; Philippe et al. 2007).

Significant incompatibilities will be frequent in such situations of evolutionary asynchronicity. Compatibility may, or may not, be attained by parallel evolution in the very long term. This will depend on the structure of the network of possible evolutionary pathways, in particular whether the pathways ultimately converge, or lead to distinct (and potentially incompatible) endpoints (see Discussion).

EXTENDING ORR's MODEL TO THE CASE OF FINITE LOCI

Here we present an extension of Orr's (1995) analysis, which considers an infinite number of loci. Our simulations obviously involve a finite number of loci (specifically, _cM_2, or between 10 and 100 nonzero loci for the experiments shown here). In addition, Orr does not allow more than one substitution at a locus, whereas we do. We now extend Orr's (1995) model to the case of finite loci with multiple substitutions permitted per locus and present an algorithmic solution (rather than a closed form solution) for how L, the hybrid load, should increase with K, the number of substitutions.

Let D be the total number of diverged loci between two populations. The number of potentially incompatible pairs is D(_D_− 1)/2. In the case that the number of loci, Λ, is finite, D can increase to at most Λ, and it will grow more slowly than K (the number of substitutions) due to the chance of repeated substitutions at the same locus. Modifying equation (1), we derive as the hybrid load when D loci have diverged. We numerically compute D, and hence L, for the case of finite Λ, as follows. We create an array of Λ integers, all initially set to zero. When the first random substitution occurs at a particular locus, the corresponding array element is set permanently to one: the locus has irreversibly diverged. At any point, we can sum the elements of the array to obtain the total number of diverged loci, D. (This approach neglects the possibility of back mutation, as did Orr's original model. If back mutation, or alternatively, parallel evolution, were permitted then D could decrease) We also introduce an additional complication by partitioning the loci into “incompatibility groups” such that only pairs within a group can be incompatible. Let G be the number of groups, assumed to be of equal size. For example, for the case of Λ= 1000 and _G_= 10, we take 10 groups of 100 loci each, combining their fitness penalties multiplicatively. In Figure 7, we plot L versus K for ; on the left is _G_= 1 (a single group), and on the right is _G_= 10. Ten replicate runs were done for each of Λ= 1,000,000, Λ= 1000, Λ= 100, Λ= 30, and Λ= 10; hence there are 10 (stochastically overlapping) lines plotted for each of the first five colors in each plot. For comparison, Orr's infinitely many loci case, equation (1), is shown as f(K); and our correction for loci grouped into _G_“incompatibility groups” in the limit of Λ→∞, is shown as . For _G_= 1, f(K) =g(K). For both values of G, the plots for Λ= 1,000,000 (red) are very near the infinitely many loci case, the black dashed line representing g(K).

7

Growth of hybrid load L with the number of substitutions K in the case of finite Λ, the number of loci. Left: one “incompatibility group”; right: 10 groups. _r_= 0.001. f(K) is Orr's infinite loci analysis and g(K) is our correction for multiple incompatibility groups. On the left, for a single group, g(K) is equal to f(K). In both plots, g(K) is very close to the red lines for Λ= 1,000,000.

On the left of Figure 7 the plots for Λ= 100 (blue) approach complete speciation (_L_= 1.0) more slowly than Orr's infinitely many loci case (f(K), black line), due to the possibility of multiple substitutions at the same locus. More dramatically, the plots for Λ= 30 (pink) and Λ= 10 (cyan) on the left of the Figure do not approach 1.0 at all, but plateau around values of 0.35, and 0.04, respectively; L cannot continue to grow after all loci have diverged. We compute that L will asymptotically approach as K approaches infinity and all loci diverge. Thus, for example, , matching the pink lines on the left of the figure. On the right side of the figure, for example, matches the blue lines.

Clearly, finitely many loci can strongly limit the growth of incompatibility, as on the left of Figure 7; this is also reflected in our gene networks simulations, below. In addition, grouping of loci in “incompatibility groups” can limit incompatibility, or substantially delay its onset, even for high values of Λ, as on the right of Figure 7. Additional volatility, and nonmonotonicity of L, can be generated if other distributions are used for r (see Orr 1995); however, to model these cases, the individual values of r drawn for each incompatible pair must be tracked (not shown).

REDUCING THE NUMBER OF INTERACTING LOCI LIMITS INCOMPATIBILITY

Returning to our gene networks simulations, we indeed find that when Λ, the number of loci, is reduced, incompatibility is restricted. Reduction of Λ is effected in our simulations by reducing c. (Recall that a fraction 1 −c of the matrix elements of the founding individuals are set to zero, and that mutation does not alter zero-valued matrix elements.) In Figure 8, we plot I against _K_2 for the case of _c_= 0.1 with _a_= 10 and a real-valued target. For this case, parallel genetic evolution did not occur, because the final element values differed significantly between the two populations (see Supporting information). Despite the lack of parallel evolution, however, all pairs of populations maintain incompatibilities of less than one percent at the end of the run. Compare Figure 8 with the right side of Figure 6, for which _c_= 1.0 and incompatibilities of 30–80% were reached. We may conclude that in this case incompatibility was limited not by parallel evolution, but because the number of interacting loci was small. Therefore, in general, either mechanism may prevent incompatibility from increasing.

Discussion

ANOTHER MODEL OF INCOMPATIBILITY

Nei et al. (1983) performed computer simulations of hybrid incompatibility in one- and two-locus multi-allele models. In the one-locus model, multiple alleles (labeled A−2, A−1, A0, A1, A2, etc.) are permitted at the locus. Homozygotes are viable, as are heterozygotes whose two alleles are only one step apart (i.e., an A_i_ A_i_− 1 heterozygote or A_i_ A_i_+ 1 heterozygote are viable), but any other heterozygote is inviable. This is not the standard formulation of BDM incompatibility, which considers incompatibility due to epistasis between (at least) two loci. Nei et al.'s two-locus case considers two independent loci of the same description, with no epistatic interaction between them. In both the one- and two-locus cases, the low fitness of heterozygotes more than one step apart greatly retards the fixation of novel alleles, producing a fundamentally different dynamic from the independent fixations of Orr's model. Nei et al. also observed an apparent delay, then a rapid acceleration of hybrid incompatibility, qualitatively similar to Orr's “snowball.” However, in the Nei model, this dynamic is due to a delay and then rapid fixation of a novel allele, rather than to a growing number of potential incompatibilities, as in the Orr case.

DYNAMICS OF INCOMPATIBILITY ON NETWORKS OF EVOLUTIONARY PATHWAYS

The “holey landscapes” metaphor of Gavrilets (2004) describes high-dimensional genotype spaces as pervaded by “fitness ridges,” nearly neutral networks on which populations may drift, constrained to remain on a ridge by stabilizing selection alone. Gavrilets' original formulation did not include sexual recombination: the spatial structure of the landscapes—the allowable stepwise movements of genotypes on them—was defined in terms of mutational transitions only. Therefore, whereas an asexual population would remain confined to a “connected component” of genotypes by the action of mutation, introducing sexual recombination would break this confinement. Gavrilets extended this model to sexual populations by neglecting polymorphism (Gavrilets and Gravner 1997). In this case, a population can be represented by its single extant genotype, and transitions from one state to the next represent allele fixations in the population. This assumes only one mutation in the entire population for the duration of each fixation event, to ensure that the population remains confined to its current connected component.

We define the network of “evolutionary pathways” that population replicates may traverse as the union of the neutral networks spanning genotype (monomorphic population state) space at various fitness levels, together with the rare adaptive transitions connecting them. The evolutionary outcome depends stochastically upon the topology of this network. In this context, we may reinterpret Orr's (1995) result as follows. Although it is true that two populations proceeding along divergent pathways will tend to have increasing genetic distance between them (and accelerating hybrid incompatibility, by equation [2]), genetic distance nevertheless may decrease for two reasons: (1) under neutrality, a population need not move in a consistent direction along a pathway; and (2) parallel genetic evolution may occur, with two pathways converging to an identical genotype by different routes (e.g., the same fixations but in different order), or asynchronously. A set of pathways diverging from the same point may in general have different probabilities of being evolutionarily “chosen” by selection or drift. Orr (2005) discusses how this increases the probability of parallel evolution.

Because genetic distance may decrease, we replace K, the number of substitutions (which can never decrease), with D, the genetic distance measured by the number of divergent alleles (which may decrease in certain circumstances). Equation (1) then becomes:

8

In addition, should not be assumed to be constant between different pairs of pathways. Orr's result, reinterpreted via equation (8), does not predict an inevitable increase of hybrid load, although, if two populations do diverge genetically, increasing (indeed, initially accelerating) hybrid load is expected. However, the rate of this increase will depend on the particular between the two pathways along which the two populations are moving.

THE EVOLUTION OF GENE EXPRESSION

A classic paper by King and Wilson (1975) is widely cited for popularizing the notion that much phenotype divergence between species may be due to changes in the regulation of gene expression, as distinct from changes in protein-coding DNA. Although the role of structural mutations clearly remains important in evolution (Hoekstra and Coyne 2007), it has been shown in yeast (Borneman et al. 2007) and in primates (Gilad et al. 2006b) that gene expression levels diverge more strongly between closely related species than their nucleotide sequences. For example, between humans and chimpanzees, it was found that 10% of genes differed in expression in at least one brain region (Khaitovich et al. 2004).

The role of natural selection in the evolution of gene regulation is hotly debated. Many researchers favor the idea that the evolution of gene expression is primarily driven by selection (Nuzhdin et al. 2004; Gilad et al. 2006a,b; Haygood et al. 2007; Mustonen and Lassig 2007; Hutter et al. 2008), with some studies advocating stabilizing selection or directional selection as more “important.” However, a significant minority urge that the role of neutral evolution must not be neglected (Lynch 2007; Fay and Wittkopp 2008). It is in the context of this debate that we present our results on the simulated evolution of gene networks in a constant environment.

We have focused on constant environments in this article; however, our gene networks also show rapid growth of hybrid incompatibility in the presence of (1) divergent selection (selection toward different optima in the two populations), and (2) directional, nondivergent selection (selection toward the same optimum in both populations). That incompatibility would evolve under divergent selection is intuitive: if two populations evolve toward two different optima, it is not surprising that a hybrid might be less fit in both environments. That incompatibility increases under directional, nondivergent selection is less intuitive; however, it can be understood if one considers that there are multiple genetic means to the same end (This is “divergent DSE despite parallel phenotypic evolution,” as we described above.), and that selected phenotypic change may drive a faster pace of genetic change than would proceed by drift.

Although speciation is rapid under both divergent and nondivergent directional selection, we have here focused on constant environments, demonstrating several dynamics, including strong growth of incompatibility, and waxing and waning incompatibility. Our conclusion is that the role of neutral evolution in the growth of incompatibility in biological gene regulatory networks, and hence in speciation, should not be dismissed out of hand.

Conclusions

We have extended Orr's model in two ways to account for the complex dynamics of incompatibility in the gene networks model: (1) by considering genetic distance, (which may decrease by back mutation or parallel genetic evolution) rather than number of substitutions (which does not decrease), and (2) by considering a finite, rather than an infinite, number of interacting loci.

We have shown that directional selection is not required to produce hybrid incompatibility in a model inspired by biological gene networks. Rather, the growth of hybrid incompatibility may proceed by any of several dynamics in a constant environment, including asynchronous parallel evolution, and developmental system evolution with no change in phenotype. The important point is that many processes can cause genetic change, some of which can increase hybrid incompatibility. Divergent directional selection clearly will promote phenotypic divergence between two populations, which necessarily implies genetic divergence between them. Nondivergent directional selection will promote phenotypic change, but not necessarily phenotypic divergence; hence it may or may not result in genetic divergence. Likewise, neutral phenotypic change may or may not result in genetic divergence. Our models suggest that two mechanisms, (1) phenotypically neutral changes in gene regulation, and (2) asynchronous parallel evolution due to rarity of adaptive mutation, may also lead to speciation in biological populations.

Associate Editor: J. Feder

ACKNOWLEDGMENTS

The authors thank M. Desai, L. Lehmann, J. Masel, and S. Ramachandran for important discussions and suggestions. We are also indebted to A. Porter (in particular for the suggestion that incompatibility could be limited in the finite loci case) and an anonymous reviewer.

LITERATURE CITED

Bateson

,

W.

1909

. Heredity and variation in modern lights. Pp.

85

–

101

in

A. C.

Seward

, ed.

Darwin and modern science: essays in commemoration of the centenary of the birth of Charles Darwin and of the fiftieth anniversary of the publication of the Origin of species

.

Cambridge Univ. Press

, Cambridge.

Bergman

,

A.

, and

W. M.

Feldman

.

2003

. On the population genetics of punctuation. Pp.

81

–

100

in

J. P.

Crutchfield

, ed.

Evolutionary dynamics: exploring the interplay of selection, accident, neutrality, and function

.

Oxford Univ. Press

, New York.

Borneman

,

A. R.

,

T. A.

Gianoulis

,

Z. D. D.

Zhang

,

H. Y.

Yu

,

J.

Rozowsky

,

M. R.

Seringhaus

,

L. Y.

Wang

,

M.

Gerstein

, and

M.

Snyder

.

2007

.

Divergence of transcription factor binding sites across related yeast species

.

Science

317

:

815

.

Cooper

,

T. F.

,

D. E.

Rozen

, and

R. E.

Lenski

.

2003

.

Parallel changes in gene expression after 20,000 generations of evolution in Escherichia coli

.

Proc. Natl. Acad. Sci. USA

100

:

1072

.

Cooper

,

V. S.

, and

R. E.

Lenski

.

2000

.

The population genetics of ecological specialization in evolving Escherichia coli populations

.

Nature

407

:

736

.

Coyne

,

J. A.

, and

H. A.

Orr

.

1989

.

Patterns of speciation in Drosophila

.

Evolution

43

:

362

.

Coyne

,

J. A.

, and

H. A.

Orr

.

1997

.

Patterns of speciation in Drosophila revisited

.

Evolution

51

:

295

.

Darwin

,

C.

1859

.

The origin of species by means of natural selection

.

John Murray

, London.

Dobzhansky

,

T.

1936

.

Studies on hybrid sterility. II. Localization of sterility factors in Drosophila pseudoobscura hybrids

.

Genetics

21

:

113

–

135

.

Fay

,

J. C.

, and

P. J.

Wittkopp

.

2008

.

Evaluating the role of natural selection in the evolution of gene regulation

.

Heredity

100

:

191

.

Gavrilets

,

S.

2004

.

Fitness landscapes and the origin of species

.

Princeton Univ. Press

, Princeton, NJ.

Gavrilets

,

S.

, and

J.

Gravner

.

1997

.

Percolation on the fitness hypercube and the evolution of reproductive isolation

.

J. Theoretical Biol.

184

:

51

.

Gilad

,

Y.

,

A.

Oshlack

, and

S. A.

Rifkin

.

2006a

.

Natural selection on gene expression

.

Trends Genet.

22

:

456

.

Gilad

,

Y.

,

A.

Oshlack

,

G. K.

Smyth

,

T. P.

Speed

, and

K. P.

White

.

2006b

.

Expression profiling in primates reveals a rapid evolution of human transcription factors

.

Nature

440

:

242

.

Haygood

,

R.

,

O.

Fedrigo

,

B.

Hanson

,

K. D.

Yokoyama

, and

G.

Awray

.

2007

.

Promoter regions of many neural- and nutrition-related genes have experienced positive selection during human evolution

.

Nat. Genet.

39

:

1140

.

Hoekstra

,

H. E.

, and

J. A.

Coyne

.

2007

.

The locus of evolution: evo-devo and the genetics of adaptation

.

Evolution

61

:

995

.

Hutter

,

S.

,

S. S.

Saminadin-Peter

,

W.

Stephan

, and

J.

Parsch

.

2008

.

Gene expression variation in African and European populations of Drosophila melanogaster

.

Genome Biol.

9

:

R12

.

Johnson

,

N. A.

, and

A. H.

Porter

.

2000

.

Rapid speciation via parallel, directional selection on regulatory genetic pathways

.

J. Theor. Biol.

205

:

527

.

Johnson

,

N. A.

, and

A. H.

Porter

.

2007

.

Evolution of branched regulatory genetic pathways: directional selection on pleiotropic loci accelerates developmental system drift

.

Genetica

129

:

57

.

Khaitovich

,

P.

,

B.

Muetzel

,

X.

She

,

M.

Lachmann

,

I.

Hellmann

,

J.

Dietzsch

,

S.

Steigele

,

H.-H.

Do

,

G.

Weiss

,

W.

Enard

, et al.

2004

.

Regional Patterns of Gene Expression in Human and Chimpanzee Brains

.

Genome Res.

14

:

1462

–

1473

.

King

,

M. C.

, and

A. C.

Wilson

.

1975

.

Evolution at two levels in humans and chimpanzees

.

Science

188

:

107

.

Lijtmaer

,

D. A.

,

B.

Mahler

, and

P. L.

Tubaro

.

2003

.

Hybridization and postzygotic isolation patterns in pigeons and doves

.

Evolution

57

:

1411

.

Lynch

,

M.

2007

.

The frailty of adaptive hypotheses for the origins of organismal complexity

.

Proc. Natl. Acad. Sci. USA

104

:

8597

.

Mendelson

,

T. C.

2003

.

Sexual isolation evolves faster than hybrid inviability in a diverse and sexually dimorphic genus of fish (Percidae: Etheostoma)

.

Evolution

57

:

317

.

Moyle

,

L. C.

,

M. S.

Olson

, and

P.

Tiffin

.

2004

.

Patterns of reproductive isolation in three angiosperm genera

.

Evolution

58

:

1195

.

Muller

,

H. J.

1942

.

Isolating mechanisms, evolution, and temperature

.

Biol. Symp.

6

:

71

–

125

.

Mustonen

,

V.

, and

M.

Lassig

.

2007

.

Adaptations to fluctuating selection in Drosophila

.

Proc. Natl. Acad. Sci. USA

104

:

2277

.

Nei

,

M.

,

T.

Maruyama

, and

C. I.

Wu

.

1983

.

Models of evolution of reproductive isolation

.

Genetics

103

:

557

.

Nuzhdin

,

S. V.

,

M. L.

Wayne

,

K. L.

Harmon

, and

L. M.

McIntyre

.

2004

.

Common pattern of evolution of gene expression level and protein sequence in Drosophila

.

Mol. Biol. Evol.

21

:

1308

–

1317

.

Orr

,

H. A.

1995

.

The population genetics of speciation: the evolution of hybrid incompatibilities

.

Genetics

139

:

1805

–

1813

.

Orr

,

H. A.

1996

.

Dobzhansky, Bateson, and the genetics of speciation

.

Genetics

144

:

1331

.

Orr

,

H. A.

2005

.

The probability of parallel evolution

.

Evolution

59

:

216

.

Orr

,

H. A.

, and

M.

Turelli

.

2001

.

The evolution of postzygotic isolation: accumulating Dobzhansky-Muller incompatibilities

.

Evolution

55

:

1085

.

Philippe

,

N.

,

E.

Crozat

,

R. E.

Lenski

, and

D.

Schneider

.

2007

.

Evolution of global regulatory networks during a long-term experiment with Escherichia coli

.

Bioessays

29

:

846

–

860

.

Porter

,

A. H.

, and

N. A.

Johnson

.

2002

.

Speciation despite gene flow when developmental pathways evolve

.

Evolution

56

:

2103

.

Presgraves

,

D. C.

2002

.

Patterns of postzygotic isolation in Lepidoptera

.

Evolution

56

:

1168

.

Sasa

,

M. M.

,

P. T.

Chippindale

, and

N. A.

Johnson

.

1998

.

Patterns of postzygotic isolation in frogs

.

Evolution

52

:

1811

.

Siegal

,

M. L.

, and

A.

Bergman

.

2002

.

Waddington's canalization revisited: Developmental stability and evolution

.

Proc. Natl. Acad. Sci. USA

99

:

10528

–

10532

.

True

,

J. R.

, and

E. S.

Haag

.

2001

.

Developmental system drift and flexibility in evolutionary trajectories

.

Evol. Dev.

3

:

109

.

Van Nimwegen

,

E.

,

J. P.

Crutchfield

, and

M.

Huynen

.

1999

.

Neutral evolution of mutational robustness

.

Proc. Natl. Acad. Sci. USA

96

:

9716

–

9720

.

Wagner

,

A.

1996

.

Does evolutionary plasticity evolve?

Evolution

50

:

1008

–

1023

.

© 2009, Society for the Study of Evolution