The role of recombination, niche‐specific gene pools and flexible genomes in the ecological speciation of bacteria

Abstract Bacteria diversify into genetic clusters analogous to those observed in sexual eukaryotes, but the definition of bacterial species is an ongoing problem. Recent work has focused on adaptation to distinct ecological niches as the main driver of clustering, but there remains debate about the role of recombination in that process. One view is that homologous recombination occurs too rarely for gene flow to constrain divergent selection. Another view is that homologous recombination is frequent enough in many bacterial populations that barriers to gene flow are needed to permit divergence. Niche‐specific gene pools have been proposed as a general mechanism to limit gene flow. We use theoretical models to evaluate additional hypotheses that evolving genetic architecture, specifically the effect sizes of genes and gene gain and loss, can limit gene flow between diverging populations. Our model predicts that (a) in the presence of gene flow and recombination, ecological divergence is concentrated in few loci of large effect and (b) high rates of gene flow plus recombination promote gene loss and favor the evolution of niche‐specific genes. The results show that changing genetic architecture and gene loss can facilitate ecological divergence, even without niche‐specific gene pools. We discuss these results in the context of recent studies of sympatric divergence in microbes.

however, does not help to uncover the evolutionary causes of diversification.
Recently, there have been growing efforts to frame species and speciation of Bacteria and Archaea in population genetic theory equivalent to theories for sexual eukaryotes. A major topic is whether recombination and gene flow between diverging populations plays an important role in bacterial speciation. Emphasizing the importance of ecological adaptation and clonal reproduction, Cohan and co-workers propose that adaptation to distinct ecological niches is the principal cause of bacterial species (which they call ecotypes).
Periodic selective sweeps reduce variation across the whole genome of a clonal niche-specialist while driving the accumulation of divergence between ecotypes (Cohan and Perry 2007), leading to a pattern of discrete genetic clusters. Recombination rates are proposed to be too low to produce maladapted "hybrid" genotypes and erode differences between ecotypes, and therefore no special mechanisms are needed to prevent gene flow between diverging populations.
Bacteria reproduce clonally but recombination occurs via several mechanisms that are decoupled from reproduction. Homologous recombination (HR) involves "copy and paste" of DNA from the environment that requires tracts of similar sequence at either end of the pasted DNA. Rates vary widely by several orders or magnitude, up to more than 10× the per nucleotide effect of mutation (Vos, Hesselman, Beek, Passel, & Eyre-Walker, 2015), which is sufficient to constrain multilocus divergence under neutral conditions (Fraser et al. 2009). Furthermore, several genomic studies report gene-specific sweeps during ecological divergence, which are only possible with recombination (Cordero and Polz, 2014;Shapiro et al. 2012), rejecting the ecotype model of divergence with strictly clonal evolution.
Akin to theories of speciation in sexual eukaryotes, it has been proposed therefore that recombination rates and gene flow between populations might be high enough in some bacteria to constrain multilocus divergence. Consequently, ecological speciation might require mechanisms to prevent gene transfer, analogous to reproductive isolating mechanisms in plants and animals (Polz et al. 2013).
Several mechanisms have been proposed for restricting HR between bacterial species. First, rates of HR decline as sequence divergence increases, which could act as a barrier to gene flow between populations (Fraser et al. 2009). This has been dismissed as a general explanation because HR depends less stringently on genetic similarity in Archaea (Polz et al. 2013), which nonetheless show similar patterns of genetic clustering (Cadillo-Quiroz et al., 2012). Second, there might be specific mechanisms to reduce gene exchange, for example, incompatible pheromone systems (Carrolo, Pinto, Melo-Cristino, & Ramirez, 2009) or incompatible restriction enzymes (Jeltsch, 2003).
While reported in some bacteria (e.g., Streptococcus), the general prevalence of such isolating mechanisms remains unclear.
A popular recent hypothesis is that restriction of gene transfer within habitat patches might restrict gene flow and permit ecological divergence, if populations adapting to distinct niches are only exposed to gene pools in their local habitat patch rather than globally across habitat types (Polz et al. 2013). This mechanism is analogous to habitat-based assortative mating in models of eukaryote speciation (Diehl & Bush, 1989). It relies on general features of the environment rather than specific genetic mechanisms and, hence, applies broadly across Bacteria and Archaea. It is also consistent with observed reductions of HR between ecologically distinct species (Polz et al. 2013). Melendrez et al. (2016) also found that genespecific selective sweeps occur within ecologically divergent species of Synechococcus, but in rebuttal of this hypothesis, they argued that separate gene pools were a consequence of ecological divergence rather than a cause. In theory, strong enough selection against recombinant genotypes could permit divergence even without nichespecific gene flow at the outset (Bürger, 2014;Gavrilets, 2004).
Resolving the role of recombination in bacterial speciation, and of niche-specific gene pools, in particular, depends critically on genetic architecture. Modelling of ecological divergence with gene flow in sexually reproducing organisms shows that gene flow has the greatest inhibitory effect on divergence when divergence involves many genes of small effects rather than few genes of large effect (Gavrilets, 2004;Kondrashov, 1986). With many genes of small effects, niche-specific assortative mating greatly enhances the progress of divergence (Diehl & Bush, 1989;Felsenstein, 1981). Friedman, Alm, and Shapiro (2013) extended this work to clonal bacteria with varying rates of HR. They showed that the constraining effect of HR between diverging populations increased with the number of loci, as in sexual models. They proposed reducing the number of loci needed for speciation, as well as niche-specific gene pools, as possible mechanisms to overcome this, without modelling them explicitly, and therefore, we construct models to explore those processes here.
We use simulations of a mathematical model to investigate the interplay between genetic architecture and niche-specific gene pools in facilitating ecological divergence of bacteria. Specifically, we ask whether evolving architecture provides an additional mechanism to aid ecological divergence. We investigate two ways in which genetic architecture can evolve; (a) through changes in the effects of loci on a trait affecting fitness, and (b) through changes in the number of loci through gene gain and loss. Previous models of sexual eukaryotes showed that evolving effect sizes of loci can protect locally adapted genotypes from the homogenizing effects of gene flow because the genetic basis of ecological traits becomes concentrated into fewer genes of large effect (Yeaman & Whitlock, 2011).
We further evaluate this idea to the case of bacteria with HR, which is not completely equivalent to meiotic recombination, as it tests a single maladapted gene region arriving by gene flow in an otherwise fully adapted genome.
A second aspect of genetic architecture that needs to be considered in modelling the effects of recombination on bacterial divergence is gene content. In contrast to eukaryotes, gene gain and loss is ubiquitous in bacteria. Among related strains, this divides the genome into a "flexible" part that differs between strains, and a "core" part that is shared (Cordero and Polz, 2014).
The rate of turnover sometimes exceeds thousands of gene gain and loss events per 1% amino acid divergence (Nowell, Green, Laue, & Sharp, 2014). Variation in gene content is often related to ecological traits. For example, within the marine cyanobacterium Prochlorococcus, both Atlantic and Pacific populations carried low-frequency genes as part of their "flexible" genomes, most of which were found on average in only one in four isolates. However, 29 genes associated with phosphorus uptake and metabolism had markedly different frequencies between populations; nearly absent in the Pacific, but close to fixation in the Atlantic (Coleman and Chisholm, 2010).
Clearly, gene gain is important for adaptive novelty and allowing bacteria to adapt to new conditions. Many accounts consider how flexible genomes aid colonization of new habitats, the spread of antibiotic resistance, and adaptation to fluctuating selection (Niehaus, Mitri, Fletcher, & Foster, 2015). Here, we focus instead on the effect of flexible genomes on divergence with gene flow, which has not previously been investigated in models of bacterial or eukaryote divergence. We hypothesize that niche-specific traits encoded by differences in gene content, that is, nonhomologous loci are protected against gene flow even if gene flow via HR occurs at sufficient rate to homogenize the core genome. Hence, both evolution of the effect sizes of loci and niche-specific genes might provide further isolating mechanisms against gene flow via HR between locally adapted populations, which we evaluate in models with or without niche-specific gene pools.

| Motivation
Our model of divergence between two subpopulations occupying distinct habitats but connected by DNA flow is inspired by the study of Shapiro et al. (2012). They reported on population divergence in Vibrio cyclitrophicus, a marine bacterium adapting to two niches in sympatry while undergoing high rates of HR. The ecotypes occupy different sized particles in the water column that are found in close proximity. The study found that genetic differentiation between the two forms (called L and S strains) was restricted to 11 discrete genome regions, with >80% of Single Nucleotide Polymorphism (SNPs) found in three regions. In contrast, the rest of the genome was intermingled between the two forms, with overall rates of HR estimated at 1000× mutation rates. One of their three questions for future investigation was "what are the barriers to gene flow between sympatric ecological populations?", which we investigate here. In order to focus on the effects of shared versus separate gene pools in our model, we assume that cells do not move between the two habitat types (i.e., particle types), whereas DNA is either able to move between the niches (global gene pool version) or not (local gene pools version) in separate versions. Our model belongs to a class of models considering the joint effects of selection, recombination and gene flow (often called migration but referring solely to transfer of DNA here, rather than movement of individuals) between two habitats (Lenormand, 2002;Gavrilets, 2004; Supporting information Appendix S1).

| Reproduction and selection
We implement a discrete time Wright-Fisher model of drift in a haploid population of size N, in which each iteration is equivalent to a generation (e.g., as applied to bacteria by Fraser, Hanage, & Spratt, 2007). Simulations follow a single ecological trait determined by L loci. A locus can be thought of as a single gene, or a set of tightly linked genes such as an integron or a plasmid (Yeaman and Whitlock 2011). Due to computational limitations, N is 2,000 individuals, 1,000 in each subpopulation, and L is set to ten loci.

| Mutation
Mutation occurs at a rate m = 0.005, specified per extant gene copy ( Figure 1). We implemented two versions. In the model with fixed effect sizes, each locus has two alleles with value +a or −a, and symmetrical mutation between them. For example, with 10 loci, all a il = +0.1 or −0.1. In the model with evolving effect sizes, as in Yeaman and Whitlock (2011), we assumed that mutation changes the effect size of a locus according to a normal distribution centered on the current allele's effect size with variance 2 m .

| Homologous recombination
Homologous recombination occurs at rate r and is assumed to be a process of gene conversion, whereby a gene copy (chosen in the same manner as the gene copy to be mutated) is replaced by an allele from the same locus. The replacing allele is chosen at random from a pool containing all extant alleles at that locus (i.e., it is thus possible for a gene copy to recombine with itself). In different versions, we assumed either that gene pools extend globally across the whole population or separately within each subpopulation, that is, from a local, niche-specific gene pool.

| Gene loss and gain
Gene loss (rate λ) is implemented in the same way as mutation ( Figure 1), but results in the matrix position a il losing its numerical value: that is, when a locus is lost, its effect size is assumed to be zero. Such loci are exempt from mutation or homologous recombination: we assume that HR requires tracts of homology with incoming DNA, but these are absent after gene loss (Retchless & Lawrence, 2007). In gene gain models, such "empty" locus slots can be re-occupied by a copy of the equivalent locus from another individual, at a rate γ, which is specified per empty locus. If no other individual carries a copy of that particular locus, the locus has gone extinct.
Note that we model gene gain as an independent process from HR.
In reality, gene gain might occur through HR at flanking regions to the lost gene. The rate of this event would be lower, however, than the rate of HR with the gene present, because there is less homologous sequence present for incoming DNA to match to (Retchless & Lawrence, 2007). Even if the same genetic mechanisms are involved, it is, therefore, appropriate to model gene gain as a different process from HR into a genotype already possessing the gene.

| Parameter values
We ran models with three values of recombination rate. r = 0.05 is 10× the mutation rate, which is lower than estimates in the Vibrio study of Shapiro et al. (2012), but above the level expected to generate quasi-sexual dynamics in neutral models (Fraser et al. 2007), and within the top 21% of estimates of r/m across a range of bacteria including other Vibrio (Vos & Didelot, 2009 is well known to depend on the strength of selection relative to recombination and gene flow (Bürger, 2014;Gavrilets, 2004): we decided to hold selection fixed and vary HR and patterns of gene flow to explore different scenarios. We used a deterministic model to choose parameter values that yielded phenotypic divergence in the clonal model but noticeably reduced phenotypic divergence in a model with recombination and a global gene pool (Supporting information Appendix S1: Figure S1). These conditions allow us to identify which extra features of our simulations are able to promote further divergence. With 2 w = 2 in Equation 1, the relative fitness of a genotype with one maladapted allele from the other subpopulation is 0.99 with 10 loci and equal effect sizes (e.g., a phenotype of 0.9-0.1 = 0.8 in the niche with optimum phenotype = 1), falling to 0.78 if the effect is shared equally between just two loci and to 0.37 if the effect is concentrated in one locus. Note that our models only consider symmetrical patterns of selection and gene flow between niches: asymmetry can lead to additional outcomes such as gene swamping (Lenormand, 2002) but our concern here is understanding the impact of gene flow in scenarios where divergence is otherwise able to occur.

| Simulations
Each run started with a fully homogenous population with an optimal phenotype for niche 1. After a burn-in of 1,000 generations of neutral reproduction and mutation (chosen to optimize computation times across all of the simulations), the ancestral population was split into two populations, and the simulation continued for 5,000 generations (flow diagram in Figure S2). Individuals cannot migrate between subpopulations and rates of mutation, homologous recombination, gene gain or gene loss did not differ between niches.  Table 1. Code development, profiling using profvis (Chang and Luraschi, 2016), and data analysis were performed in R 3.2.3 (R Core Team, 2015), while simulations ran on Imperial College London's high performance computing cluster in R 2.13.0 (R Core Team, 2011).

| Differentiation and summary statistics
We measured phenotypic differentiation as the magnitude of (mean phenotype population 1-mean phenotype population 2). Linear models were used to evaluate treatment effects. For simulations in which genome content was variable, the number of extant loci, and the number of extant gene copies per surviving locus were counted at the end of the simulation run.
In models of evolving effect sizes, we followed Yeaman and Whitlock (2011) to capture the degree of divergence among loci. As multiple alleles can exist within a subpopulation, only the most frequent or "leading" allele was considered. Between subpopulations, divergence at a given locus was calculated as the absolute difference in effect size of the leading alleles, d = |a 1 −a 2 | (where a 1 and a 2 are the leading alleles in each subpopulation). d was calculated for every locus, and then further summarized as the divergence of the locus with the largest divergence, d max , and the mean divergence, d mean , which were both recorded each generation. The greater the difference between d max and d mean , the more a single locus (namely that of largest effect) contributes to the divergence between both subpopulations. For the simulations involving variation in genome content, gene absences were treated as loci of effect size zero, and can be the "leading allele" when most individuals in a subpopulation have lost that locus.

| Fixed effect sizes and genome content
As a baseline for comparison with models of evolving effect sizes, we first ran a set of models assuming fixed effect sizes (i.e., all a il = +0.1 or −0.1) and no gene loss or gain. We first report results   These results confirm expectations from earlier models (Gavrilets, 2004;Friedman et al., 2013; and results for a deterministic two-locus model in the Supporting information Appendix S1: Figure S1) and show that our choice of parameter values provides suitable variation to explore effects of evolving genetic architecture in the next section: the extent of divergence varies both with respect to recombination rate and to global versus local gene pools.

| Fixed effects with gene loss and gain
The simulations were repeated with models including either just gene loss or gene loss plus gain in turn, still with fixed and equal ef- Appendix S1: Figure S3). In turn, gene loss protects local adaptation from the homogenizing effects of gene flow.
With the local gene pool model, the effects of gene loss are reduced and constitute a small decline in phenotypic divergence irrespective of the rate of recombination. In common with the clonal model, the only effect is the loss of some loci during initial divergence: recombination marginally reduces the rate of loss, but the effect is not significant (Figure 4).

| Phenotypic divergence with evolving effect sizes
The simulations were next repeated for the model of evolving effect sizes. To allow for the longer timescale to reach stationary conditions, these models ran for 500,000 generations (as in Yeaman With selection and a global gene pool, the optimal divergence of two phenotypic units is attained in most parameter combinations.

| Genetic divergence with evolving effect sizes: fixed genome content
To further understand these patterns, we looked at how gene effect sizes changed across loci and among treatments, first under As a locus can have positive as well as negative effects, the only constraint selection imposes is that the gene effects sum to a value close to the optimum. With a high mutation rate and many loci, a mildly deleterious mutation in one locus can be easily offset by a mutation elsewhere, resulting in a nearly neutral accumulation of mutations (see also Yeaman & Whitlock, 2011). Consequently, the locus of largest divergence increases steadily in effect throughout the simulation (Supporting information Figure S4, black line), and is compensated by other loci when its divergence exceeds the optimum value.
With recombination and a global gene pool, divergence is con-  Figure S5). This is because selection more efficiently removes maladapted genes arriving into each population: no intermediate phenotypes arise that are less strongly selected against but carry maladapted alleles (Supporting information Appendix S1). The benefit increases at higher r (Supporting information Figure S5 Figure S6). On average, it took each replicate simulation roughly 200,000 ± 40,000 generations to reach the equilibrium level, with a range of roughly 40,000-400,000 generations. This occurs because the populations become "stuck" at a local optimum of two loci, each with two alleles locally adapted to each niche, that is, locus A with allele effect sizes +0.5 and −0.5, and locus B with allele effect sizes +0.5 and −0.5. Any initial single mutation to progress from that genotype to the optimal genotype (i.e., a single A locus with allele effect sizes +1 and −1) leads to a loss of mean fitness: only a few combinations of 2 or 4 mutations lead to increased mean fitness (Supporting information Figure S7). In contrast, populations do not become stuck at the local optimum when r = 0.05 because (a) there is weaker selection against single mutants, (b) high recombination leads to faster production of beneficial combinations and (c) high recombination rate leads to stronger selection pressure for concentrating gene effects into a single locus (Supporting information Figure S7).
In contrast, recombination and local gene pools (Figure 6c) lead to a diffusion of divergence throughout the genome similar to that observed in strictly clonal populations. The differences between the local and global gene pool models mirror predictions in sexual eukaryotes: allopatric speciation results in even divergence across the genome whereas speciation with gene flow leads to divergence concentrated in genomic islands under selection (Ellegren et al. 2012).
Note that, although there is no gene flow between niches in the local gene pool model, divergence across the genome is lower at the higher recombination rate. Runs with r = 0.005 display divergence patterns very similar to the clonal model: the locus with greatest divergence steadily increases in divergence over time with compensation at other loci (Supporting information Figure S4b, dark gray line).
In contrast, with r = 0.05, divergence at the locus with the greatest divergence slows down over time. Once the optimum phenotypic divergence is reached, r = 0.05 selects against new mutations with large effects and subsequent compensating mutations because these are incompatible with the prevailing genetic background, which they are tested against more frequently than when r = 0.005. This demonstrates additional effects of HR within populations on the evolution of genetic architecture.

| Genetic divergence with evolving effect sizes: gene loss and gain
With gene loss and gene gain, by the end of the simulation on average, 88.6% of genes are lost from both populations, irrespective of treatment. As above, these results are to be expected under our model of additive gene effects and high mutation rate, whereby one locus can evolve to replace the phenotypic effect of all other loci. As before, gene gain had little effect.
Upon closer inspection, 71% of simulation runs involving loss ended with only two loci remaining, and only one run ended with as many as four surviving loci. Each surviving locus was found on average in 1,019 cells. This implies that for the majority of simulations, only a single locus survives in each subpopulation (of size 1,000 individuals), and that these are usually different loci. Among the simulations involving local gene pools, however, there were three runs that had only a single locus remaining (i.e., the same locus in both populations).
Again, these results can be interpreted using a simpler model of two biallelic loci (Supporting information Appendix S1: Figure S5).
Comparing scenarios with global gene flow and with and without gene loss, the optimal divergence is obtained when effects are partitioned between loci, that is, cells in population 1 have a locus with alleles a and A with effect sizes +1 and 0, whereas population 2 has a different locus with alleles b and B with effect sizes −1 and 0. No transfer of maladapted alleles can occur because HR cannot occur between different loci. At low recombination rates, the benefit over the single locus, specialist allele solution that evolves in the absence of gene loss is minimal, but it becomes greater at higher recombination rates: hence the lower phenotypic divergence without gene loss than with gene loss when r = 0.05. Testing these qualitative predictions will require more systematic studies of bacterial divergence in nature and in the lab.
Unfortunately, there has been surprisingly little research into the genetic architecture of bacterial traits. Although there is some evidence for the sort of additive gene effects we modeled (Abe & Benedetti, 2016;Hunter & Keener, 2014), many traits display more complex interactions (Brbić et al., 2016). However, we believe our qualitative predictions are likely to hold at the sequence level as well.
For example, in the motivating study for our model, Shapiro et al. (2012) characterized the genomes of two recently diverged populations of the marine bacterium Vibrio cyclitrophicus. They found evidence of high homologous recombination (enough to obscure any signal of clonality), and a strong trend for recent recombination to happen mainly within, rather than between populations. The single nucleotide polymorphisms (SNPs) that supported ecological differentiation between the two populations were spread among 11 loci, of which three accounted for over 80% of SNPs. Although this distribution is consistent with the qualitative predictions of both our model and that of Friedman et al. (2013), it is also possible they reflect other mechanisms influencing the evolvability of these loci.
Selection experiments of divergence in laboratory bacteria, with and without gene flow, would also be a useful approach for testing the evolution of different processes described here. Although we have emphasized recombination and gene flow here, evaluation should consider the range of recombination rates observed across bacteria-many of which are indeed effectively clonal and ecotype models of divergence should apply.
Our models also show that gene loss can further protect locally adapted genes from an influx of maladapted genes by leading to private loci adapted to each environment. Gene loss is a well-known phenomenon in bacteria and is often interpreted in terms of specialization: generalists carry many genes to enable metabolic flexibility, but specialists lose costly genes that are unnecessary in their specific environment (Koskiniemi, Sun, Berg, & Andersson, 2012).
Especially in bacteria with little or no horizontal gene transfer, phenotypic diversity between strains can originate through differential loss of genes, as found in the highly clonal Mycobacterium tuberculosis species complex (Bolotin & Hershberg, 2015). Our model highlights a new potential force on gene loss in bacteria, namely that differential gene loss could be intensified when two differentiating specialists maintain the potential to exchange DNA.
The comparison of patterns of divergence between E. coli and Salmonella enterica across genome regions by Retchless and Lawrence (2007)  region associated with production of the protective envelope of the heterocyst, and is associated with the deletion of two genes. While gene loss might be favoured through direct effects on the heterocyst phenotype, reduction of gene flow in this genome region in order to maintain ecologically adapted genotypes could also play a role, as predicted by our model.
Our gene-loss mechanism requires that multiple, exchangeable loci contribute to variation in a trait under differential selection between populations. For example, changes in the utilization of a particular carbon source could arise through multiple changes in underlying enzymes involved in its metabolism and the regulatory pathways that control them. If divergent selection acted instead on a set of multiple traits each coded by separate loci that are all essential for cellular function, then clearly gene loss could not aid divergence.
A more realistic version of our model could incorporate loci targeted by different components of selection that are either shared or divergent between populations. It remains to be seen if the predicted interaction between high rates of HR plus gene flow and gene loss will hold as more empirical studies on the genomics of bacterial divergence in sympatry become available. Cases where the selective advantage of having a given locus is outweighed by the disadvantage of high gene flow at that locus may turn out to be rare in nature.
However, in the lab, bacteria can quickly evolve compensatory mutations and restore growth, even after the loss of a core glycolysis gene (Charusanti et al., 2010).
In contrast to gene loss, gene gain had little impact in our simulations. This is consistent with the idea that gene gain is only rarely selectively advantageous and largely plays a role in adaptation to a novel niche. For example, Niehus et al. (2015) showed that horizontal transmission of genes dominates a population (or community) only insofar as there is an influx of nonadapted strains lacking an ecologically relevant gene. In the absence of migration, the selective advantage of individuals with the gene in question is large enough for transmission to be almost exclusively driven vertically. The limited effects of gene gain also demonstrate that the emergence of private genes in our gene-loss models was not an artifact of our assumption that genes cannot re-enter the population once lost.
Further mechanisms not considered here could be incorporated into future models. In interpreting our model, each "locus" could be either one or multiple genes. Thus, the concentration of gene effects we report can be seen either as a consequence of cumulative mutations in one gene, or the concentration of multiple previously dispersed genes into, for example, a single operon or integron (see also Yeaman & Whitlock, 2011). A future modeling effort may wish to distinguish between these two processes; we postulate that a process involving genome rearrangements would result in similar concentration of gene effects over time.
A final mechanism for reducing maladaptive gene flow between These predictions need to be tested empirically across bacteria with a wide range of recombination rates.

ACK N OWLED G M ENTS
Thanks to Reuben Nowell, Alberto Pascual and Tin-yu Hui for discussion, James Rosindell for advice on using the computing cluster, and Ricardo Azevedo and 3 reviewers for comments.

AUTH O R CO NTR I B UTI O N S
Both authors devised the study, programed the model, ran simulations, performed analyses and wrote up the manuscript.

DATA ACCE SS I B I LIT Y
All scripts for generating and analyzing simulation outputs: Dryad