Epigenetic induction may speed up or slow down speciation with gene flow.

Speciation is less likely to occur when there is gene flow between nascent species. Natural selection can oppose gene flow and promote speciation if there is variation in ecological conditions among the nascent species' locations. Previous theory on ecological speciation with gene flow has focused primarily on the role of genetic variation in ecological traits, largely neglecting the role of nongenetic inheritance or transgenerational plasticity. Here, we build and analyze models incorporating both genetic and epigenetic inheritance, the latter representing a form of nongenetic inheritance. We investigate the rate of speciation for a population that inhabits two patches connected by migration, and find that adaptively biased epigenetic induction can speed up or slow down speciation, depending on the form of the map from genotype and epigenotype to phenotype. While adaptively relevant epigenetic variation can speed up speciation by reducing the fitness of migrants and hybrids, it can also slow down speciation. This latter effect occurs when the epialleles are able to achieve adaptation faster than the genetic alleles, thereby weakening selection on the latter.

For spatially sub‐divided populations, phenotypic divergence occurs less readily when there is greater gene flow. Gene flow acts to homogenize populations, while disruptive natural selection drives populations phenotypically apart. Provided selection is sufficiently strong relative to migration, these opposing forces can lead to a migration‐selection balance, an equilibrium state where alleles favoured in each habitat are both maintained (Charlesworth & Charlesworth 2010, p.145). As the evolution of reproductive isolation can accompany ecological divergence (i.e., ecological speciation, Nosil 2012), both the probability and rate of speciation may be reduced by gene flow. Despite growing appreciation of the role of transgenerational plasticity and epigenetic inheritance in adaptive evolution (Klironomos et al. 2013, Kronholm & Collins 2015, Stajic & Jansen 2021), little attention has been given to its potential role in population divergence and speciation (but see ideas explored verbally in Pfennig & Servedio 2013). In this study, we develop a simulation model to explore how epigenetic variation may change the rate at which phenotypic divergence and speciation occur.

The role of nongenetic inheritance in adaptation, including that which is mediated by epigenetic inheritance of DNA methylation states, chromatin modifications, or small RNAs, has been largely overlooked until recently (Bonduriansky & Day 2018, Smith & Ritchie 2013). Accumulating evidence suggests that phenotypes can be inherited through nongenetic means and that epigenetic markers can be adaptively relevant and heritable (Herman & Sultan 2016, Kelley et al. 2021a, Shirai et al. 2021). Pointing to a possible role in population divergence, epigenetic variation has been found to correlate with environmental variables, a sign of local adaptation (Herrera & Bazaga 2010, Heckwolf et al. 2020, Kelley et al. 2021b, McNew et al. 2017, Platt et al. 2015, Richards et al. 2012, ,2017, Robertson et al. 2017, Thorson et al. 2017, Wilschut et al. 2016). Whether these patterns are due to within‐generation plastic responses associated with epigenetic changes or transgenerational inheritance of epigenetic factors, however, is often unclear, but compelling examples exist. In Polygonum plants, drought conditions led parents to produce offspring with longer roots, an adaptive phenotype that has been shown to be mediated by methylation (Herman & Sultan 2016); in stickleback fish, populations from marine environments with different salinities exhibit differences in methylation markers and laboratory experiments that vary salinity have shown that a subset of these same methylation differences are inducible and exhibit transgenerational plasticity, suggesting a role of adaptively inducible epigenetic variation in nature (Heckwolf et al. 2020); in Daphnia, stress‐induced epigenetic modifications occurred and persisted up to four generations, with modification occurring especially in stress‐related genes, suggesting a possible role of heritable epigenetic change in adaptation (Feiner et al. 2021).

In addition to the potential role of epigenetic variation in enhancing local adaptation, transient interactions between genetic and epigenetic variation may be influential during the speciation process in determining the tempo with which genetically based ecological divergence or speciation occurs, even if the signature of these interactions is not detectable in studies of most contemporary species. By building and analyzing simulation models, we can provide first steps in understanding the potential role epigenetic variation plays in the speciation process and, further, provide hypotheses that might be tested empirically by studies conducted in incipient or recently diverged species (Stankowski & Ravinet 2021) where signatures of transient epigenetic influences may still be visible.

It is difficult a priori to predict the impact of adaptively relevant epigenetic variation on population divergence or speciation, as intuition might lead us to hypothesize that epigenetic variation could either increase or decrease the probability or rate of speciation. On the one hand, because epigenetic markers likely have higher (epi‐)mutation (i.e., induction) rates (van der Graaf et al. 2015), epigenetic variation could kick‐start the evolution of reproductive isolation in a manner similar to phenotypic plasticity (Pfennig et al. 2010). For example, in darter fish that are thought to be in the early stages of speciation, divergence for methylation markers was greater than divergence for genetic markers (Smith et al. 2016). Intriguingly, in this system, there also appears to be a correlation between epigenetic differentiation and reproductive isolation mediated by behavioral differences.

On the other hand, epigenetic variation, by allowing a population to adapt quickly, could interfere with adaptation through genetic evolution (Kronholm & Collins 2015), a process of interference akin to Hill–Robertson interference (Hill & Robertson 1966), which has been found to occur in yeast (Stajic et al. 2019). If epigenetic changes are unstable, ecologically relevant phenotypic divergence wrought primarily through them may not represent a stable long‐term evolutionary phenomenon and should, therefore, not be termed “speciation.” Using theoretical modeling, we seek to clarify expectations about whether epigenetic variation and adaptively biased epimutation promote or interfere with genetically based population divergence and speciation.

Materials and Methods

MODEL OVERVIEW

We consider an individual‐based stochastic model of a haploid, sexual population divided between two patches, each of fixed population size, K, that are connected by migration (see Table 1 for a summary of parameter definitions). The patches differ in the direction of selection that acts on an ecologically relevant trait, thereby providing an opportunity for ecological divergence and local adaptation. This trait is influenced by genetic loci and epigenetic loci (Fig. S1), which we also refer to as epigenetic markers. All loci, genetic and epigenetic, are freely recombining with Mendelian segregation, and have two alternative alleles or “epialleles,” respectively.

Table 1

Model parameters, definitions, and values used. For parameters with a range of values considered, those values in boldface are the ones used in Figs. 1, 2, 3, 4, 5, S3, and S5–S10

Parameter	Definition	Values
m	Migration probability between the patches.	0,136,118,112,19,16,13,12
n_ecol;g	Total number of genetic loci influencing the ecological trait.	10
n_ecol;e	Total number of epigenetic markers influencing the ecological trait.	10
n T	Number of genetic loci for the ornament.	20
n P	Number of genetic loci for the preference.	20
K	Fixed population size in each patch.	2500
α	Strength of sexual selection.	5
β	Cost of male ornament.	1
γ	Cost of female preference.	0.02
μ_epi;A	Rate of adaptive induction at adaptively‐biased epigenetic loci.	10^−5, 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.3, 0.4, 0.5
μ_epi:M	Rate of maladaptive induction at adaptively‐biased epigenetic loci.	10^−5
μgen	Mutation rate at genetic loci.	10^−5
σ	Width of natural selection function.	0.35 , 0.4
τ	Persistence time for epigenetic induction.	1,2,3

Figure 1

Sample dynamics for 35 replicate runs under our “core” model (additive genotype/epigenotype to phenotype map and hard selection). Parameter values are those highlighted in Table 1 with , and . (a) Mean frequency of genetic 1 alleles for ecological adaptation in patch 0 (black) and patch 1 (grey). (b) Final score of female preference with the average across replicates shown in grey. (c) The frequency of hybrids. (d) The fitness of migrants under natural selection.

Figure 2

Time to speciation is generally shorter for higher rates of adaptively biased epimutation () across a range of migration rates (m) and epigenetic stabilities (τ) under our “core” model. Points show mean speciation time when epigentic induction is adaptively biased (black) or unbiased (blue; adaptively biased epimutation rate and maladaptive epimutation rate are both equal to 10−5 for left‐most points); error bars denote standard error. Each point represents the average over 35 replicate model runs. The dotted line shows the time to speciation for a model with 20 genetic loci and no epigenetic loci, and its standard error is displayed as a bar on the left side. Parameters values are listed in Table 1 with reset times for epimutations, τ, equal to 1 (light grey), 2 (dark grey), and 3 (black). Note the differing scales for the y ‐axes among panels.

Figure 3

Migrant fitness due to natural selection (averaged through time at evenly spaced time points) varies with the rate of adaptively biased epimutation () across a range of migration rates (m) and epigenetic stabilities (τ) under our core model. Here, migrant fitness is simply calculated from simulations as the ratio , where is the frequency of migrants before the natural selection step and is the fraction of migrants directly after natural selection. Points show mean fitness when epigentic induction is adaptively biased (black) or unbiased (blue; adaptively biased epimutation rate and maladaptive epimutation rate are both equal to 10−5 for left‐most points); error bars denote standard error (and are too small to see in this figure). Each point represents the average over 35 replicate model runs. The dotted line shows fitness for a model with 20 genetic loci and no epigenetic loci, and its standard error is displayed as a bar on the left side. Parameters values are listed in Table 1 with line shading again denoting reset time for epimutations; (light grey), (dark grey) and (black).

Figure 4

Time to speciation varies with the rate of adaptively biased epimutation () across a range of migration rates (m) and epigenetic stabilities (τ) under a soft selection scenario. Points show mean speciation time when epigentic induction is adaptively biased (black) or unbiased (blue; adaptively biased epimutation rate and maladaptive epimutation rate are both equal to 10−5 for left‐most points); error bars denote standard error. Each point represents the average over 35 replicate model runs. The dotted line shows the time to speciation for a model with 20 genetic loci and no epigenetic loci, and its standard error is displayed as a bar on the left side. As in our core model, we assume an additive genotype/epigenotype to phenotype map. Parameters values are listed in Table 1 with (light grey), (dark grey), and (black). Note the differing scales for the y ‐axes among panels.

Figure 5

Time to speciation varies with the rate of adaptively‐biased epimutation () across a range of migration rates (m) and epigenetic stabilities (τ) under a “redundant” genotype/epigenotype to phenotype map. Points show mean speciation time when epigentic induction is adaptively biased (black) or unbiased (blue; adaptively biased epimutation rate and maladaptive epimutation rate are both equal to 10−5 for left‐most points); error bars denote standard error. Each point represents the average over 35 replicate model runs. The dotted line shows the time to speciation for a model with 20 genetic loci and no epigenetic loci, and its standard error is displayed as a bar on the left side. As in our “core,” we assume hard selection. Parameters values are listed in Table 1 with (light grey), (dark grey) and (black). Note the differing scales for the y ‐axes among panels.

A single generation comprises selection, mate choice, sexual reproduction, mutation, epigenetic induction (also referred to as epimutation) and migration, in that order (Fig. S2) (although we also briefly investigate the effect of having migration precede mutation and epigenetic induction). We describe these processes in more detail below.

To be concrete, consider the following brief example of a scenario we are modeling. We begin with one species that is genetically and epigenetically uniform, and adapted to one particular habitat. Then, following a disruptive vicariance event, the population becomes separated into two subpopulations (connected with some amount of migration) with different ecological conditions, one similar to the original habitat and the other not. Genetic change at multiple loci or epigenetic induction at various epigenetic markers are both routes towards adaptation to the novel habitat. We are interested in the speed of genetically based ecological divergence between the subpopulations and concomitant speciation owing to the evolution of reproductive isolation.

EPIGENOTYPE AND GENOTYPE TO PHENOTYPE MAPPING

An individual's ecologically relevant phenotype is determined by its combination of genetic alleles and epialleles (Fig. S1). We use the term epiallele to refer to the state of an epigenetic marker, and these states can be reached by unbiased epimutation or adaptive environmental induction (i.e., with an elevated rate of production of locally adaptive epialleles). When epialleles reset after a single generation, this process of adaptively biased epigenetic induction models within‐generation plasticity, whereas when these epialleles take multiple generations to reset, this process models a system of transgenerational plasticity. We consider both scenarios here. Empirical evidence for adaptively relevant epigenetic variation comes from studies of plants (Wilschut et al. 2016), unicellular organisms (Kronholm et al. 2017), and animals (Hu & Barrett 2017). Recent evidence also suggests that epimutations can experience positive natural selection in the form of selective sweeps (Shirai et al. 2021).

We initially assume that the map from genetic and epigenetic states to phenotype is additive (termed the “additive” map), but we also consider an alternative map, as noted below. According to our “additive” map, the ecological phenotype, x (with ), of an individual is where represents the contribution of genetic allele i and represents the contribution of epiallele i. Each genetic and epigenetic allele assumes values of either 0 or 1. We assume that an epigenetic induced state gives a value of 1 and produces an identical phenotypic effect to a genetic 1 allele, and that these epialleles or genetic alleles are favored in the novel habitat (which we denote patch 1). Likewise, uninduced epigenetic markers have an equivalent effect on phenotype as genetic alleles of value 0, which are favoured in the ancestral habitat (which we denote patch 0). We denote the total number of genetic loci by , the total number of epigenetic loci by , and the total number of genetic or epigenetic loci by .

In the Supporting Information (S3), we present an alternative map, in which we introduce interference between the genetic and epigenetic systems by having genetic alleles redundant with epigenetic markers. The purpose of this alternative map, which we call the “redundant” map, is to investigate a situation in which genetic evolution or epigenetic evolution compete to “accomplish” adaptation. For comparison, we also present results from a purely genetic model with n ecol genetic loci and no epigenetic markers.

Ornaments and preferences, the basis for mate choice in our model, are determined by and genetic loci, respectively (Fig. S1). We assume that ornament quality in a male may be further influenced by that male's physical condition, as determined by the match between his ecological trait and the environment. This condition‐dependent form of sexual selection is premised on the assumption that mutations that affect an individual's fit to its environment also affect the quality of his ornamentation and, thereby, attractiveness. In theoretical studies, female preferences for condition‐dependent male ornaments evolve and condition‐dependent sexual selection can promote speciation with gene flow without requiring so‐called “magic” traits (i.e., traits that both contribute to local adaptation and serve as cues for assortative mating) (Proulx & Servedio 2009, Proulx 2001, Reinhold 2004, Schindler et al. 2013, van Doorn et al. 2009, Veen & Otto 2015). By building a model in which speciation with gene flow is likely, we can conveniently study how epigenetic induction affects the speed with which speciation occurs.

A male's ornament, s, is determined by the mean of his alleles (each taking a value of 0 or 1) at n T genetic loci, referred to as T ‐loci, and, because ornaments are condition‐dependent, on his ecological fitness. Specifically, the ornamentation of a male in patch k with ecological trait x who possesses j 1 alleles at his T loci is given by , where denotes his ecologically based fitness (calculation of this quantity is described below). Essentially, this multiplication means that the male gains units of ornamentation for every unit of ecological fitness.

A female's preference strength, a, is determined by her alleles at n P genetic P ‐loci (each taking a value of 0 or 1) and, for a female with i 1 alleles, is calculated as . The T ‐loci have no effect in females and the P ‐loci have no effect in males.

NATURAL SELECTION

We consider two patches, denoted patch 0 and patch 1, that differ in the direction of selection on an ecological trait. As described above, an individual's ecological trait is a function of its alleles and epialleles at genetic and epigenetic loci. We assume that allelic/epiallelic values of 0 and 1 are optimal in patch 0 and patch 1, respectively. We then calculate the ecological fitness of individual i with phenotype x in patch k as where is the optimum phenotype in patch k (0 in patch 0; 1 in patch 1) and σ measures the strength of natural selection within each patch (larger σ corresponds to weaker selection).

In addition, mate choice traits (preference for females and ornament for males) can impose fitness costs. Specifically, the total fitness of a female with i preference alleles of value 1 is equal to , where and the total fitness of a male with j ornament alleles of value 1 is equal to , where . Here, γ and β denote the rate of decrease in fitness as an individual allocates more toward mate choice traits (e.g., more 1 alleles at the relevant mate choice loci). Note that this cost only affects a male's survival probability and not his ornament quality (and, therefore, not his likelihood of being chosen as a mate) if he survives.

Natural selection occurs via survival; an individual's survival probability is given by its total fitness. We refer to this selection regimen as “hard selection” (Wallace 1975). In the Supplementary Material (S4), we also describe an alternative “soft selection” regime where, during selection, total population size is held constant. We use this scenario to investigate the importance of population bottlenecks that emerge under the hard selection scenario.

MATE CHOICE AND REPRODUCTION

Surviving males express ornaments and surviving females choose among these males according to their preferences (although our model applies equally to species with male preferences and female ornaments). Mate choice occurs according to a fixed relative‐preference scheme (Kirkpatrick 1982). Under such an assumption, a particular male's likelihood of being chosen by a given female depends on his attractiveness to her relative to all other males.

We assume that every female mates, so there are no emergent fecundity costs associated with being choosy. In patch k, the probability that a female with preference mates with a male with ornament is then where α is the strength of sexual selection, is the frequency of males with ornament in patch k, and is the total attractiveness of males in the population, ensuring that probabilities sum to 1. We assume equal reproductive fitness (i.e., fertility) among females so that, within each patch, the number of offspring produced by each mother is drawn from a multinomial distribution, , where K is the fixed population size of each patch. We assume free recombination among all loci and epigenetic markers and random segregation in a transient diploid stage.

MUTATION AND EPIMUTATION

At all genetic loci (ecological trait, ornament, and preference), mutations (0 mutates to 1 and 1 to 0) occur without bias and with probability μgen, per locus. At epigenetic loci, epigenetic induction (0 changes to 1) and resetting (1 reverts back to 0) occurs. These changes could correspond to transitions between methylated and unmethylated states. Note that because the 1 epiallele is favored in patch 1 and the 0 epiallele is favored in patch 0, epigenetic induction is adaptive in patch 1 and maladaptive in patch 0, while epigenetic resetting is adaptive in patch 0 and maladaptive in patch 1. In patch 1, adaptively biased epimutation occurs stochastically at rate , while in patch 0, maladaptive epimutation    occurs at a rate, . We consider cases with ; when the epigenetic induction process is adaptively biased and when it is unbiased. Resetting occurs deterministically after τ generations, allowing us to consider within‐generation plasticity () alongside transgenerational plasticity (). Assumptions of short duration and high mutation rates for epigenetic markers is supported by recent research on plants and animals (Beltran et al. 2020, van der Graaf et al. 2015), albeit for processes of adaptively unbiased epimutation. As we allow for adaptive bias in epigenetic induction, our approach to modeling epimutation differs from that of existing models of adaptation with unbiased epimutation (Klironomos et al. 2013, Kronholm & Collins 2015). While definitive empirical evidence of adaptively biased epimutation is currently lacking, our approach to modeling adaptively biased epimutation is supported by empirical data that show transgenerational plasticity (Agrawal et al. 1999, Beemelmanns & Roth 2017, Galloway & Etterson 2007, Herman & Sultan 2016, Hales et al. 2017, Heckwolf et al. 2020, Sobral et al. 2021), which may be mediated by epigenetic markers, such as methylation.

MIGRATION

Each individual migrates to the other patch with probability m. When , no migration occurs, and when , migration fully mixes the two patches every generation.

ANALYSIS

We initialize model runs with genetic and epigenetic loci all fixed at 0, meaning individuals are perfectly adapted to the “ancestral” habitat (patch 0). We run simulations for 4 × 104 generations and record when speciation occurs. We use a speciation criterion (described fully in the Supplementary Material) that is based on the degree of bimodality in the distribution of genetic values of the ecological trait. Our metric yields a value of 1 for an initially unimodal population and then, as genetic differentiation increases, the metric decreases and we can conclude that speciation has occurred once it falls below a threshold. We base our calculation of this metric entirely on the genetic loci so that “speciation” means that significant divergence has occurred at genetic loci. We also confirm that female preferences for condition‐dependent traits evolve to become stronger and that correlations between ecological traits and female preferences develop among replicate simulations, both indicators of speciation. The model we have implemented here has previously been shown to lead to ecological speciation (van Doorn et al. 2009), measured via phenotypic divergence and the evolution of female preferences for locally adapted traits. Our interest here is not in asking whether or not speciation occurs but, instead, in identifying how the rate of speciation depends on adaptively biased or unbiased epigenetic induction.

We calculate average speciation time across 35 replicate model runs. We judge speciation to be faster for one parameter set compared to another when the mean speciation time for those replicates is lower. For runs in which speciation does not occur, we cannot conclude that speciation will not occur if we run the model longer. For these cases, we record their “speciation” times as the total number of generations, in order to be able to make qualitative comparisons of speciation time among parameter sets. Simulations were written in R (R Development Core Team 2008).

Results

OVERVIEW

We explore how adaptively biased or unbiased epigenetic induction affects the evolution of ecological divergence and speciation. It is important to note that what we call speciation here is based on a bifurcation at the genetic loci only (and not epigenetic loci). Because they revert, epigenetic changes, and their contribution to phenotypic divergence, do not yield stable phenotypes capable of long‐term persistence. For this reason, we reserve the term “speciation” to describe genetic changes. Moreover, as a modeling convenience, defining speciation in this way makes comparisons among scenarios with different rates of epigenetic induction straightforward, because we are only comparing dynamics at the genetic components where mutation rates are common across all cases.

We confirmed previous findings that the combination of divergent natural selection and condition‐dependent sexual selection promotes speciation (van Doorn et al. 2009). We base this determination on the genetic bifurcation of the ecological trait (Fig. 1a), as well as the evolution of stronger female preferences for locally‐adapted (i.e., condition‐dependent) male traits (Fig. 1b), reduction in the frequency of hybrids (Fig. 1c), reduction in the fitness of migrants (Fig. 1d), and correlation among simulation replicates between female preference and ecological adaptation over time (Fig. S3), indicating a process of ecological speciation in which ecological divergence and sexual selection‐mediated isolation are linked. The correlations become largest when migration is intermediate. When correlations do not build up because divergence of the ecological trait occurs rapidly without the facilitation of sexual selection, and when m is high (e.g., ) correlations do not build up readily because of the influx of maladapted ecological alleles linked to female preference alleles.

SPECIATION RATE UNDER CORE MODEL

We observe significant variance in the speed with which speciation occurs. Under our core model (Fig. 2, additive phenotype map and hard selection), unbiased epimutation slows down speciation compared to a purely genetic model (i.e., the blue points in Fig. 2, which correspond to cases where the adaptive epimutation rate and maladaptive epimutation rate are equal, are all above the dotted line, which denotes the speciation time for a scenario without epimutation). With adaptively biased epimutation, increasing the rate of adaptively biased epimutation speeds up speciation (Fig. 2). For sufficiently high rates of adaptively biased epimutation, speciation reliably occurs faster than it would under a purely genetic model (curves in Fig. 2 drop below the dotted line). Higher stability of the epigenetic markers also speeds up speciation (darker lines in Fig. 2 are below paler lines), except when there is no migration and epimutation is unbiased (i.e., the rate of adaptively biased epimutation in patch 1, , equals the rate of maladaptive epimutation in patch 0, ). In line with expectations based on earlier work, a higher rate of migration between patches slows down speciation (Fig. 2), as does weaker selection (compare Fig. 2 and Fig. S4 for ; the speed up with weaker selection when appears to be due to the effect of bottlenecking under “hard” selection, explained in more detail below).

To understand how increasing the rate of adaptive epigenetic induction speeds up speciation, we measure its effects on the fitness of migrants and hybrids. We track survival and reproductive contributions of migrant individuals and find that the fitness of migrants declines with a higher rate of adaptive epigenetic induction both under natural selection (Fig. 3) and sexual selection (Fig. S5). Increasing the rate of adaptive epigenetic induction renders an individual more adapted to its current environment, which lowers its survival and its chance at procuring a mate in the alternative habitat. This ultimately lowers the effective migration rate between the two patches, facilitating speciation. A lower effective migration with greater epigenetic induction is also manifest in the decline in frequency of hybrids (individuals with one parent from patch 0 and one parent from patch 1; Fig. S6).

SPECIATION RATE UNDER SOFT SELECTION REGIME

In models with epigenetic variation, speciation is faster with greater adaptive epigenetic induction bias even when there is no migration. This indicates that there is an effect of adaptive epigenetic induction beyond the migration effects described above. In our core model, selection occurs through survival (i.e., what we term “hard” selection), meaning that every generation a potential bottleneck occurs before the production of offspring. With greater adaptive epigenetic bias, fitnesses tend to be higher and, consequently, the size of such bottlenecks smaller (Fig. S7). Less severe bottlenecks lower the effects of genetic drift and, thus, with greater adaptive epigenetic bias, ecological divergence between the patches at the genetic loci is faster. To validate the above logic, we ran our model under a “soft” selection scenario, which elimates bottlenecking completely (described in Supporting Information S4). Indeed, we find that speciation time no longer declines with increasing rate of adaptive epigenetic induction when there is no migration (Fig. 4a). In fact speciation times become longer with greater epigenetic induction, suggesting that epigenetic and genetic systems are weakly interfering with one other in achieving local adaptation. However, with non‐zero migration rates, increasing adaptive bias in epimutation consistently leads to faster speciation (Fig. 4b–d). We conjecture that this same process of recurrent bottlenecking leads to slightly “faster” speciation when hard selection is weak and there is no migration (compare Fig. 2a to Fig. S4A, the former showing results for stronger selection): with weaker selection on epigenetic markers, the bottleneck experienced by genetic alleles is smaller, allowing them to diverge faster, an effect of drift that is obscured when there is migration (i.e., comparing Fig. 2b–d and Fig. S4B–D, weaker selection slows down speciation when ). When migration rates are very small but non‐zero, weaker selection still slightly speeds speed up speciation (, compare Figs. S11A and S12A), but this effect disappears as migration rates increase (e.g., , compare Figs. S11B and S12B).

SPECIATION RATE UNDER REDUNDANT PHENOTYPIC MAPPING

Although increasing the rate of adaptive epigenetic induction can increase the rate of speciation, as described above, increasing the rate of such induction may also slow down speciation. By providing an alternative route to adaptation, adaptive epigenetic induction could allow populations to rapidly gain fitness in novel environments without the need for genetic evolution, whose pace is limited by modest genetic mutation rates. While this effect was not obviously visible in the previously examined core model, it becomes evident when we instead consider a “redundant” map from genotype and epigenotype to phenotype which strengthens competition between genetic and epigenetic loci. In this scenario (described in Supporting Information S3), each genetic locus is paired with and redundant to one epigenetic marker. This assumption means that perfect adaptation can be achieved entirely by changes at epialleles. The rapid increase in fitness achieved by induced epialleles subsequently relaxes selection on genetic alleles, thereby slowing the rate of fixation of beneficial genetic alleles. Because we measure speciation only at genetic loci, epigenetic interference with genetic evolution can slow down speciation, as defined here, both under hard (Fig. 5) and soft (Fig. S8) selection. Similarly, when migration occurs before epigenetic induction (instead of after, as above), adaptively biased epigenetic induction can slow speciation even under an “additive” map (Fig. S9). This effect was anticipated by Thibert‐Plante & Hendry (2011), who showed that less reproductive isolation evolved when plasticity was expressed after migration. In both cases, when plasticity or adaptive epigenetic induction occurs after migration, nongenetic changes help to align an individual's phenotype with its selective environment, thereby weakening selection at genetic loci, which is what drives speciation.

ALTERNATIVE SPECIATION METRIC

To buttress our assessment of speciation that is based on ecological trait bifurcation, we also considered an alternative speciation metric that quantifies the degree to which mating occurs within groups of genetically similar individuals, with higher values indicating greater levels of assortative mating, with respect to the ecological trait. We explain this metric in more detail in Supporting Information S2. We find that our conclusions are largely supported when using this alternative metric: higher assortative mating evolves early on (Fig. S10) for parameter sets exhibiting faster speciation (Fig. 2).

LOWER EPIMUTATION RATES

By considering a range of adaptive epigenetic induction rates, we explore the array of possible effects of epimutation on speciation. However, the strength and prevalence of adaptive epigenetic induction in nature remains difficult to quantify. Thus, we also conduct a set of model runs where epigenetic inducation is relatively weaker than for those results presented so far. We find that many of the trends that are visible between adaptively biased epimutation rates of 0 and 0.1, in the aforementioned figures (Figs. 2, 4, 5, S4, S8, S9) remain for even lower rates (Figs. S13– S18, respectively) , indicating that even if adaptive epimutation rates are slight, they may still have similar impacts on adaptation and speciation rates.

Discussion

We have shown that unbiased epigenetic induction either slows or has little effect on the rate of ecological speciation compared to a purely genetic model. Inclusion of adaptive bias to epigenetic induction usually increases the rate of speciation relative to the unbiased case. Further increases of the rate of adaptive epimutation often speeds up speciation as well, but may slow it down in some cases.

Because biodiversity tends to increase with the rate of speciation, our results suggest that nongenetic inheritance may play a role in generating and maintaining biodiversity. Adaptively induced epigenetic variation is a mechanism for producing phenotypic plasticity, either within‐generation plasticity, when markers last only a single generation, or transgenerational plasticity, when markers can persist across multiple generations. Although we have focused here on epigenetic inheritance, our work should also apply to other nongenetic mechanisms that enable parents to transmit phenotypes to offspring (e.g., cytoplasmic inheritance). We emphasize that our model is designed to qualitatively guide intuition about how epigenetic variation may affect speciation and not to make quantitative predictions. Hence our model serves as a proof‐of‐concept test of verbal hypotheses of biological processes (Servedio et al. 2014).

Within‐generation nonheritable phenotypic plasticity has previously been suggested to promote speciation (Pfennig et al. 2010). Theory has shown that, whether plasticity promotes or inhibits speciation can depend on whether plasticity is expressed before or after migration (Thibert‐Plante & Hendry 2011). When plasticity occurs after migration, it helps move individual phenotypes towards their local adaptive peak, which can reduce the strength of selection on genetic alleles and thereby slow speciation. Conversely, plasticity before migration can strengthen selection on genetic alleles (Thibert‐Plante & Hendry 2011). Here, we also find that the order of these events matters: when adaptive epigenetic induction occurs prior to migration it promotes speciation to a greater extent than when induction occurs after migration. Our model does not include adaptive habitat choice, which can also promote speciation (Nonaka et al. 2015). Exploring adaptive habitat choice in the context of speciation would be an interesting avenue for future work.

Theoretical work has generated hypotheses about how epigenetic variation may evolve (Geoghegan & Spencer 2012, 2013), when epigenetic systems for transgenerational plasticity are most likely to evolve (Furrow & Feldman 2013, Greenspoon & Spencer 2018), and how the contribution of epimutation to adaptation in a peripheral population depends on the rate of epimutation and its adaptive bias (Smithson et al. 2019). However, in nature the broad‐scale importance of epigenetic variation in evolution remains uncertain. There are a number of situations where epigenetic variation does appear to play an important role, including heritable phenotypic variation (Blewitt et al. 2006, Cubas et al. 1999, Duempelmann et al. 2019), stress‐induced heritable epigenetic variation (Verhoeven et al. 2010), environmentally induced epigenetic traits that can be heritable over multiple generations (Klosin et al. 2017, Rechavi et al. 2014), and fitness effects that are caused by heritable epigenetic changes (Herman & Sultan 2016, Shirai et al. 2021). Yet our understanding of these phenomena is still in its infancy, and current data on epigenetic variation and transgenerational plasticity do not permit estimation of the rate of adaptively biased epimutation, although estimates of unbiased epimutation have been made (van der Graaf et al. 2015). Owing to this uncertainty, we considered a range of rates and found that even slight adaptive bias in epimutation can have impacts on the rate of speciation.

Our work underscores the importance of understanding how genetic and epigenetic factors combine and map onto phenotype and fitness. We show that redundant maps, which create greater opportunity for interference between epigenetic and genetic loci are likely to hinder speciation compared to nonredundant maps. Because adaptively relevant epigenetic variation can arise faster than genetic variation, epigenetic loci can quickly improve an individual's fit to its environment and, consequently, reduce the strength of selection on genetic loci. Determining how genetic and epigenetic information combine in determining phenotype is an important open question for empirical work.

Comparison between our hard and soft selection scenarios suggests that adaptively relevant epigenetic variation may interact with genetic drift to affect genetic evolution. By rapidly enhancing adaptation, epigenetic variation reduced the opportunity for bottlenecking and the effect of drift, thereby promoting genetically based ecological adaptation. Future theoretical work should focus specifically on understanding this effect.

Theory generally predicts that genetic differentiation between incipient species will be higher in genomic regions with lower recombination, implying that these regions will play a greater role in the process of speciation (Nachman & Payseur 2012). To keep our model computationally‐manageable and comparable to related earlier work (e.g., van Doorn et al. (2009)), we have made the simplifying assumption that recombination occurs freely among all loci. Although we might speculate that speciation would be faster with lower rates of recombination among loci, we conjecture that the qualitative effects of epigenetic variation that we have identified here should not depend critically on the rate of recombination, as recombination should affect genomes with or without epigenetic variation in similar ways. Nevertheless, this hypothesis would be worth examining.

Previous work has shown that epigenetic inheritance can alter the tempo of adaptation on landscapes with a single fitness peak, either speeding it up or slowing it down (Klironomos et al. 2013, Kronholm & Collins 2015). Our work has shown that ecologically relevant epigenetic variation can also speed up or slow down the rate of evolutionary diversification. Empirical understanding of the extent and nature of heritable epigenetic variation across the tree of life may, thus, prove crucial for understanding rates of diversification and levels of biodiversity.

DATA ARCHIVING

Simulation code and data can be found on Dryad (https://doi.org/10.5061/dryad.573n5tb9d) (Greenspoon et al. 2022).

CONFLICT OF INTEREST

The authors have declared no conflict of interest.

Associate Editor: R. Unckless

Handling Editor: T. Chapman

Figure S1: For each individual of the population, genotype and epigenotype is mapped onto phenotype.

Figure S2: The life cycle of the population. For information about each stage, see the Materials and Methods section.

Figure S3: The correlation between the strength of the female preference and the locally adapted genotype over time in patch 1.

Figure S4: Time to speciation is generally shorter for higher rates of adaptively‐biased epimutation (μ epi;A) when selection is weak.

Figure S5: Mean migrant fitness due to sexual selection on males (averaged through time at evenly spaced time points) varies with the rate of adaptively‐biased epimutation (μ epi;A) under our core model.

Figure S6: Mean hybrid frequency (averaged through time at evenly spaced time points) is generally lower for higher rates of adaptively‐biased epimutation (μ epi;A) under our core model.

Figure S7: Mean population size after selection (averaged through time at evenly spaced time points) varies with adaptively‐biased epimutation rate (μ epi;A), demonstrating the bottlenecking that recurrently occurs under hard selection.

Figure S8: Time to speciation with soft selection under a redundant genotype/epigenotype to phenotype map.

Figure S9: Time to speciation under our core model when genetic and epigenetic mutation happen after migration.

Figure S10: The assortative mating index (averaged through time over the first 1000 time steps at evenly spaced time points) defined in Section S2 under our core model.

Figure S11: Analog of Fig. 2, but with lower migration rates. Parameters values are listed in Table 1 with τ = 1 (light grey), τ = 2 (dark grey) and τ = 3 (black).

Figure S12: Analog of Fig. S4, but with lower migration rates.

Figure S13: Analog of Fig. 2, but with smaller values of the adaptively biased epimutation rate.

Figure S14: Analog of Fig. 4 (soft selection scenario), but with smaller values of the adaptively biased epimutation rate.

Figure S15: Analog of Fig. 5 (redundant map), but with smaller values of the adaptively biased epimutation rate.

Figure S16: Analog of Fig. S4 (weak selection case), but with smaller values of the adaptively biased epimutation rate.

Figure S17: Analog of Fig. S8 (redundant map and soft selection), but with smaller values of the adaptively biased epimutation rate.

Figure S18: Analog of Fig. S9 (mutation after migration), but with smaller values of the adaptively biased epimutation rate.

Click here for additional data file.

Supporting information

Click here for additional data file.