Proceedings of the Royal Society B: Biological Sciences
You have accessResearch articles

The phoenix hypothesis of speciation

Ryo Yamaguchi

Ryo Yamaguchi

Department of Advanced Transdisciplinary Science, Hokkaido University, Sapporo, Hokkaido 060-0808, Japan

Department of Zoology & Biodiversity Research Centre, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z4

[email protected]

Contribution: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

,
Bryn Wiley

Bryn Wiley

Department of Zoology & Biodiversity Research Centre, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z4

Contribution: Formal analysis, Investigation, Visualization, Writing – review & editing

Google Scholar

Find this author on PubMed

and
Sarah P. Otto

Sarah P. Otto

Department of Zoology & Biodiversity Research Centre, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z4

Contribution: Conceptualization, Funding acquisition, Methodology, Supervision, Visualization, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Published:https://doi.org/10.1098/rspb.2022.1186

Abstract

Genetic divergence among allopatric populations builds reproductive isolation over time. This process is accelerated when populations face a changing environment that allows large-effect mutational differences to accumulate, but abrupt change also places populations at risk of extinction. Here we use simulations of Fisher's geometric model with explicit population dynamics to explore the genetic changes that occur in the face of environmental changes. Because evolutionary rescue leads to the fixation of mutations whose phenotypic effects are larger on average compared with populations not at risk of extinction, these mutations are thus more likely to lead to reproductive isolation. We refer to the formation of new species from the ashes of populations in decline as the phoenix hypothesis of speciation. The phoenix hypothesis predicts more substantial hybrid fitness breakdown among populations surviving a higher extinction risk. The hypothesis was supported when many loci underlie adaptation. With only a small number of potential rescue mutations, however, mutations that fixed in different populations were more likely to be identical, with such parallel changes reducing isolation. Consequently, reproductive isolation builds fastest in populations subject to an intermediate extinction risk, given a limited number of mutations available for adaptation.

1. Introduction

Biological diversity emerges from the interplay between speciation and extinction. Historically, these two factors have been studied as independent processes, but recent speciation studies have emphasized the need to consider population persistence for new taxa [13]. In addition, extreme environmental changes may affect the development of reproductive isolation as well as the risk of extinction, a possibility that we study through simulations. Here we explore how extinction risk itself could impact the likelihood of speciation.

Empirical studies have suggested that adaptation to a new environment facilitates the origin of species [4]. Yet, an abrupt environmental change poses a risk of extinction and species loss. Adaptation to a rapidly changing environment enables populations to accumulate phenotypically large-effect mutations, leading to severe fitness reduction in hybrids referred to as hybrid breakdown [57]. Although large-effect mutations were seen as unlikely to contribute to adaptation because small-effect changes are much more likely to be beneficial [8], there is a growing list of studies that have identified large-effect loci contributing to adaptive divergence, particularly in response to environmental shifts, such as Mimulus monkeyflower [9], Peromyscus oldfield mice [10], maize [11], Ectodysplasin (Eda) allele in three spine sticklebacks [1214], Heliconius Müllerian mimicry patterns [15], cave adaptation (Oca2) in Mexican tetras [16], and beak shape (Alx1) and beak size (Hmga2) in Darwin's finches [17,18]. For populations at risk of extinction, if they fail to adapt, those populations that do survive may be especially likely to accumulate large-effect mutations [19]. Here we explore the effect of evolutionary rescue on the distributions of fixed mutations and the consequences for speciation.

Experimental evolution studies have also confirmed that large-effect mutations contribute to adaptation in a variety of species adapting to a new environment. In these studies, large-effect mutations typically fix early during adaptation, followed by a series of many smaller effect substitutions (e.g. [2022]). Using Chlamydomonas, Collins & de Meaux [21] revealed that populations adapting to an abrupt environmental change showed evidence of substitutions of large-effect contributing to the earlier phase of adaptation, while populations adapting to a slowly changing environment showed patterns consistent with adaptation via mutations of small effect. Following an environmental change, Saccharomyces cerevisiae accumulates large-effect mutations and develops strong intrinsic incompatibilities, even when adapting to similar ecological changes [23,24]. These empirical studies suggest that environmental changes are likely engines of speciation for allopatric populations.

For populations that must adapt or face extinction (‘evolutionary rescue’ scenario, see Bell & Gonzales [25]), environmental stress initially causes the population to decline, which limits the potential for small-effect mutations to spread because selection is less effective in shrinking populations [26], leading to enrichment of large-effect substitutions among populations that can adapt and persist in the new environment [19]. This enrichment of large-effect mutations may thus drive higher rates of speciation. The influence of population dynamics and extinction risk on the evolutionary rate of reproductive isolation has not, however, been explored.

We investigate the hypothesis that extinction risk may elevate the chance of speciation, which we refer to as the phoenix hypothesis of speciation. The hypothesis is named after the Greek myth of the phoenix, which burns to ashes from which the next generation emerges. The phoenix hypothesis predicts that there is a greater degree of reproductive isolation among populations experiencing extinction risk, conditional on the populations successfully adapting and persisting.

Mutations producing large phenotypic changes may, however, be more likely to fix in parallel among different populations, because there is a limited set of possible genetic pathways to adaptation (as seen empirically, e.g. [27], and predicted theoretically, e.g. [28]; see also [29] and [30]). If the number of potentially beneficial mutations in a given new environment is small, the probability of evolving similar phenotypes based on shared genetic mechanisms (i.e. parallel evolution) is expected to be higher when the environment changes abruptly. Parallel genetic adaptation interferes with the development of reproductive isolation [31]. Thus, adaptation to extreme environmental shifts is expected to involve large-effect mutations, which contribute to reproductive isolation, but also parallel mutations, which do not. It thus remains unclear what the level of incompatibility is between populations adapting to the same environmental change that do survive extinction.

Thus, depending on the extinction risk of the parent population, we have contrasting expectations for reproductive isolation, with large but more parallel mutations accumulating during evolutionary rescue. To test the phoenix hypothesis, we extend Fisher's geometric model to simulate mutation-order speciation among diploid populations facing an extinction risk. When many genes could potentially contribute to rescue, we confirm the phoenix hypothesis and show that more reproductive isolation evolves in populations facing higher extinction risks (figure 1a). When the number of loci available for adaptation is finite and parallel mutations are more likely, a moderate extinction risk results in more extensive reproductive isolation among surviving populations (figure 1a). When facing a moderate extinction risk, populations are more likely to survive and adapt via different, large-effect mutations, generating substantial incompatibilities among populations.

Figure 1.

Figure 1. (a) The phoenix hypothesis of speciation. When extinction risk is high, large-effect mutations are likely to accumulate (green). Conversely, when extinction risk is low, small-effect mutations are likely to accumulate. Under the infinite-sites assumption (i.e. many loci in the panel), no common mutation accumulates among populations (darkest purple). A greater degree of reproductive isolation develops under a high extinction risk. When the number of loci available for adaptation is finite, large-effect mutations are more likely to be common among populations (lighter purple). Given the contrasting relationships with mutational effect sizes (green), an intermediate extinction risk can promote speciation. (b) An adaptive walk on the fitness landscape. Dots represent the position of the initial phenotype (white) and hypothetical phenotype at a certain point of adaptation (black) after accumulating four mutations (black arrows), shown both on the plane at the bottom (the phenotypic space) and then projected up onto the fitness landscape. (c) Genomic architectures of an individual are explored in this work. Two extreme examples of the number of loci are shown (L = 50 and 10 000). Please see the main text for the notation of parameters and variables. (d) Conceptual illustration of phenotypic parallelism. Here we show the case of a two-dimensional phenotype with a new optimum at the centre of the blue region, illustrating the adaptive paths taken by two independent populations (black and white). A white dot and grey dots represent the position of initial phenotype and hypothetical hybrids' phenotype, respectively. Red arrows are the mutations that are shared between the populations. (Online version in colour.)

2. Model

We conduct individual-based simulations inspired by Fisher's geometric model with explicit population dynamics of a sexually reproducing diploid organism with non-overlapping generations, using SLiM3.3 [32]. Table 1 summarizes the notation and gives the typical parameters used in this study.

Table 1. Notation and parameters used. Default values are given in square brackets, except where noted in the text. We run simulations for different values of L and u when their product (genome-wide mutation rate) is fixed as L × u = 5 × 10−3. Specifically, (L, u) = (50, 104), (100, 5 ×105) and (10 000, 5 ×107).

notation description
n number of phenotypes [4]
zi = {zi,1, zi,2, … ,zi,n} phenotype of ith individual
o = {o1, o2, … , on} optimal phenotype (o1 = 2 and otherwise 0)
q parameter translating the phenotypic distance from the optimum into the strength of selection [1]
k degree of epistasis [2]
v mutation size, drawn from an exponential distribution with mean λ= [0.4]
ξj, Δω random number drawn from a standard normal distribution
Δzi = {Δzi,1, Δzi,2, … , Δzi,n} phenotypic change of ith individual due to mutation
s, w strength of selection, fitness
N number of individuals
N0 initial number of individuals in a population [2000]
u mutation rate [10-4, 5 × 105, 5 × 107]
L genome size (number of sites) [50, 100, 10 000]
m recombination rate between consecutive sites [0.01]
r intrinsic growth rate (no selection) [0.8, 1.0, 1.4]
β constant birth and death rates [0.1]
K carrying capacity [2000]
c degree of pleiotropic side effect in mutation library [0, 0.5,1.0]

(a) Fitness landscape

Fisher's geometric model assumes multivariate stabilizing selection around a single optimum. Similar to Yamaguchi & Otto [7], we model this by defining the phenotype of an individual as an n-dimensional vector, where n is the number of traits relevant to adaptation in a given environment. The phenotype of the ith individual is denoted as zi = {zi,1, zi,2, … , zi,n}, where zi,j is the phenotypic value for trait j.

The fitness of an individual depends on the distance of its phenotypic values from the optimum o = {o1, o2, … , on}, measured as the Euclidean distance zio j=1n(zi,joj)2. We assume that the optimal phenotype is common to all individuals and populations. An individual's relative fitness is given by

w(zi)=exp(qziok),2.1
where q and k determine the decay rate and the curvature (the degree of epistasis) of the fitness landscape, respectively [33]. It should be noted that isotropic selection is assumed throughout the paper; meaning that there is equally strong selection on all n traits at a given distance to the optimum.

(b) Mutation and genetic architecture

To construct trait values from the accumulated mutations, we assume that each mutation additively contributes to the phenotypic value. Thus, an individual's phenotype is determined by the sum of effect sizes of all mutations affecting the traits plus the original phenotype (taken to be at the origin). While each mutation has an additive effect on the phenotype, the effect converted to fitness is not additive but depends on the location of the set of phenotypic values on the fitness landscape given by equation (2.1) (figure 1b).

For homozygous diploids, we define the effect of a mutation occurring in the ith individual by Δzi={Δzi,1,Δzi,2,,Δzi,n}, with Δzi,j describing the mutational effect on trait j. We assume that a single mutation affects each phenotype independently, with the change Δzi,j, proportional to the value drawn from a standard normal distribution, denoted as ξj [34,35]. The overall effect of a single mutation (its length across all phenotypic axes) was set to v, drawn from an exponential distribution with mean λ.

We also consider how the degree of pleiotropic side effects impacts the likelihood of parallel adaptation and the development of reproductive isolation. To vary the strength of pleiotropy, we allowed mutations to have a main effect on only one axis (weighted by an amount c), as well as a randomly drawn mutation vector (weighted by 1–c). Specifically, for each mutation, we randomly choose a focal trait k and set Δzi,k = ±cv. In addition, we let the mutation randomly affect each of the traits, j (including k), by an amount Δzi,j=(1c)vξj/l=1nξl2. This mutation scheme is exactly the same as the original Fisher's geometric model when c = 0 (isotropic mutation), while there is no pleiotropy if c = 1. The phenotypic value of the mutant individual was then set to zi + Δzi. The above describes the properties of homozygous mutations. Heterozygous diploids are assumed to be shifted half as far in phenotype space, Δzi/2 (i.e. we assume additivity on a phenotypic scale, although dominance for fitness depends on the mutation vector and position relative to the optimum).

Mutations can occur with equal probability u at each locus. The total number of biallelic loci is L. Recombination is modelled by assuming that cross-over events occur uniformly at random at rate m across loci (figure 1c). We prepare a mutation library to assign Δzi to each locus for each potential mutation that may occur in any of the populations during a given replicate simulation. By not assuming an infinite-sites model in which all novel mutations are unique, adaptive mutations can sometimes be shared between two populations evolving independently.

(c) Population dynamics

We consider a population of diploid individuals with a time-dependent population size N(t). We first calculate the extinction risk analytically, modelling the population dynamics as a birth–death process with birth and death rates b(t,N(t)) and d(t,N(t)), respectively:

N(t)N(t)+1:b(t,N(t))N(t)N(t)N(t)1:d(t,N(t))N(t).2.2

The expected change of N(t) over a small-time interval Δt is E[ΔN|N(t)] = r(t, N(t))N(tt, where r(t, N(t)) = b(t, N(t)) − d(t, N(t)). Following Uecker & Hermisson [36], we consider the process of mutation accumulation with a population that follows logistic growth (or decline) until it has reached its carrying capacity K. We assume that a decreasing availability of resources per individual leads to a lower birth rate and that the death rate depends on the strength of selection given by s(t, N(t)). Similar to eqns 20 in Uecker & Hermisson [36], the birth and death rates in the current model are given by

b(t,N(t))=β+r(1N(t)K)d(t,N(t))=β+s(t,N(t)).2.3

The parameter β is the constant birth and death rates of a population at carrying capacity if all individuals are at the optimum (no selection), and r is the intrinsic growth rate of the population (no selection). We define the average strength of selection as s(t,N(t))=1w¯(t), where w¯(t) represents the mean fitness of the population at the tth generation: w¯(t)=(1/N(t))i=1N(t)w(zi). If the entire population is at the optimum, then w¯(t)=1 holds, and the population size approaches the carrying capacity, K. Deviations from the optimum, however, induce selection, which reduces the steady-state population size. In our model, 0w¯(t)1 is always satisfied, and thus 0 ≤ s(t, N(t)) ≤ 1 also holds, which ensures that the birth and death rates remain positive. Under the deterministic approximation, the total population size thus changes according to

dNdt=rN(t)(1N(t)K)N(t)s(t,N(t)).2.4

(d) Extinction probability without adaptation

In the absence of adaptation due to de novo mutations, the population size is determined based on the initial mean fitness (i.e. mean phenotypic distance from the optimum), w¯(t)=w0, and approaches the following equilibrium:

N=Kr(w0+r1).2.5

If N* is close to 0, the population may become extinct by chance due to demographic stochasticity. To quantify the risk of extinction, we use the probability of loss in a general birth–death model starting from an initial state with N0 individuals:

ρ0,N0=1k=0N01γkk=0γk,2.6
where γi=j=1id(N(t)=j)/b(N(t)=j) (for i > 0) and γ0 = 1 (e.g. [37]).

(e) Simulation methods

We performed simulations by initially sub-dividing a genetically monomorphic species into two finite populations, each comprising K individuals with no migration between them (i.e. in allopatry). These populations are genetically identical initially, having just separated geographically.

To perform the evolutionary simulation with explicit population dynamics using SLiM3.3 [32], the demography described in the birth–death process above needs to be approximated using the Wright–Fisher model with discrete generations. In a single population at time t, mating pairs are formed by randomly drawing individuals with replacement according to their relative fitness, and each pair produces a single offspring. This process is repeated N(t + 1) times, where N(t + 1) is the total number of offspring in the next generation, t + 1. To determine N(t + 1), we consider the following stochastic differential equation to calculate the expected number of offspring in the next generation: ΔN=M(N(t))Δt+V(N(t))Δω, where Δω is a random variable with mean 0 and variance equal to Δt. The stochastic dynamics can be handled by considering the systematic change M(N(t))=limΔt0E[ΔN|N(t)]/Δt and the variance generation rate V(N(t))=limΔt0Var[ΔN|N(t)]/Δt as follows:

M(N(t))=rN(t)(1N(t)K)N(t)s(t,N(t))andV(N(t))=rN(t)(1N(t)K)+N(t)(s(t,N(t))+2β).2.7

Thus, we can calculate the total number of offspring in the next generation: N(t+1)=N(t)+ΔN with Δt set to 1.

We followed the populations for 4000 generations, unless otherwise stated (as illustrated in electronic supplementary material, figure S1, 4000 generations was sufficient to reach equilibrium). For the analysis of hybrid fitness reduction, we focus on simulations in which the final fitness of both parental populations is higher than 0.95 to concentrate on those populations that had successfully adapted and avoid asymmetries in parental fitness that can affect hybrid fitness (i.e. these were considered to be still at risk of extinction and not included, see electronic supplementary material, table S1 for the number of simulations excluded). Note that our individual-based simulations allow polymorphic sites to contribute to adaptation and hybrid breakdown and allow for transient overdominance, which is often found in Fisher's geometric model with diploids (e.g. [38]). For example, the fixation probability of a mutation can depend on the individual in which it arises in a polymorphic population, which is not true in models that simply track the series of fixed substitutions.

(f) Phenotypic parallelism

Two parental populations on independent evolutionary paths may undergo adaptation by using some common genetic loci. To quantify the degree of parallel evolution between two populations, we define ‘phenotypic parallelism’ as the fraction of the distance from the ancestral population to the optimum spanned by shared mutations (figure 1d). At the end of each simulation, for alleles with a frequency of 95% or more in each population, we calculate phenotypic parallelism. In the limit of the infinite-sites model, it is impossible to share the exact same mutation between populations, and thus phenotypic parallelism is zero. However, if two populations complete their adaptation with exactly the same set of genetic variants (with the same phenotypic effects), phenotypic parallelism is 1. In rare cases, phenotypic parallelism can be negative, which means that only maladaptive mutations (those that move away from the adaptive peak) are shared among the populations.

3. Results

Following an environmental change in the optimal phenotype from {0, 0, 0, 0} to {2, 0, 0, 0} at time t = 0, the populations were initially poorly adapted, with mean fitness of w0 ≈ 0.018 for the parameters assumed (table 1; with c = 0 for full pleiotropy). Without any subsequent adaptation, population size was expected to decline rapidly from the carrying capacity (K = 2000) towards a new equilibrium given by equation (2.5) (e.g. N* = 36 when r = 1). At such small population sizes, stochastic extinction is likely to occur, which was frequently observed when r was roughly lower than 1 (electronic supplementary material, figure S2). Evolutionary rescue could, however, occur if sufficiently large adaptive mutations arose soon enough (electronic supplementary material, figure S1). As expected, the population was more likely to be rescued when there were more sites under selection (electronic supplementary material, figure S2). Hereafter, we consider three scenarios with three different extinction risks: high extinction risk (r = 0.8), moderate extinction risk (r = 1.0) and low extinction risk (r = 1.4). See electronic supplementary material, figure S2 for examples of extinction frequencies associated with these r values.

As expected, populations at higher risk of extinction fixed larger effect adaptive mutations, given that the population survived (figure 2a). With a large number of mutable sites within the genome (L = 10 000, which we consider to be close to an ‘infinite-sites model’), we find that phenotypic parallelism is almost always zero except for a small fraction of populations that incorporated the same mutations by chance (electronic supplementary material, table S1), and the fixation of larger effect mutations leads to greater reproductive isolation in populations that face a higher extinction risk, consistent with the phoenix hypothesis (figure 2e).

Figure 2.

Figure 2. Characteristics of mutations that accumulate during the process of adaptation and the reduction of F2 fitness (left panels for L = 10 000 and right panels for L = 50). (a,b) Surviving populations at higher risk of extinction are more likely to adapt via large-effect mutations on the first phenotypic axis (the focal trait under environmental change). Here, extinction risk is low, moderate or high for populations with intrinsic growth rates of r = 1.4, 1 and 0.8, respectively. (c) With a large number of mutable sites within the genome, parallel adaptation rarely occurs. (d) However, when the number of loci available for adaptation is finite, parallel mutations are more likely under a high extinction risk. (e) F2 fitness decreases monotonically as extinction risk increases, which is expected in the phoenix hypothesis. (f) In the case of L = 50, the net result is that populations facing an intermediate risk of extinction exhibit the lowest F2 hybrid fitness. The black dots and background violin plots represent the mean ± s.d. calculated from 1000 simulations. The 1000 simulations were divided into cases with zero and non-zero phenotypic parallelism values, shown in blue and red, respectively. Likewise, solid points and lines indicate mean ± s.d. in right panels. The small dots are the values for each simulation run, and for the F2 fitness value, it is the average value of 2000 F2 hybrid individuals created. Percentage differences of the mean values from left to right data are shown at the top of each panel. Parameters used were the default values (table 1), with c = 0. Asterisks indicate the significance level for Kruskal–Wallis tests. ***p < 0.001; **p < 0.01; *p < 0.05; NS, not significant. (Online version in colour.)

While such large-effect mutations are expected to increase the extent of reproductive isolation among populations, with a limited number of loci available for adaptation (L = 50), the overall pattern of phenotypic parallelism also becomes higher as extinction risk increases (figure 2d), because surviving populations are more likely to re-use the same ‘rescue’ mutations. Note that, if we measure genetic parallelism by counting the number of shared mutations rather than phenotypic parallelism, the number is slightly smaller when extinction risk increases even if the distance spanned by the shared mutations is expected to be larger (electronic supplementary material, table S1). Because of these contrasting effects between phenotypic effect size and parallelism of the mutations, we find that hybrid fitness is lowest for populations at intermediate risk of extinction (figure 2f), with different large-effect mutations incorporated in different populations. By contrast, when there is a low extinction risk, more mutations that accumulate during adaptation are unique, but those mutations are relatively small-effect and do not promote speciation. Therefore, moderate extinction risk results in the highest possibility of speciation due to the balance between mutation size effect and mutation uniqueness. This pattern predicts a greater degree of reproductive isolation among populations experiencing a high extinction risk, but not so high that rescue frequently requires the same mutation(s).

We next asked if these results were robust to changing assumptions. First, we increased the number of sites (from L = 50 to 100, 500 and 1000). Similar results were obtained, except that phenotypic parallelism did not display a monotonic decline with increasing probability of population survival (i.e. with increasing r). Specifically, populations at low and at high extinction risk showed significantly more phenotypic parallelism than those populations at moderate risk of extinction (electronic supplementary material, figure S3 and S4, focusing on red points where some phenotypic parallelism was observed). The phenotypic parallelism observed at higher growth rates (r) was only slightly higher and likely reflects the maintenance of more genetic variation when the population size remained larger, favouring the same highly beneficial mutations especially when the number of loci was greater. Consistent with this interpretation, fewer mutations fixed during the course of adaptation (electronic supplementary material, figure S5) when the extinction rate was high (fixing fewer and larger mutations) and when it was low (higher competition for the most beneficial mutations, see electronic supplementary material, figure S6c), compared to populations at moderate extinction risk.

We also explored changing the Fisher's geometric model to relax the assumption of complete pleiotropy, i.e. that each mutation affects each phenotype independently. Specifically, we allowed mutations to have a primary effect on a specific axis (no pleiotropy with probability 1 − c), as well as a randomly drawn mutational effect with complete pleiotropy. With moderate pleiotropy (c = 0.5), we again observed that a moderate extinction risk results in more extensive reproductive isolation (electronic supplementary material, figure S3 for L = 100 and electronic supplementary material, figure S7 for L = 50). The number of fixed mutations during adaptation also showed qualitatively similar results (electronic supplementary material, figure S5). These results support our conclusion that moderate extinction risks cause the highest chance of speciation, whether reproductive isolation results primarily from pleiotropic side effects (c = 0) or from overshooting effects (i.e. a combination of large-effect mutations on the focal phenotypic axis overshoots the optimum in hybrids) (c = 0.5). Furthermore, with no pleiotropy (c = 1), the average degree of parallel evolution between the two populations is high but similar for extinction risks considered (electronic supplementary material, figure S7d), and we observed greater loss in hybrid fitness among those populations facing the greatest risk of extinction, as expected by the phoenix hypothesis (electronic supplementary material, figure S7f).

Finally, we explored another scenario where there is just an increase in the strength of selection but no change in the optima (electronic supplementary material, text S1). Extinction was possible in these simulations, as long as mutation rates were high enough and the stringency of selection increased enough that few individuals could survive in the new environment. In this case, again, adaptation in different populations often involved different sets of compensatory mutations giving rise to parental populations very close to the optimum. Crosses between these populations, however, often led to low hybrid fitness, especially when extinction risk was high, again supporting the phoenix hypothesis (electronic supplementary material, figure S8).

In summary, whether we see the greatest degree of reproductive isolation in populations facing the highest extinction risk (the original phoenix hypothesis) or in populations facing an intermediate extinction risk (so that mutations did not fix in parallel too often) can vary depending on the genetic architecture underlying adaptation, including the number of loci and the degree of pleiotropy.

4. Discussion

Allopatric speciation cannot be observed unless geographically isolated populations persist for the thousands or millions of years required for reproductive isolation to build [39]. During this period, environmental change can reduce the size of populations, increasing the likelihood of extinction due to demographic stochasticity [40]. In extreme cases, populations may be so maladapted that they decline in size to extinction unless beneficial mutations arise and lead to evolutionary rescue [25]. In this study, we argue that populations that persist through such periods of rapid adaptation are more likely to accumulate the large-effect genetic changes likely to cause reproductive isolation. We call this process the phoenix hypothesis of speciation, because new species arise from the near extinction of populations that have undergone evolutionary rescue.

Reproductive isolation, as measured by a reduction in fitness of hybrids, is strongest when adaptation involves different large-effect mutations in different parental populations. Using Fisher's geometric model with the assumption of an infinite-sites model and constant population size, Yamaguchi & Otto [7] showed that large-effect mutations can promote strong reproductive isolation even when populations are adapting to similar environmental changes. This previous work did not, however, consider population dynamics explicitly, and so did not explore the risk of extinction. Here we find an increase in the fixation of large-effect adaptive mutations in populations at high extinction risk and a greater reproductive isolation among those populations. Thus, the phoenix hypothesis was supported under the infinite-sites model. The study by Yamaguchi & Otto [7] also did not consider the potential for parallel adaptation, which is high when few mutations can rescue a population (figure 1). Parallel genetic evolution increases with the strength of natural selection and is particularly likely to occur for genes with large phenotypic effects [41]. We find that populations that persist despite a high risk of extinction are more likely to accumulate large-effect adaptations, which more often generate reproductive isolation, but these adaptations are also more likely to be parallel, reducing the amount of reproductive isolation (figure 2). Given these contrasting effects, we find that populations at an intermediate extinction risk develop the most reproductive isolation (figure 2).

Because the magnitude of these contrasting effects depends on the parameters, the precise conditions most favourable to the phoenix hypothesis may not be simple (figure 1a). The number of loci underlying adaptation varies from a few to hundreds depending on the species and environments encountered (as seen in experimental evolution with yeast, e.g. [42]; see also [43] for review). Because most genomic studies focus on detecting specific large-effect mutations or genomic regions, the overall number of mutations involved can be larger when considering many small-effect mutations and compensating mutations. In addition, Conte et al. [29] revealed that the probability of gene reuse (parallel change) was surprisingly high when the age of the common ancestor of compared taxa was young. This may imply that a harsh environment change to closely related populations may often promote parallel adaptations in nature. Furthermore, we have observed non-monotonic trends for both phenotypic parallelism (electronic supplementary material, figures S3 and S4) and the number of fixed mutations at the end of adaptation (electronic supplementary material, figure S5), results that might reflect stronger competitive replacement when many alternative mutations co-occur in large populations (electronic supplementary material, figure S6).

Changing population size is a key feature of a related model of speciation. Founder effect speciation considers the case where a small number of individuals colonize a new site, with genetic drift causing considerable genetic differentiation from the ancestral population [44,45]. The founder effect can cause a rapid peak shift through an extreme bottleneck, but this differs substantially from our model, which focuses on adaptation to a single peak, with environmental change driving genetic change rather than drift.

Although we mainly focused on large-effect beneficial mutations and adaptation, fixation of weakly deleterious mutations by drift is likely to occur in a small population given enough time and a sufficiently high mutation rate. Deleterious mutations can also contribute to hybrid fitness breakdown because a certain fraction of the progeny will be expected to carry more deleterious mutations than either parent [46,47]. On the other hand, masking of these mutations in diploid hybrids can offset reproductive isolation, leading to hybrid vigour [48,49]. Although we did not include unconditionally deleterious mutations, there is a chance that linked deleterious alleles hitchhike to fixation along with beneficial rescuing alleles [50,51]. Such ‘undesirable hitchhikers' could also contribute to fitness differences between populations. Future work that includes such hitchhiking would be valuable to determine the extent to which the fixation of different deleterious mutations in different parental populations facilitates introgression and counteracts the phoenix hypothesis.

It would also be worthwhile exploring the effect of standing genetic variation or migration in future studies. Standing genetic variation enables the parental populations to incorporate the same alleles when the environment changes, increasing phenotypic parallelism [31,52]. Similarly, migration between parental populations allows the spread and fixation of the same favourable alleles. A comparative genome study of nine- and three-spined stickleback populations revealed that small effective population sizes and restricted gene flow limit the potential for parallel local adaptation [27]. Furthermore, the history of a population can impact reproductive isolation, with genetic surfing at the range of a population leading to ‘allele surfing’ to high frequency of deleterious recessive alleles, which again facilitates introgression when secondary contact occurs [49].

5. Conclusion

The phoenix hypothesis predicts that populations are more likely to speciate if they have faced a history of environmental change so challenging that rapid adaptation was required to avoid extinction. While this hypothesis predicts a coupling between extinction risk and speciation, we emphasize that the new species formed by this process are not equivalent to the many lineages lost to extinction. Each lineage lost carries with it a long history of genetic adaptations. It has been estimated that new species arise, on average, at a rate of once per 2 Myr [53], and extinction eliminates all of the adaptations that allowed that species to survive and thrive since its last common ancestor with extant species. By contrast, newly formed species under the phoenix hypothesis can have minimal genetic differences, which makes them more prone to species collapse through introgression or competitive exclusion in the future [54].

With the accumulation of studies highlighting the role of large-effect mutations and adaptation, we find that rapid environmental change typically facilitates the evolution of reproductive isolation despite the higher extinction risk. Those lucky few populations that survive severe environmental change by evolving might, like the mythical phoenix, rise from the ashes reborn, taking the first baby steps to form a new species.

Data accessibility

All code necessary to repeat the simulation described in this study have been made available. SLiM source codes of our model for speciation dynamics and data that we used to produce figures are available from the Dryad Digital Repository [55].

The simulation results are provided in the electronic supplementary material [56].

Authors' contributions

R.Y.: conceptualization, formal analysis, funding acquisition, investigation, methodology, validation, visualization, writing—original draft and writing—review and editing; B.W.: formal analysis, investigation, visualization and writing—review and editing; S.P.O.: conceptualization, funding acquisition, methodology, supervision, visualization, writing—original draft and writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This work was supported by funding from the JSPS KAKENHI (21K15160 and Overseas Research Fellowships to R.Y.) and NSERC Discovery Grant RGPIN- 2022-03726 to S.P.O.

Acknowledgements

We would like to thank the Otto laboratory, Amy Forsythe, Dolph Schluter, Kiarash Bahari and two anonymous reviewers for their helpful feedback on previous versions of the manuscript.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.6261899.

Published by the Royal Society. All rights reserved.

References