Evolutionary Rescue Is Mediated by the History of Selection and Dispersal in Diversifying Metacommunities

Rapid evolution can sometimes prevent population extirpation in stressful environments, but the conditions leading to “evolutionary rescue” in metacommunities are unclear. Here we studied the eco-evolutionary response of microbial metacommunities adapting to selection by the antibiotic streptomycin. Our experiment tested how the history of antibiotic selection and contrasting modes of dispersal influenced diversification and subsequent evolutionary rescue in microbial metacommunities undergoing adaptive radiation. We first tracked the change in diversity and density of Pseudomonas fluorescens morphotypes selected on a gradient of antibiotic stress. We then examined the recovery of these metacommunities following abrupt application of a high concentration of streptomycin lethal to the ancestral organisms. We show that dispersal increases diversity within the stressed metacommunities, that exposure to stress alters diversification dynamics, and that community composition, dispersal, and past exposure to stress mediate the speed at which evolutionary rescue occurs, but not the final outcome of recovery in abundance and diversity. These findings extend recent experiments on evolutionary rescue to the case of metacommunities undergoing adaptive diversification, and should motivate new theory on this question. Our findings are also relevant to evolutionary conservation biology and research on antimicrobial resistance.


INTRODUCTION
Biodiversity is organized across multiple scales, from populations of individual species to communities of multiple species, that evolve and move across a range of spatial scales, from local ecosystems to entire continents. Situated within this hierarchy are metacommunities, defined as a set of local communities that are connected by dispersal (Leibold et al., 2004). Metacommunities face increasing pressure from human-induced environmental degradation, often leading to population decline, local extinctions, and biodiversity loss (Pereira et al., 2010;Haddad et al., 2015). Species occupying metacommunities can respond to these pressures by moving to other communities, by adapting in situ, or they can undergo extinction. Under some circumstances, populations can rapidly evolve resistance to stressors and persist in severely degraded environments (Bürger and Lynch, 1995;Hufbauer et al., 2015;Reid et al., 2016). This phenomenon is described by the theory of evolutionary rescue (Gomulkiewicz and Holt, 1995;Bell and Gonzalez, 2009;Bell, 2017), which explains how populations can recover from extreme environmental stress through rapid adaptation.
Evolutionary rescue occurs when stress-resistant individuals are selected within a perturbed population, allowing abundance to recover and a viable population to be maintained in conditions that would otherwise cause population extirpation in the absence of evolution. Resistant types responsible for evolutionary rescue may already be present in a population before rescue is needed (due to past in situ evolution or to the immigration of resistant types from connected habitats), or evolve de novo after the onset of extreme stress. Theory and laboratory experiments have described key drivers of evolutionary rescue for singlespecies populations. First, evolutionary rescue is more likely if populations were previously exposed to lower levels of the same stressor (leading to an increase in the abundance of stressresistant genotypes through selection). Second, large population sizes favor evolutionary rescue by reducing the risk of stochastic extinction (Gomulkiewicz and Holt, 1995;Bell and Gonzalez, 2009;Gienapp et al., 2013). Third, if populations are spatially connected to others, the dispersal of stress-resistant genotypes among local populations can enable rescue even in environments that were not previously contaminated by the stressor (Bell and Gonzalez, 2011;Carlson et al., 2014). Understanding the conditions that favor evolutionary rescue has clear implications for both biodiversity conservation and for the management of pests and pathogens evolving resistance to biocides (Alexander et al., 2014).
A few recent studies have expanded the scope of evolutionary rescue from populations to diverse assemblages of species and examined evolutionary rescue in entire communities confronted with severe stress (Fussmann and Gonzalez, 2013;Low-Décarie et al., 2015;Bell et al., 2019;Fugère et al., 2020). Community rescue occurs when the populations of multiple species recover rapidly following exposure to levels of stress that were lethal to the community in its ancestral form, allowing a community to recover both abundance and diversity in severely degraded environments (Low-Décarie et al., 2015). Some key factors promoting evolutionary rescue of populations were also found to favor community rescue (Low-Décarie et al., 2015;Bell et al., 2019;Fugère et al., 2020). First, a history of stress exposure increases the relative frequency of stress-resistant individuals in communities, which provides a correlated advantage at a higher dose of stress and thus facilitates rescue. Second, just as intraspecific genetic diversity can promote the evolutionary rescue of individual populations (Carlson et al., 2014), a greater diversity of species can also favor rescue, as communities holding a more diverse set of species are more likely to contain resistant types (Low-Décarie et al., 2015). Finally, when local communities are spatially connected, thus forming a metacommunity, dispersal was found to favor evolutionary rescue in local communities by moving resistant genotypes within a heterogeneous metacommunity (Low-Décarie et al., 2015). Nonetheless, only a few recent studies have tested the conditions that promote evolutionary rescue in communities of multiple species, and much remains to be uncovered.
Here, we expand on previous work testing evolutionary rescue in communities, and ask whether a history of stress, greater biodiversity, and the presence of dispersal favor evolutionary rescue in metacommunities across a heterogeneous landscape. In contrast to previous community rescue experiments that have used microbial assemblages with existing intra and interspecific variation, we examined evolutionary rescue in an experimental system in which all diversity and stress resistance is generated de novo through rapid evolution. We also manipulated not only the presence of dispersal prior to exposure to extreme levels of stress, but also its spatial structure, contrasting global vs. local dispersal in metacommunities. These two modes of dispersal were shown to have distinct effects on the likelihood of evolutionary rescue in metapopulations of the yeast, Saccharomyces cerevisiae (Bell and Gonzalez, 2011), but how the spatial structure of dispersal affects evolutionary rescue in metacommunities remains unknown. Global dispersal connects all communities to each other, mixing individuals from the whole metacommunity. This mode of dispersal brings in diversity upon which natural selection will act and potentially enable evolutionary rescue. At the same time, the resulting migrants could be maladapted to their new habitat, thus hampering adaptation and subsequent evolutionary rescue by migration load (Bell and Gonzalez, 2011;Schiffers et al., 2013;Carlson et al., 2014;Hao et al., 2015). On the other hand, local dispersal [also known as "stepping-stone" dispersal (Bell et al., 2019)] is a mode of dispersal where communities are only connected by directional migration up a gradient of environmental stress. Local dispersal is expected to favor evolutionary rescue, because it will more likely bring betteradapted individuals as they move up the stress gradient (Bell and Gonzalez, 2011). Over longer time-scales in communities undergoing adaptive diversification and speciation, dispersal may also favor the generation of diversity and productivity (Venail et al., 2008), which would also promote the likelihood of rescue following environmental degradation.
To address how stress and dispersal interact to modulate the adaptation and diversification of metacommunities across a heterogeneous landscape, and subsequent evolutionary rescue, we used the plant symbiont Pseudomonas fluorescens SBW25 as a model system exposed to the antibiotic streptomycin. We build on previous studies with this organism (Rainey and Travisano, 1998;Kassen et al., 2000;Massin and Gonzalez, 2006;Perron et al., 2006;Venail et al., 2008Venail et al., , 2010Ramsayer et al., 2013). This bacterium shows rapid and repeatable in vitro diversification when grown in a heterogeneous environment. As they grow, individuals of this aerobic bacterium compete for oxygen, thus creating a vertical gradient of oxygen in the liquid medium. The resulting environment favors mutants able to colonize the air-liquid interface (Rainey and Travisano, 1998) via the formation of a biofilm composed of cellulosebased polymer (McDonald et al., 2009;Lind et al., 2018). Diversification can be recorded at the phenotypic level by growing these bacteria on solid media where they display striking differences in colony morphology that relates to their niche preference. These morphotypes are easy to detect and are heritable (Rainey and Travisano, 1998). In our experiment, the ancestral cells were derived from a single isomorphic colony ("smooth" opaque morph) which then grew asexually. This means that mutations are only transferred to the next generation by single cell division. Consequently, adaptive mutations that arise in one morphotype are independent from the evolutionary pathway of other morphotypes. Quantifying the emerging morphological diversity allows us to track this evolution occurring in vitro. The species concept does not readily apply to bacteria (Rosselló-Mora and Amann, 2001;Riley and Lizotte-Waniewski, 2009), but the subsequent rapid diversification, niche specialization, and growth by asexual reproduction in P. fluorescens allows us to consider different morphotypes analogous to different species, and to consider our diversified bacterial assemblages as a model for communities of multiple species. Previous studies have documented the capacity of P. fluorescens to rapidly evolve resistance to the antibiotic streptomycin (Ramsayer et al., 2013), a versatile and widely used antibiotic (World Health Organization, 2015). This work, combined with the tendency for P. fluorescens to adaptively radiate (MacLean and Bell, 2002;Barrett and Bell, 2006), makes it an excellent system to study the factors promoting adaptation to stressors and evolutionary rescue in the context of rapidly diversifying communities.
Following previous experiments (Bell and Gonzalez, 2011;Low-Décarie et al., 2015;Bell et al., 2019;Fugère et al., 2020), we conducted an experiment which proceeded in two phases. In the first phase, communities evolved and diversified across a gradient of streptomycin. We created replicated four-patch metacommunities with one of three modes of dispersal: local dispersal, global dispersal, and a control with no dispersal. We recorded the dynamics of adaptation by quantifying the growth and morphological diversification of each community across the gradient of stress. In the second phase, we transferred each community regardless of its history of stress to a severe dose of streptomycin which was established to be lethal to the majority of ancestral organisms, following the method employed by previous community rescue experiments (Low-Décarie et al., 2015). In this second phase, dispersal was ceased such that any recovery of abundance and diversity could be attributed to local eco-evolutionary processes -and not demographic rescue due to dispersal. We quantified the trajectory and outcome of evolutionary rescue in Phase 2 of the experiment and linked these responses to three potential drivers of rescue manipulated in Phase 1: the history of exposure to sublethal doses of streptomycin, the mode of dispersal within the metacommunity, and the morphotype diversity of the community (i.e., the outcome of diversification occurring during Phase 1). Based on previous experiments (Bell and Gonzalez, 2011;Low-Décarie et al., 2015;Fugère et al., 2020) we expected that: (1) a history of streptomycin exposure in Phase 1 would facilitate evolutionary rescue in Phase 2; (2) the presence of dispersal in heterogeneous metacommunities would increase local morphotype diversity and would spread resistant genotypes during Phase 1, both of which would facilitate adaptation to, and rescue from, severe antibiotic stress-especially in communities naive to the stressor; and (3) that local dispersal would have a greater influence on the likelihood of evolutionary rescue than global dispersal.

Bacterial Cultures
We used the ancestral strain of Pseudomonas fluorescens SBW25 (Rainey and Bailey, 1996) cultured in basic growth medium was King's B (KB) medium (20.00 g/L Proteose Peptone (Difco no.3), 15 mL/L glycerol, 1.50 g/L K 2 HPO 4 , 1.50 g/L MgSO 4 , distilled water). Populations in stressful treatments were supplemented with streptomycin (MilliporeSigma: Montreal, Canada). We initiated bacterial cultures from a single isolate clone of P. fluorescens SBW25 in a 125-mL glass vial supplied with 50 mL of King's B medium, grown for 24 h at 28 • C and shaken at 150 rpm. One percentage of this culture was transferred to 96-well plates supplied with 200 µL KB and grown for a further 24 h at 28 • C, to initiate separate experimental populations. These were not shaken to allow diversification into morphologically diverse communities of P. fluorescens morphotypes over the course of the experiment (Rainey and Travisano, 1998). We thus refer to the "ancestral cells" at the onset of the experiment, and then to "communities" once microwells contained a diverse assemblage of morphotypes.
Abundance was measured spectrophotometrically as the optical density at 590 nm (OD 590 ), using a microplate reader (BioTek: Winooski, USA). Correspondence between absorbance readings and cell density was verified by growing bacterial cultures with known OD 590 values on at least 3 replicate KB-Agar plates and by counting the number of colonyforming units. The relationship between cell density and OD 590 (R² = 0.6) was: y (cells/mL) = − 3 × 10 8 + 5 × 10 9 × OD 590 (Supplementary Figure 1). Population size of the ancestral population was recorded after 24 h of growth in benign media, where it reached an abundance of ∼12 × 10 9 cells/mL (cell density recorded on agar, after 10 6 -fold dilution). This population served to initiate the experimental metapopulations with a standardized initial abundance of 10 5 cells/mL in 96-well plates supplied with 200 µL of KB medium (see below). Then, throughout the experiment, we tracked community abundance by measuring the OD 590 of all plates every 24 h.
To determine the susceptibility of the ancestral bacteria to streptomycin, we inoculated the ancestral bacteria of P. fluorescens SBW25 in densities of 10 5 cells/mL in wells with 200 µL of KB medium supplemented with 10 different concentrations of streptomycin, ranging from 0 to 500 µg/mL. Eight replicate populations were treated with one of 10 streptomycin concentrations for 24 h (incubated at 28 • C), during which OD 590 was read every 30 min (Supplementary Figure 2). This served to identify the dose of streptomycin that was lethal to the majority of ancestral cells, the precondition required for evolutionary rescue. Note that, just as in previous experiments that tested population and community rescue, this dose is not required to be lethal to all individuals (Bell and Gonzalez, 2009;Low-Décarie et al., 2015).

Experimental Design
We randomly assigned metacommunities to streptomycin or control treatments crossed with three modes of dispersal (global, local, and none). Each metacommunity consisted of four communities. Streptomycin-treated metacommunities were composed of four communities exposed to concentrations of 0, 100, 200, and 400 µg/mL, replicated 4 times ( Figure 1A). These concentrations partially inhibited growth and represented a selection pressure on P. fluorescens during the first part of the experiment, which is expected to generate a correlated genetic response conferring some degree of tolerance to severe stress before it was actually experienced in Phase 2 of the experiment. Control metacommunities had the same dispersal treatments but were unexposed to streptomycin. This factorial design resulted in 96 local communities arrayed into 24 metacommunities distributed evenly on 4 separate 96-well plates.
During Phase 1 of the experiment, 1% of each culture was transferred to a new plate with fresh medium every 24 h to maintain growth. Dispersal treatments within metacommunities occurred simultaneously with the transfers. In the two dispersal treatments, 2 µL of grown culture was moved from each well to a dispersal pool, to match the rate of 1% dispersal used in previous studies with this bacterial strain (Venail et al., 2008). In the local dispersal treatments, this pool contained a contribution from the well with the next lower streptomycin concentration on the old plate: 2 µL of this dispersal pool was transferred to the next higher level of streptomycin, moving migrants up the stress gradient in the case of stressed communities. In the global dispersal treatments, the pool contained equal contributions from all wells of the metacommunity, and 2 µL of this was distributed to all wells of the metacommunity. In the no dispersal treatment, each well was inoculated exclusively with 1% from the corresponding well on the old plate, so that no cross-well transfer occurred.
Phase 2 of the experiment started 24 h after the seventh transfer. We diluted the grown communities to 10 5 cells/mL in KB (matching densities of 10 5 cells/mL for which we had assayed the tolerance of P. fluorescens SBW25 to streptomycin) supplemented with 500 µg/mL streptomycin, a concentration that was lethal to the great majority of ancestral cells (Supplementary Figure 2). Bacteria were incubated at 28 • C for 4 days, with no transfer or dispersal events. After the transfer, the abundance of each community was recorded by absorbance readings (OD 590 ) after 24 h (day 8), 30 h, 48 h (day 9), 72 h (day 10), and 120 h (day 12). Morphological diversity was scored after plating cells on KB-Agar on two occasions: on day 8 (24 h after the beginning of Phase 2) and on day 12 (at the end of Phase 2).

Measuring Morphological Diversity
Counts of morphotypes were scored every day during Phase 1 and on two instances during Phase 2, by plating on KB-Agar after dilution in KB. Each community was sampled and grown on two replicate Petri dishes. The morphotypes of all colonies were scored visually after 3 days of growth at 28 • C; this corresponded to 50-500 colonies per replicate community. The two values obtained from replicate Petri dishes were then averaged to give the composition of each community. In some cases, when colonies failed to grow or when colonies grew very quickly and fused into a continuous biofilm, morphotype composition could not be estimated reliably. Thus, some replicates had to be discarded in analyses of composition and diversity (see below).
In keeping with previous work (Rainey and Travisano, 1998), we identified three morphs: smooth morph (the planktonic, ancestral morph), wrinkly spreader (which colonizes the airliquid interface by forming a biofilm within 1 to 2 days) and fuzzy spreader (which colonizes the bottom of the wells after 4 days). We further divided the smooth morphotype into five subclasses -all had a smooth appearance but differed in opacity and pattern: "smooth morph translucent, " "smooth morph opaque" ( Figure 1B); "shiny smooth morph;" "eclipse smooth morph;" and "radial smooth morph." We also divided the wrinkly spreader morphotype into three subclasses -all contained "wrinkles" (distinctive asperities on the surface of the colony), but differed in the extent to which these wrinkles covered them: "wrinkly spreader 1" colonies were completely covered in wrinkles, "wrinkly spreader 2" colonies had a wrinkly center and smooth edge, and "wrinkly spreader 3" colonies had a smooth center and a wrinkly edge. This diversity observed within the wrinkly spreader morph is consistent with previous work with P. fluorescens (Hodgson et al., 2002;Massin and Gonzalez, 2006;Bantinaki et al., 2007). Finally, fuzzy spreader colonies displayed a distinctive blurry edge. Diversity counts thus included a total pool of nine morphological types ( Figure 1B). Diversity in a local community was computed as the exponent of the Shannon index: where S is the number of morphotypes and p i is the relative abundance of morphotype i in the community. The exponent of Shannon is known as the effective number of species (Jost, 2006) -or in this case the effective number of morphotypes (henceforth diversity).

Statistical Analyses
All statistical analyses were conducted in R version 3.4.3 (R Core Team, 2017). We used linear models (ANOVAs) to test the effect of dispersal mode (three levels), streptomycin concentration (four levels, in the case of exposed metacommunities only), and their two-way interaction on the abundance (OD 590 ) and diversity of communities at key time points of the experiment: at the end of Phase 1 (day 7), the beginning of Phase 2 (day 8) and the end of Phase 2 (day 12). We evaluated the response of control metacommunities separately from exposed metacommunities. We also included experimental replicates (the 4 micro-well plates) as a blocking factor in all models. The number of replicates for each day of the experiment is referenced in the Supplementary Material (Supplementary Table 1 for  abundances and Supplementary Table 2 for morphological diversity replicates). When main effects were statisticallysignificant, we conducted post hoc pairwise comparisons of factor levels using Tukey's "Honest Significant Difference" method (function "TukeyHSD" in R). Response variables were log-transformed when it improved the normality and homogeneity of model residuals.
To visualize differences in morphotype composition among treatments, we performed non-metric multidimensional scaling (NMDS) analyses separately at the end of Phase 1 and Phase 2, combining control and exposed metacommunities (R-package "vegan" version 2.5-6, function "metaMDS"). These NMDS ordinations and all other multivariate analyses Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 1 | Experimental design and morphotypes. (A) Experimental design for Phase 1, where each circle represents a single microwell, blue shades and numbers represent streptomycin concentrations (in µg/mL), and arrows represent the mode and direction of dispersal between communities. This design was replicated 4 times on separate multi-well plates. (B) Photographs of typical bacterial colonies grown on agar. We identified nine morphotypes that emerged during diversification over the course of the experiment: smooth morph opaque (SMO), Eclipse smooth morph, Radial smooth morph, smooth morph translucent (SMT), Shiny smooth morph, fuzzy spreader (FS), wrinkly spreader (WS) submorphs 1, 2, and 3. Colonies were photographed under a 10-fold magnifying lens (Leica), illuminated from below (except for WS2). The size of colonies at this stage averaged 5 mm in diameter.
used the Bray-Curtis dissimilarity index and proportional (relative abundance) data when calculating distance matrices. We tested whether morphotype composition varied among treatments using permutational multivariate analysis of variance (pMANOVA) implemented with the function "Adonis" in vegan. Separate pMANOVA models were fitted for control and streptomycin-treated metacommunities, using as a grouping factor "dispersal" (for both exposed and unexposed metacommunities) or "streptomycin concentration" (for exposed metacommunities only). We also conducted two analyses to reveal treatment effects on beta diversity, the variance in morphotype composition among local communities. We first used permutation-based tests of multivariate homogeneity of groups dispersions (the "PERMDISP" procedure implemented in function "betadisper" in vegan) to compare beta diversity at the "treatment" scale (e.g., variation in composition among all communities grown in concentrations of 100 vs. 200 µg/mL of streptomycin). These tests used "plate" and "dispersal" as grouping factors (for control metacommunities), or "plate, " "dispersal, " "streptomycin concentration, " and their two-way interaction as factors (for exposed metacommunities). We then computed a metacommunity-scale measure of beta diversity, calculating, for each metacommunity of four microwells, the mean multivariate distance of the four local communities to their group centroid. Metacommunity-scale beta diversity was analyzed with ANOVA using "dispersal" (local vs. global vs. isolated), "metacommunity type" (control vs. streptomycintreated), their two-way interaction, and "plate" as factors. This analysis only included the subset of metacommunities in which morphotype composition could be reliably estimated for all four local communities.
Finally, to test the hypothesis that treatment effects on community composition and alpha diversity affected the trajectory and outcome of evolutionary rescue, we used independent linear regressions to link alpha diversity and community composition on day 7 (end of Phase 1) with OD 590 on day 8 and day 12. OD 590 on day 8 indicates the initial potential of communities to grow at a dose of stress lethal to the ancestral cells, while OD 590 on day 12 indicates final community abundance at this lethal dose of stress, i.e., the outcome of evolutionary rescue of multiple morphs. The measure of community composition used in this analysis were scores on the first axis of the NMDS ordination described above.

Ancestral Organisms
The ancestral bacterial culture was isomorphic and composed of Smooth Morph Opaque (SMO) exclusively. In these bacteria, growth was impeded by streptomycin concentrations as low as 50 µg/mL (ANOVA: F = 2,680, p < 0.0001), and no growth was recorded for concentrations of 400 µg/mL or more after 24 h (Supplementary Figure 2).

Phase 1
Growth was initially negatively affected by streptomycin, but communities displayed rapid adaptation to streptomycin during Phase 1 (Figure 2). During the first 24 h (∼30 generations Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 2 | Growth dynamics throughout the experiment. Each point represents the OD 590 of a single community, and lines represent the mean OD 590 across replicate communities of a given treatment combination. Shades of blue indicate Phase 1 streptomycin exposure while line types represent the type of metacommunity (control vs. streptomycin-exposed). The dashed red line represents the 95th quantile of abundance of control communities 24 h after the transfer of all communities to lethal concentrations of streptomycin. The duration of Phase 2 is shown by the pink polygon.
in benign conditions), cultures of P. fluorescens reacted to streptomycin in the same way as the ancestral cells: growth was significantly lower in cultures with higher concentrations of streptomycin in the environment (ANOVA: F = 238.6, p < 0.0001). However, by the end of Phase 1, while streptomycin still negatively affected growth at the highest concentrations (Tukey HSD between 0 vs. 200 or 400 µg/mL: p < 0.0001), communities that had evolved in concentrations of 100 µg/mL grew to abundance levels that were not significantly different than in benign environments (Tukey HSD: p = 0.91). By the end of Phase 1, OD 590 of communities exposed to streptomycin at 200 and 400 µg/mL were increased by 2 or 3-fold compared with the beginning of Phase 1. In contrast to streptomycin exposure, the dispersal treatment did not significantly affect Phase 1 abundance in neither control metacommunities (ANOVA: F = 0.812, p = 0.451), nor streptomycin-exposed metacommunities (F = 2.386, p = 0.108). The interaction between dispersal mode and streptomycin concentration was not significant in streptomycinexposed metacommunities (F = 1.151, p = 0.356).
Streptomycin exposure also influenced the dynamics of diversification. Starting from isomorphic, clonal cells, bacteria diversified into nine different morphotypes during the experiment (Figures 1B, 3A). In benign conditions, communities displayed a repeatable diversification pattern where SMO first dominated the community, coexisting with shiny smooth morph and SMT in lower abundances, followed by a rise within the first 2 days in the abundance of wrinkly spreader 1 (WS1), which replaced SMO as the dominant morph throughout the rest of Phase 1. As a result, the effective number of morphotypes increased in the first few days of Phase 1, reaching a maximum of 3.6 at the end of the 2nd day (averaged across control communities, s.d. = 0.72), and subsequently decreased as WS1 then dominated the community (Figure 3B). These repeatable diversification dynamics were altered in FIGURE 3 | Diversification dynamics of metacommunities throughout the experiment. Time series of morphotype relative abundance (A) and alpha diversity measured as Shannon exponent (B) in the different treatments. Plots are ranked from left to right according to streptomycin exposure during Phase 1 (values indicated at the top of the figure in µg/mL). "Controls" refer to control metacommunities in which none of the wells received streptomycin. Phase 2 is indicated by a red shading after the 7th transfer. In (A), we show fitted values of "LOESS" smoothing functions interpolating relative abundance over time. In (B), lines represent the average among replicate communities within each dispersal treatment, indicated by colors. Error bars = 1 standard deviation. communities exposed to streptomycin. In environments supplemented with streptomycin in concentrations of 100 and 200 µg/mL, communities displayed a consistently different diversification pattern where SMT rapidly became dominant, replacing SMO as the dominant morph within the first 2 days, and subsequently co-occurring with SMO and WS1 (and WS2 in highly stressful environments) throughout the rest of Phase 1. Communities exposed to streptomycin in concentrations of 400 µg/mL did not display such a consistent trend in morphological diversification, in part because the dispersal treatment modulated the increase in the relative abundance of SMT ( Figure 3A): SMT was abundant in concentrations of 400 µg/mL only in spatially connected communities.
Dispersal clearly influenced diversification dynamics and the generation of alpha diversity during Phase 1 (Figure 3). By the end of Phase 1, morphotype diversity was significantly higher in communities connected through dispersal, in both control metacommunities (ANOVA: F = 3.82, p = 0.0337) and exposed metacommunities (ANOVA: F = 5.24, p = 0.0148). However, the spatial structure of dispersal did not influence final Phase 1 diversity, as communities linked by global and local dispersal had comparable diversity, in both control metacommunities (Tukey HSD: p = 0.97) and exposed metacommunities (Tukey HSD: p = 0.73). In exposed metacommunities, streptomycin exposure had a weak but non-significant effect on diversity (ANOVA: F = 2.78, p = 0.068), while the interaction between streptomycin and dispersal was not significant (F = 0.75, p = 0.62).
The dispersal and streptomycin treatments influenced both the mean composition of communities and the variance in composition among local communities (Figure 4). In control metacommunities, dispersal did not have a statistically significant effect on morphotype composition (pMANOVA: p = 0.06), but had a very strong effect on heterogeneity among local communities receiving the same dispersal treatment ( Figure 4A; PERMDISP: p < 0.0001). Local dispersal reduced heterogeneity more than global dispersal. In streptomycinexposed metacommunities, the antibiotic gradient had a much stronger influence on composition than did dispersal. Indeed, while streptomycin exposure influenced the mean composition of communities ( Figure 4B; pMANOVA: p = 0.002), the main effect of dispersal on composition and the two-way interaction between dispersal and streptomycin concentration were not statisticallysignificant (pMANOVA: p = 0.26 and 0.66, respectively). Nonetheless, as in control metacommunities, dispersal reduced heterogeneity among local communities (PERMDISP: p = 0.03). Streptomycin exposure had the opposite effect of increasing heterogeneity among communities at higher doses ( Figure 4B; PERMDISP: p = 0.0003). At the metacommunityscale, variance in community composition among the four local communities forming a metacommunity was highly reduced by dispersal (Figure 4C; ANOVA, main effect of dispersal: F = 24.95, p = 0.001). This effect however disappeared in exposed metacommunities (Figure 4C; ANOVA, main effect of metacommunity-scale streptomycin treatment: F = 8.96, p = 0.02; two-way interaction effect of streptomycin and dispersal: F = 5.87, p = 0.04).
In summary, dispersal and streptomycin exposure over Phase 1 jointly affected the diversification of evolving Pseudomonas metacommunities, resulting in a gradient of morphotype diversity, composition, and history of stress across treatments, all of which were hypothesized to influence subsequent recovery of abundance and diversity in Phase 2.

Phase 2
In Phase 2, all communities were transferred to a concentration of streptomycin lethal to the ancestral cells (500 µg/mL). In the 24 h following exposure to 500 µg/mL of streptomycin, communities that had never been exposed to streptomycin in Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 4 | of NMDS ordination of community composition, indicating the location of local communities (points) and morphotypes (words) in multivariate space. For visualization purposes, results are shown separately for control metacommunities (A) and streptomycin-exposed metacommunities (B). Symbols denote dispersal treatments while colors indicate streptomycin exposure. Convex hulls regroup communities having received the same dispersal (A) or streptomycin (B) treatment. Morph acronyms are defined in Figure 1. (C) Metacommunity-scale beta diversity (mean distance to group centroid of 4 local communities comprising a metacommunity) as a function of metacommunity type and dispersal treatment. Symbols represent metacommunities while horizontal lines correspond to treatment means. S-, control metacommunities. S+, streptomycin-exposed metacommunities. I, Isolated; G, global dispersal; L, local dispersal.
Phase 1 underwent a 5-fold decline in abundance compared with the end of Phase 1 (Figures 2, 5A), confirming that this dose of streptomycin is highly stressful to the majority of cells even after in vitro diversification. In contrast, local communities that had been exposed to at least 100 µg/mL of streptomycin during Phase 1 recovered an abundance two to three times higher than naïve communities in the first 24 h of Phase 2 ( Figure 5A; ANOVA, main effect of streptomycin exposure: F = 14.8, p < 0.0001). However, dispersal altered this general trend (Figure 5A; ANOVA, interaction effect of streptomycin exposure and dispersal: F = 2.97, p = 0.02), as the beneficial effect of exposure to a high dose of streptomycin in Phase 1 was diminished in isolated communities. The spatial structure of dispersal did not influence these responses (Tukey HSD between global vs. local dispersal: p > 0.1 at all Phase 1 streptomycin concentrations). Contrary to our hypothesis, diversity levels achieved at the end of Phase 1 did not have a significant effect on the abundance recovered by communities at the onset of Phase 2 ( Figure 5B; linear regression: R² = 0.002, p = 0.71). Community composition at the end of Phase 1 did however have a significant effect on the abundance recovered in the first 24 h of Phase 2 ( Figure 5C; linear regression: R² = 0.25, p = 0.0028). More specifically, communities where SMT occurred in higher abundances at the end of Phase 1 recovered a significantly higher abundance following the first 24 h of Phase 2 (p < 0.0001). This morph eventually came to dominate most communities throughout Phase 2 (Figure 3A), while SMO, WS1 and WS2 occurred at lower abundances, and other rare morphotypes that existed in Phase 1 were lost entirely. For example, "radial" smooth morph was only detected in one of 143 communities throughout Phase 2, and "fuzzy spreader" in only five communities. Conversely, one morph ("wrinkly spreader 3") that had never been detected in Phase 1 was detected in 17 of 143 communities throughout Phase 2. However, in contrast to Phase 1 where differences in diversity were noted across treatments, in the first 24 h of Phase 2, communities displayed similar levels of morphological diversity with no significant effect of the treatment they had received during Phase 1 (ANOVA: F = 1.05, p = 0.4215). Importantly, all communities lost a large fraction of their cells at the onset of Phase 2 regardless of which morphs they contained (even communities dominated by SMT), suggesting that no morph was fully resistant to this dose of antibiotic and that within-morph evolutionary rescue was required for persistence in Phase 2.
Eventually, by the end of Phase 2, all communities recovered similar levels of abundance and diversity, with no significant effect of their history of stress and dispersal on final abundance and diversity (Figures 2, 3B; ANOVA for OD 590 : F = 1.07, p = 0.3997; for diversity: F = 0.91, p = 0.5544). Effects of Phase 1 treatments on community composition had also vanished by the end of Phase 2 (pMANOVA: p > 0.05 for all factors).
In sum, evolutionary rescue of multiple morphs occurred in all communities, allowing abundance and diversity to be maintained at a dose of antibiotic lethal to the majority of the ancestral, isomorphic cells. Thus, dispersal and the history of streptomycin exposure influenced the trajectory of evolutionary rescue (i.e., the initial decline in abundance and growth following lethal stress) but not the final outcome of rescue.

DISCUSSION
Since the original formulation of evolutionary rescue theory for isolated populations of single species (Gomulkiewicz and Holt, 1995), extending the theory to conditions closer to those at play in natural communities has been a major challenge in the study of evolutionary rescue (Fussmann and Gonzalez, 2013;Osmond and de Mazancourt, 2013). This challenge is now being addressed through careful experimentation (Low-Décarie et al., 2015;Fugère et al., 2020;Scheuerl et al., 2020). Here we studied evolutionary rescue in metacommunities undergoing adaptive diversification across a gradient of environmental stress. We have shown that the history of exposure to stress, and the presence and spatial structure of dispersal in metacommunities of the rapidly evolving bacterium P. fluorescens SBW25 can affect adaptation, diversification and evolutionary rescue following severe environmental degradation caused by antibiotic stress.
During the experiment, the isomorphic ancestral cells diversified in vitro into a total of nine morphotypes, forming phenotypically diverse communities that recovered both abundance and diversity following exposure to initially lethal levels of stress, demonstrating evolutionary rescue in these metacommunities. Mutations are the source of this phenotypic variation, allowing bacteria to colonize different ecological niches, in particular the surface of the liquid medium (Spiers et al., 2002;McDonald et al., 2009). Despite extensive genetic and phenotypic variation, however, the bacteria still belong to the same species, P. fluorescens SBW25. The species concept is blurry in bacteria (Rosselló-Mora and Amann, 2001;Riley and Lizotte-Waniewski, 2009) and challenging to quantify (Staley, 2006), and the case of rapid diversification in P. fluorescens is no exception. Nonetheless, asexual growth combined with niche specialization and morphological diversity supports the analogy between the different morphotypes and different species in a community, and ultimately allows us to interpret the recovery of diverse and abundant populations of morphotypes in terms of evolutionary rescue of multiple species forming a metacommunity. Our results show that: (1) dispersal and antibiotic stress jointly influenced diversification dynamics. At the end of this diversification FIGURE 5 | represent the median and quartiles of the OD 590 of replicate communities, colored by their mode of dispersal during Phase 1. In (B,C), each dot represents a single local community and the line represents the linear response of OD 590 on day 8 to diversity and community composition on day 7, respectively. Shapes represent the mode of dispersal, and shades of blue represents the concentrations of streptomycin communities were exposed to during Phase 1 (in µg/mL).
process, variation in community composition, but not diversity, subsequently influenced evolutionary rescue; (2) the history of exposure to stress was a strong predictor of the trajectory of evolutionary rescue, but not the outcome of rescue. Communities with a history of antibiotic exposure had greater fitness at the onset of severe stress than naive communities, even though all communities eventually rescued; and (3) the presence of dispersal, but not its spatial structure, modulated diversification in Phase 1 and evolutionary rescue in Phase 2. We further elaborate on each of these results below.

Effect of Antibiotic Stress and Dispersal on Diversification and Adaptation
During the selection phase, growth, diversity and community composition changed over time but followed different dynamics depending on local environmental conditions, corroborating results from previous studies with this model system (Rainey and Travisano, 1998;Kassen et al., 2000;Massin and Gonzalez, 2006). At the beginning of the selection phase, both growth and diversity were reduced in higher concentrations of streptomycin. This effect on diversity could relate to the nature of the system. Diversification is driven by competition for oxygen in rapidly growing cultures of the aerobic bacteria P. fluorescens, but because streptomycin hampers growth, lower densities at sublethal concentrations could lead to greater oxygen availability, reduced competition, and slower diversification rates. For example, SMT rapidly dominates harsh environments in Phase 1 only in spatially connected communities, while its emergence is much slower in highly stressed, isolated communities. By the end of the selection phase, stressed communities evolved and grew to abundances and diversity levels close to that of benign environments, demonstrating a striking adaptation to high doses of streptomycin within just a few days. Further, we have shown that streptomycin markedly altered community dynamics and dominance patterns among morphotypes. For example, while rare in benign environments, the morph SMT was dominant at intermediate doses of streptomycin (100 or 200 µg/mL) and at the highest dose of streptomycin (400 µg/mL) when communities were also linked by dispersal, potentially indicating a mass effect (Leibold et al., 2004). Other morphotypes such as SMO and WS1, which were dominant in benign environments, occurred in lower abundances in streptomycin-exposed environments, while yet other rare morphs (e.g., fuzzy spreader) went locally extinct after exposure to lethal levels of streptomycin.
These shifts in dominance patterns during the experiment suggest that the potential for resistance evolution is different among morphs. Alternatively, streptomycin might alter the relative competitive abilities of the different morphs as they simultaneously evolve within patches (Rainey and Travisano, 1998) and across the gradient of antibiotic stress between patches (Osmond and de Mazancourt, 2013). Streptomycin normally kills bacteria by binding to the 30S ribosomal subunit, inhibiting mRNA translation and thus protein synthesis (Biswas and Gorini, 1972). Mutations conferring resistance to streptomycin may trade-off with growth potential in benign environments, making resistant individuals less abundant in the absence of streptomycin. Yet another hypothesis for the differences in diversification dynamics concerns the mechanisms by which different morphs emerged in different environments. For example, the transition from SMO-dominated communities to SMT-dominated communities with streptomycin exposure could be the result of both genetic evolution (streptomycininduced selection on standing variation or novel mutations) and phenotypic plasticity. Antibiotic application is known to alter other phenotypes (e.g., biofilm formation) in populations of the closely related Pseudomonas aeruginosa that also acquire antibiotic resistance (Drenkard and Ausubel, 2002). Nonetheless, a strong genetic contribution to phenotypic differentiation has previously been shown for some P. fluorescens morphs, where the emergence of a wrinkly spreader morph from a monoclonal isogenic population of smooth morphs was driven by novel mutations and selection (Rainey and Travisano, 1998;Spiers et al., 2002;McDonald et al., 2009). Further analyses would be necessary to distinguish plastic (e.g., expression of cellulose polymer forming biofilm, Spiers et al., 2002) from genetic influences on morphological diversification and antibiotic resistance in our experimental design, bearing in mind that it is also possible, and perhaps more realistic, that both processes could be involved (Chevin et al., 2013;Kovach-Orr and Fussmann, 2013;Lind et al., 2018;Carja and Plotkin, 2019).
At the metacommunity level, dispersal and streptomycin had opposing effects on both alpha and beta diversity generated during the selection phase. While streptomycin greatly lowered local (alpha) diversity, the effect of streptomycin on community composition also translated into greater compositional difference between local communities. This is expected as the antibiotic gradient creates habitat heterogeneity, promoting beta diversity within metacommunities (Veech and Crist, 2007;Matthiessen et al., 2010). However, dispersal countered the effects of streptomycin, as it increased local diversity in harshly stressed (200 or 400 µg/L) and otherwise depauperate communities, while at the same time homogenizing community composition within metacommunities. Even in control metacommunities without an antibiotic gradient, dispersal increased mean local diversity and reduced metacommunity beta diversity. Both of these effects (higher local diversity, lower beta diversity) of the intermediate rate of dispersal that we used (1%) are consistent with the metacommunity theory (Leibold et al., 2004;Howeth and Leibold, 2010) and with experimental evidence (Matthiessen et al., 2010). Despite these strong effects of dispersal on diversification, the spatial structure of dispersal did not have a consistent effect on diversity; both local and global dispersal resulted in similar patterns of alpha and beta diversity at the end of the selection phase.
Increased diversity through dispersal has the potential to be selected upon in the case of environmental change (Schiffers et al., 2013). We therefore predicted faster or more complete adaptation to the antibiotic stress in metacommunities linked by dispersal. Our results do not support this prediction: neither the presence nor the spatial structure of dispersal influenced community abundance reached by the end of Phase 1. However, the presence of dispersal clearly influenced how fast communities with a history of streptomycin exposure recovered their abundance during the rescue trial of Phase 2.

Drivers of Evolutionary Rescue
We observed repeated and consistent evolutionary rescue in P. fluorescens metacommunities exposed to harsh levels of streptomycin. By the end of Phase 2, viable and diverse communities grew in conditions that were lethal to the ancestral, isomorphic cells, which demonstrates evolutionary rescue of multiple morphs. Because the ancestral population lacked variation, this rescue process was ultimately driven by evolution. Adaptive evolution occurred both during Phase 1 as isomorphic bacteria underwent adaptive radiation driven by competition and modulated by antibiotic exposure, and during the rescue trial in Phase 2 when communities adapted quickly to the high dose of antibiotic. Past exposure to streptomycin was the main cause of the difference in the trajectory of recovery following exposure to severe stress, confirming previous findings from experimental studies of evolutionary rescue in microbial systems (Gonzalez and Bell, 2013;Low-Décarie et al., 2015;Bell et al., 2019;Fugère et al., 2020).
Evolutionary rescue can arise from standing variation or from mutations arising early (i.e., first few cell divisions) in the selection treatment. Our design was focused on the study of the net outcome in the race between the decline in abundance and the rate of recovery by rare resistant cells, or individuals, among different morphotypes of P. fluorescens. Evolutionary rescue is above all a question of whether the rate of adaptation to environmental change occurs on the same time scale as the demographic decline (Gienapp et al., 2013;Lindsey et al., 2013). If all communities eventually recovered similar levels of abundance and diversity, they did so at varying rates.
Following exposure to harsh dose of streptomycin in Phase 2, communities that had been exposed to at least 100 µg/mL of streptomycin in Phase 1 had greater fitness (i.e., population growth averaged across morphotypes) than communities naive to the antibiotic. This effect likely arose because genetic and plastic changes conferring increased resistance at an intermediate dose of stress also provided a correlated advantage at a higher dose of stress. This result supports the conclusion that historical selection facilitates adaptation to a deteriorating environment (Gonzalez and Bell, 2013;Samani and Bell, 2016).
By the end of Phase 2, all communities recovered similar levels of abundance and diversity, perhaps owing to the evolved diversity that was absent in the ancestral cells. Indeed, communities that were allowed to evolve for 7 days even in benign conditions had a much higher morphological diversity than the ancestral cells, variation upon which natural selection could act in Phase 2. However, we did not find a significant effect of diversity on either the trajectory or the outcome of evolutionary rescue. Instead, our results suggest that community composition, rather than diversity, determined the trajectory of evolutionary rescue. Indeed, the abundance of SMT at the end of the selection phase significantly increased the abundance in communities 24 h after the onset of Phase 2 -although only in communities that had also experienced streptomycin in Phase 1. That is, both historical selection and a high relative abundance of SMT were necessary for a swift recovery in Phase 2 (Figure 5C), suggesting that ecological processes (morphotype sorting in favor of the relatively more resistant SMT morph) and evolutionary processes (past adaptation of the SMT morph to streptomycin) both promoted community recovery. Higher abundance of SMT before the rescue trial of Phase 2 was a result of both past exposure to streptomycin and presence of dispersal. Some isolated communities exposed to high doses in Phase 1 had a low relative abundance of SMT, while some control communities with dispersal had a high abundance of SMT -and both these types of communities collapsed to low abundances at the onset of Phase 2. Therefore, our results suggest that community composition influenced the speed of community recovery, as a result of the interplay between past dispersal events and previous exposure to stress. However, as for diversification, the presence of dispersal mattered most, not its spatial structure.
Here, we quantified the outcome of rescue based on the recovery of abundance following exposure to severe stress. By the end of Phase 2, communities recovered abundance and diversity in levels comparable to the end of the selection phase, indicating that evolutionary rescue of multiple morphs occurred in all communities by the end of Phase 2 of the experiment. Abundance or diversity are both aggregate measures of community recovery. Our analysis of community composition showed that communities at the end of Phase 2 were very different in composition compared with Phase 1 communities, indicating a shift to different compositional states despite the recovery of similar levels of abundance and diversity.

CONCLUSIONS AND IMPLICATIONS
Here, we studied the dynamics of evolutionary rescue occurring in evolving communities. We found that the history of stress and dispersal promotes the incidence of evolutionary rescue (Low-Décarie et al., 2015;Bell et al., 2019;Fugère et al., 2020). Our results have management implications, although a few caveats should be acknowledged before extrapolating to natural systems. While dispersal could in theory increase the likelihood of adaptation and subsequent recovery in communities across a fragmented, heterogeneous landscape, other ecological factors such as population size and generation time are crucial to the recovery of populations. Here, the very large populations responsible for striking evolutionary rescue in P. fluorescens may be common for microbial species but are rarely so for metazoan species that typically exist at much lower densities or have slower generation times and lower reproductive rates (vander Wal et al., 2013) although cases of evolutionary rescue in metazoans have been observed (Ozgo, 2014;Reid et al., 2016). This suggests that evolutionary rescue may be less likely in conservation contexts but more likely in agroecosystems or clinical settings where microbes respond to high doses of biocides (Alexander et al., 2014).
We found that dispersal fosters evolutionary rescue following environmental degradation, which corroborates previous findings (Perron et al., 2007;Bell and Gonzalez, 2011;Low-Décarie et al., 2015;Gokhale et al., 2018). For pathogens and invasive species, globalization is increasing the potential for dispersal to contribute to the evolution of resistance across common selective environments, including the widespread application of common antibiotics and pesticides (Thanner et al., 2016;Hudson et al., 2017). For example, streptomycin is used in agriculture worldwide for the control of plant pathogenic bacteria (e.g., in orchards). Resistance is often observed in multiple bacterial pathogens, including Pseudomonas sp. in these orchard settings (Vanneste and Voyle, 2002;Sundin and Wang, 2018), although this may not influence the abundance and diversities of major bacteria taxa soil communities or cause convergence in community similarity (McManus, 2014;Walsh et al., 2014). It remains unclear whether the dispersal of micro-organisms between exposed sites may engender resistance across a microbial metacommunity.
Evolutionary rescue theory is also relevant for conservation biology, where the focus is not the eradication of species but their protection and restoration. Widespread habitat loss is decreasing the potential of already vulnerable populations and communities to disperse across habitat patches and entire landscapes to adapt to new disturbances such as climate change (Norberg et al., 2012). Our results suggest that, especially for depauperate communities, isolation of habitat patches or dispersal barriers could hinder the recovery of communities following severe environmental deterioration (Cheptou et al., 2017). Although we used a simple laboratory model system that lacks some of the complex dynamics and interactions that characterize natural ecosystems exposed to human impacts, our results support the conclusions from theory that spatial connectivity through dispersal is an important determinant of eco-evolutionary dynamics and persistence of diversity across changing and degraded landscapes (Crooks and Sanjayan, 2006;Norberg et al., 2012;Thompson and Fronhofer, 2019).

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.