The propagation of admixture-derived adaptive radiation potential

Adaptive radiations (ARs) frequently show remarkable repeatability where single lineages undergo multiple independent episodes of AR in distant places and long-separate time points. Genetic variation generated through hybridization between distantly related lineages can promote AR. This mechanism, however, requires rare coincidence in space and time between a hybridization event and opening of ecological opportunity, because hybridization generates large genetic variation only locally and it will persist only for a short period. Hence, hybridization seems unlikely to explain recurrent AR in the same lineage. Contrary to these expectations, our evolutionary computer simulations demonstrate that admixture variation can geographically spread and persist for long periods if the hybrid population becomes separated into isolated sub-lineages. Subsequent secondary hybridization of some of these can reestablish genetic polymorphisms from the ancestral hybridization in places far from the birthplace of the hybrid clade and long after the ancestral hybridization event. Consequently, simulations revealed conditions where exceptional genetic variation, once generated through a rare hybridization event, can facilitate multiple ARs exploiting ecological opportunities available at distant points in time and space.

The growth performance Pi of an individual i in a habitat Hi is given by exp{−φ(xi − xopt_Hi) 2 }, where φ is a parameter controlling the strength of natural selection. The growth performance Pi has two-fold effects on survival fitness. First, we assumed that individuals with higher growth performances (locally better-adapted individuals) have stronger competitive influence on other individuals and higher tolerance to competition. Additionally, second, we assumed that basal survival rate in the absence of competition increases proportional to the growth performance. Owing to the latter assumption, colonization of an unoccupied habitat required a certain degree of phenotypic specialization to the habitat unless natural selection is weak. A previous simulation study suggested that the effect of hybridization to promote adaptive radiation will be most pronounced under such conditions (i.e. when optimal phenotypes of distinct ecological niches are isolated by fitness valleys of maladaptive intermediate phenotypes and thus invasion of novel ecological niches requires discontinuous phenotypic change towards a novel optimal value). Such a situation would be realistic when a lineage adaptively radiates into novel ecological niches that are highly dissimilar to each other. Following these assumptions, we assumed that survival probability of an individual i is given by = , where MHi is the number of individuals in the habitat Hi. In the absence of competitors (MHi = 0), the denominator becomes 1 and wi equals to the growth performance Pi.
Negative effect of competition on survival rate increases with the sum of growth performance of all competitors and decreases with KHi(Pi), which is the tolerance of the individual i to competition. We assumed KHi(Pi) = PiNf/(f − 2) so that the tolerance to competition linearly increases with growth performance. Parameter f is per-capita female fecundity (f > 2) and N is a constant that determines the carrying capacity. When all individuals have the optimal trait value (and thus Pi = 1), carrying capacity for adults in a habitat becomes equal to N, because the demographic dynamics within a habitat approaches the classical Beverton-Holt model with carrying capacity N (for adult) when Pi = 1. Then, each surviving female selects a mating partner randomly from males living in the same habitat of the same patch. Females mates only once, whereas males can mate multiple times. The number of offspring per mating pair follows Po (f). Sex of newborn individuals is randomly assigned. Newborn individuals move to a neighboring patch with probability mP and to a randomly selected habitat in the same patch with probability mH.

Appendix S2. Conditions under which hybridization promotes adaptive radiation
Since we aimed to simulate repeated occurrence of hybridization-enabled adaptive radiation, we intended to use parameter sets with which hybridization is essential for rapid adaptive radiation.
Therefore, prior to investigating repeated occurrence of adaptive radiations, we explored conditions under which hybridization promotes a single episode of adaptive radiation using simulations with a single patch including five different habitats (H1,2, …,5). In this simulation, we allowed spontaneous mutations both before and after hybridization; thus, both spontaneous mutations and hybridizationinduced genetic variation can contribute to adaptive radiation. Each simulation was started by introducing 100 clone individuals of each two allopatrically evolved lineages into the patch (see Simulation of hybridization in Methods for more details) and was continued for 100,000 generations after the colonization event. We performed simulations systematically varying the strength of natural selection in each habitat φ and the length of the period where two parental lineages had independently evolved in allopatry before they hybridized (T0 parameter combinations that we explored, rapid adaptive radiation within 5000 generations was extremely unlikely to occur when hybridizing lineages were not genetically differentiated (Fig. S2 ac; T0 = 0), which is in accord with the results of Kagawa & Takimoto (2018). In contrast, adaptive radiation occurred within 5000 generations when hybridizing lineages were highly genetically differentiated ( Fig. S2 a-c; large values of T0). With parameter conditions that were not conducive to adaptive radiation within 5000 generations, species diversification was rarely observed even after 100,000 generations ( Fig. S2 d-f). These results confirmed that adaptive radiation by spontaneous mutations alone without hybridization between genetically differentiated lineages was very unlikely with parameter values that we used. Adaptive radiation was not likely to start without hybridization producing initial standing genetic variation in our model, because we assumed that ecological niches of distinct habitats were discrete and that phenotypic effects of single mutations were not very large.
As ecological niches were discrete, establishment of a population in a novel empty habitat, which is a prerequisite for natural selection to operate to cause local adaptation, required a certain level of phenotypic specialization to the habitat unless stabilizing natural selection in each habitat was extremely weak. Such phenotypes specializing to a novel habitat was not likely to evolve by accumulation of spontaneous mutations with small effects in populations inhabiting different habitats.
Accordingly, standing genetic variation, that hybridization between genetically differentiated lineages can efficiently generate, played key roles for promoting colonization of new habitats in the early stage of adaptive radiation. Kagawa & Takimoto (2018) provides more comprehensive comparisons between cases with and without hybridization regarding the likelihood of adaptive radiation in face of an ecological opportunity.

Appendix S3. Empirical basis of the default parameter values
Default values for most of model parameters were determined as in Kagawa & Takimoto (2018). The number and length of chromosomes (n and l) and the mutation rate per locus (μ) were set to n = 15, l = 2×10 8 base-pairs (bps), and μ = 10 -5 /generation, respectively, so that they fall in realistic ranges in animals and plants. Recombination rate r was set to 10 -8 /bp/generation, because 100 centi-Morgan typically corresponds to about 10 8 base-pairs in human genome (Lynch & Walsh 1998). The length of single loci was set to 5000 bps to be in the range of typical values (Rafalski & Morgante 2004) and to be consistent with the fact that mutation rate per locus is 100 to 10000 times higher than the mutation rate per base-pair (Lynch & Walsh 1998). For the number of loci that influence a quantitative trait, we used L = 400 instead of larger values (values in between 1000 and 10000) used in Kagawa & Takimoto (2018). Although the ratio of typical mutation rate per quantitative trait (10 -1 /generation) and per locus (10 -6 ~ 10 -5 /generation) indicates that each trait should be influenced by more than 10,000 potential loci (Lynch & Walsh 1998) were chosen to simulate long-distant range expansion that takes many generations. We also performed simulations varying k. In the temporally repeated AR scenario, we set the default migration rate between sub-populations during the refugial phase (mR) to 0 and performed simulations varying the value of mR. The length of allopatric evolution of parental lineages controls the amount of genetic variation generated through hybridization (Kagawa & Takimoto 2018). We selected T0 = 200,000 as the default value so that hybridization can produce sufficiently large genetic/phenotypic variation to cause adaptive radiation (Fig. S2).

Appendix S4. Methods to count the number of incipient species and genotypically distinct clusters of individuals
Method to count the number of incipient species We counted the number of genetically isolated and phenotypically distinct incipient species by using the method proposed in the Appendix S4 of Kagawa & Takimoto (2018). In order to apply the method, we recorded trait values of all individuals who had survived to maturation in two successive generations g and g + 1. We call individuals from the generation g + 1 "offspring-individuals". To count the number of species that exist in a focal patch, we sorted all individuals in the patch into 10 clusters based on phenotypic similarity using k-means algorithm. Then, we merged clusters that are connected by gene flow. The strength of gene flow was measured based on records of parent-offspring relationships in the two successive generations. Clusters i and j were merged when: more than t % of parents of offspring-individuals from the cluster-i belonged to another cluster j, and more than t % of parents of offspring-individuals from the cluster-j belonged to the cluster i. Additionally, if more than 50% of parents of offspring-individuals from the cluster-i belonged to another cluster j, cluster i was merged into the cluster j. This procedure was repeated until merging of clusters stopped. Finally, each cluster was counted as a species only if more than 90% of parents of offspring-individuals in the cluster belonged to the same cluster, because otherwise the cluster is likely to be a temporal hybrid swarm.
Also, clusters that contained less than 20 individuals were not counted as species because such small clusters are probably ephemeral. For the threshold level of gene flow to merge clusters, we used t = 2%. With this threshold value, ecological differentiation between populations in parapatric habitats was necessary for the splitting of species. In our model, newborn individuals migrate to a randomly selected habitat at a rate 0.2 every generation. Thus, in patches with five parapatric habitats, about 4% of individuals in a habitat should be immigrants from another certain habitat when population sizes in all habitats are equivalent. Therefore, spatial isolation between parapatric habitats is not enough for splitting species under our criterion. Ecological differentiation between populations could enable speciation by strengthening natural selection against immigrants in their non-native habitats (i.e. immigrant inviability), which operates as a barrier to gene flow. We note that our definition of incipient species does not indicate good species with permanent reproductive isolation. Since our model does not incorporate non-random mating and post-zygotic reproductive isolation, incipient ecological speciation of our model can collapse if spatial isolation between habitats or divergent natural selection are removed.

Method to count the number of genotypically distinct clusters of individuals
For the temporally repeated AR scenario, we counted the number of genotypically distinct clusters of individuals at the end of simulation (after 30000 generations from the hybridization). We analyzed a subset of the population consisting of 100 randomly sampled individuals to reduce the computational time. To analyze genotypic differences among individuals, we listed all polymorphic nucleotides among genomes of all sampled individuals. For each polymorphic nucleotide, genotype was expressed by the number of ancestral and derived nucleotide (0 and 2 indicate homozygotes of ancestral and derived nucleotide whereas 1 indicates heterozygote). Genome-wide genotype of each individual was expressed as a vector of which i-th element represents the number of ancestral and derived alleles at the i-th polymorphic nucleotide. We applied hierarchical clustering of individuals based on Euclidian distance between genome-wide genotype vectors. The number of genotypically distinct clusters of individuals was examined by counting clusters that are distant from each other by at least the threshold value DT. Since we aimed to find genetically distinct populations rather than the statistically optimal clustering, we used an arbitrary constant threshold value DT = 30, which roughly captured the number of genetically distinct populations ignoring small-scale population genetic structures generated by a few allelic polymorphisms. For the spatially repeated AR scenario, we counted the number of genotypically distinct clusters of individuals that were found in (i) only region 2 and (ii) both region 1 and 2 at the end of simulation (after 5000 generations from the hybridization). In this purpose, we first sampled 100 individuals randomly from each of region 1 and 2. Using the method described above, we identified genotypically distinct clusters in all 200 individuals sampled. Then, we counted the number of clusters that contain (i) only individuals from the region 2 and (ii) individuals from both regions. The algorithms to count genotypic clusters were implemented with the R language.

Appendix S5. A review of potential examples where the maintenance of hybridization-induced evolvability by temporal isolation and subsequent admixture of sub-lineages might have promoted recurrent adaptive radiation.
Hybridization generating high evolvability, followed by the maintenance of elevated evolvability by temporal isolation and secondary admixture of sub-lineages could have contributed to the repeated occurrence of adaptive radiations in archipelagos, such as Hawaiian, Galapagos, and Canary Islands.
Colonization of oceanic island from mainland will induce strong population bottleneck causing the lack of genetic variation in the island population. Thus, inter-and intra-specific hybridization prior to or subsequent to the colonization of island is thought to be an important source of genetic variation for adaptive radiations in islands ( Additionally, natural disturbances, such as volcanic eruption, storms, and see-level fluctuations, could occasionally cause a collapse of island biota, followed by re-colonization and a renewal of biota (Price & Clague 2002;Schneider et al. 2005). These characteristics of archipelago would be conducive to repeated occurrence of adaptive radiation if clades do not lose evolutionary potential in the course of adaptative evolution and progressive colonization of islands. Long-term observational studies of Galapagos finches suggest that temporal isolation and subsequent admixture of populations contributes to prevent the exhausting of standing genetic variation in local populations adapting to changing environment (Grant & Grant 2008). Moreover, in some other island taxa, populations in new islands show higher standing genetic variation compared to populations of their relatives in older islands, due to recent secondary hybridization events (Stacy et al. 2014;De Busschere et al. 2015;Hendrickx et al. 2015). Taken     was extremely unlikely to occur (T0 = 0 in panels a-c). In such cases, species diversification rarely observed even after long-term evolution from the colonization event (T0 = 0 in panels d-f), especially when stabilizing natural selection within each habitat was strong (φ ≥ 0.1). With hybridization between genetically differentiated lineages (large T0), in contrast, adaptive radiation could occur within 5000 generations unless natural selection was extremely strong. Extremely strong stabilizing natural selection within each habitat made adaptive radiation difficult by impeding invasion of new unoccupied habitats. This effect was owing to our model assumption that individuals with maladaptive phenotypes in their local environment have low survival rates even in the absence of competitors.
When natural selection was very weak (φ = 0.025), on the other hand, the number of species remained small despite that wide range of trait values have evolved. This was because immigrant inviability, which contributes to reduce gene flow between ecologically differentiated populations, was mitigated when stabilizing natural selection in each habitat was weak. These results confirmed that adaptive radiation was unlikely to occur without hybridization between genetically differentiated lineages in our model with parameter values that we used.  (Table 1). In all simulation runs, there were no genotypic clusters shared between region 1 and 2. Consistently, the number of genotypic clusters that were endemic in region 2 was equivalent to the number of reproductively isolated species in region 2. These results together imply independent evolution of phenotypically diverse species in two regions. As speciation in our model requires ecological differentiation, the number of species in region 2 was proportional to the trait range in region 2. The scarcity of species shared between two regions is reasonable in our model for three reasons. First, as a population in corridor is established by a subset of the hybrid swarm in region 1 before the loss of hybridization-induced standing genetic variation, species in region 1 and sublineage(s) in corridor(s) will fix different sets of genes. Owing to this founder effect, colonizers of region 2 will be genotypically distinct from species in region 1. Although gene flow between two regions through corridor populations can homogenize genotypes of species using the habitat 1, such genotypic homogenization has not been observed within 5000 generations from hybridization. Second, when secondary admixture between the genetically isolated sub-lineages occurs in region 2, increased genotypic diversity will lead to new genotypic combinations to be recruited into the new species in region 2. Third, even when the corridor environment imposes weak or no selection and potentially allows coexistence of multiple species, direct colonization of multiple species is difficult in our model because weak selection in the corridor will weaken reproductive isolation between species such that they will merge in the course of their range expansion. Figure S4. Effects of the carrying capacity of single habitats in corridor patch NC, and the natural selection strength within corridor patches φC in the spatially repeated AR scenario. We performed simulations with (a) a single environmentally homogeneous corridor, (b) two geographically isolated environmentally homogeneous corridors, and (c) a single corridor containing habitats with two distinct environments. Other parameters were set to default values (Table 1). Figure S5. An example of failure of geographically repeated adaptive radiations with two geographically isolated corridors. Trait value x of all individuals in two regions and corridor patches CPi,9,CPi,14,and CPi,20 (i = 1,2) are shown for each generation. Individuals who did and did not survive to participated in reproduction are shown in green and grey respectively. In this run, spread of hybridization-induced genetic variation from region 1 to 2 was prevented probably by formation of a hybrid zone in the corridor 2. Hybrid zone could be formed on a corridor when a lineage from one corridor reaches the region 2 earlier than the lineage of the other corridor and then migration from region 2 to the other corridor leads to collision of two lineages in corridor. In this simulation example, we observed high phenotypic variation at around the corridor patch CP2,14. As novel hybrid phenotypes are selected against in corridor patches, genetic variation of hybrids would not easily spread out from the hybrid zone. Consequently, adaptive radiation in region 2 did not occur in 5000 generations from the hybridization event. Parameter values: k = 20 and other parameters were set to default values (Table 1). Figure S6. The spatially repeated adaptive radiation mediated by two geographically isolated corridors was diminished if there was a time lag in the introduction of two parental lineages. In this simulation, 100 individuals of parental lineages 1 and 2 were introduced to region 1 at the generation 0 and 100, respectively. We simulated cases with (a) a single environmentally homogeneous corridor, (b) two geographically isolated environmentally homogeneous corridors, and (c) a single corridor containing habitats with two distinct environments. Additionally, we systematically varied the length of corridor k, and the natural selection strength within corridor patches φC. Other parameters were set to default values ( Table 1). Comparison of Figs 3b and S6b indicates that the time lag in introduction of two parental lineages suppresses repeated adaptive radiation mediated by two expansion corridors. In such situation, the first lineage spread across the entire system of connected patches before the introduction of the other parental lineage. Then, the arrival of the second lineage in region 1 generated a hybrid population and facilitated an adaptive radiation in region 1, which then led to formation of a hybrid-zone between the hybrid-lineage and lineage 1 in the corridors. This impeded the spread of genes from the second parental lineage toward the second region, because most of hybrid genotypes with associated phenotypes were less well adapted to the habitat in corridor patches compared to genotypes of the resident lineage 1 population. However, our additional simulations revealed that the presence of two isolated paths for range expansion still slightly elevated the likelihood of repeated adaptive radiation in long-term even if there was a time lag in the introduction of two parental lineages (Fig.   S7). Simulation results confirmed that the feasibility of repeated adaptive radiation was higher with two geographically isolated corridors than with only a single corridor. However, the time lag in the introduction of two parental lineages delayed the start of adaptive radiation in region 2 compared to the case of simultaneous introduction of two parental lineages (Fig. 3b). Additionally, the effect of the presence of two isolated corridors to promote repeated adaptive radiation was not as strong as in the case of simultaneous introduction of two parental lineages (Fig. 3). Figure S8. Effect of one-way migration in the spatially repeated AR scenario with a single environmentally homogeneous corridor. In simulations shown here, only migration in the direction from region 1 to 2 is allowed (i.e. migration from CPi,j to CPi,j+1 is allowed whereas migration from CPi,j to CPi,j-1 is forbidden). We performed simulations systematically varying the length of the corridor k and the strength of natural selection in corridor patches φC. Other parameters were set to default values (Table 1). Simulation results demonstrate that repeated adaptive radiations at geographically distant areas is possible even with a single environmentally homogeneous corridor if migration is one-way.  Table   S1. In the region 2, secondary admixture of sub-lineages from two isolated corridors generated various phenotypes including those with low fitness in the region 1, which enabled rapid evolution of four incipient species specialized to novel habitats that did not exist in the region 1.  Table S1. In the region 1, a hybrid population rapidly diverged into five incipient species with distinct phenotypes of the ecological trait 1 but with the same phenotype of the ecological trait 2. Two of them expanded their range towards the region 2 through the corridor using two distinct habitats of the corridor. Subsequently, two incipient species, which were genetically isolated in region 1 and corridor by ecological natural selection on the trait 1, admixed in the region 2.
This reestablished phenotypic variation in not only the trait 1 but also the trait 2 through transgressive segregation. This enabled adaptive radiation in the region 2 in response to divergent natural selection on the trait 2. This result supports that incipient speciation by ecological divergent selection on a single trait axis can maintain potential genetic variation for other trait axes subjected to stabilizing selection.
Secondary admixture between such incipient species can promote rapid adaptive diversification along novel trait axes that have never diversified in the previous radiation.  Table S1.    (Table 1).

Figure S13.
Simulation results with an alternative version of the temporally repeated AR scenario. In this scenario, we assumed that the environmental change led to a loss of divergent natural selection by causing shifts of optimal trait values in five habitats (xopt_H for H = H1, 2, … 5) to the same value 0. To simulate such a dynamic change of optimal trait values, we assumed that the absolute value of optimal trait value in each habitat is not allowed to exceed a certain positive value dmax and that dmax gradually declines from 10 to 0 in the period of generations 10000 -10100. When the value of dmax became smaller than the optimal value in a habitat H (xopt_H), the value of xopt_H was updated to dmax.
Subsequently, environmental condition returns to the original state in the period of generations 20000 -20100. In this period, the value of dmax gradually increased from 0 to 10. When the value of dmax exceeded the current value of xopt_H and if the original value of xopt_H was greater than or equal to the current value of dmax, the value of xopt_H was updated to dmax. The strength of natural selection was not changed throughout the simulation (i.e. φR = φB). Therefore, refugial populations experienced a stabilizing natural selection with only a single optimal phenotype value 0. Simulations of this scenario with systematically varied parameter values revealed that the reestablishment of adaptive radiation could occur only when (i) migration rates between three patches during the refugial phase was low and (ii) divergent natural selection in the period before and after the refugial phase was moderately strong (panels a -c). These results are analogous to the results in our default version of simulation ( Fig. 5). When natural selection in periods before and after the refugial phase was very weak (φB = 0.025), some extent of phenotypic diversity could have been reestablished after the collapse of radiation, but the reestablishment of species diversity was unlikely. This is because phenotypic diversification was easy under weak natural selection that produce only shallow fitness valleys between optimal phenotypes of five habitat, while immigrant inviability, the only mechanism for reproductive isolation of species in our model, could not be strong enough to separate ecologically distinct populations when natural selection was weak. Panel (e) shows an example of evolutionary dynamics in which adaptive radiation was rapidly reestablished after the environmental condition returned to the original state (φB = 0.2, mR = 0). Other parameters were set to the default values (Table   1).  Table S1. An extreme phenotype (x ≈ 10) has fixed in all three refugia during the refugial phase. However, gradual environmental change during generations 10000 to 10100 caused speciation reversal (i.e. merge of formerly ecologically isolated incipient species into a hybrid swarm) in three refugia, which enabled genetic differentiation among them through fixation of different genotypes associated with similar phenotypes. After the renewal of environmental variation, secondary admixture of three refugial sublineages reestablished phenotypic variation through transgressive segregation, thereby enabled reestablishment of incipient species specialized to each habitat type. However, reestablishment of specialists for habitats H4 and H5 was less likely compared to the case with the intermediate refugial environment (Fig. S13), due to the increased distance between the optimal phenotype in the refugial environment and in habitat H4 and H5.