
==== Front
bioRxiv
BIORXIV
bioRxiv
2692-8205
Cold Spring Harbor Laboratory

38979251
10.1101/2024.06.26.600870
preprint
2
Article
Short-term fluctuating and long-term divergent selection on sympatric Monkeyflowers: insights from decade-spanning reciprocal transplants
http://orcid.org/0000-0003-3352-4006
Dong Caroline M. 12
http://orcid.org/0000-0002-2544-4551
Rolón Bolívar Aponte 1
Sullivan Juj K. 1
http://orcid.org/0000-0003-2632-2262
Tataru Diana 1
Deleon Max 1
Dennis Rachael 1
Dutton Spencer 1
Machado Perez Fidel J. 13
Montano Lissette 1
http://orcid.org/0000-0002-9611-0401
Ferris Kathleen G. 1
1 Tulane University, Department of Ecology and Evolutionary Biology, New Orleans, LA
2 Grinnell College, Department of Biology, Grinnell, IA
3 University of California Merced, Life and Environmental Sciences Department, Merced, CA
AUTHOR CONTRIBUTIONS

KGF conceived the study and CMD, DT, BAR, and JKS contributed to experimental design. CMD, BAR, JKS, MD, RD, SD, FMP, and LM performed the reciprocal transplant experiments. CMD completed data analysis and drafted the manuscript with input from KGF. All authors contributed to and approved the final version of the manuscript.

Corresponding author: dongca@grinnell.edu
03 9 2024
2024.06.26.600870https://creativecommons.org/licenses/by-nc-nd/4.0/ This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which allows reusers to copy and distribute the material in any medium or format in unadapted form only, for noncommercial purposes only, and only so long as attribution is given to the creator.
nihpp-2024.06.26.600870.pdf
Sympatric species are often locally adapted to distinct microhabitats. However, temporal variation may cause local maladaptation and species boundary breakdown, especially during extreme climatic events leading to episodic selection. Repeated reciprocal transplants can reveal the interplay between short and long-term patterns of natural selection. To examine evolutionary trajectories of sympatric Monkeyflowers adapted to different niches, Mimulus guttatus and M. laciniatus, we performed three replicated transplants and combined them with previous experiments to leverage a dataset of five transplants spanning 10 years. We performed phenotypic selection analyses on parents and hybrids in parental habitats in Yosemite NP, CA during years of drastically differing snowpack. If there is ecological isolation, then we predicted divergent phenotypic selection between habitats in line with species’ differences and local adaptation. We found interannual fluctuations in phenotypic selection, often in unpredicted directions. However, a combined-year analysis detected longer-term divergent selection on flowering time, a key temporally isolating and adaptative trait, suggesting that selection may reinforce species boundaries despite short-term fluctuations. Finally, we found temporal variation in local adaptation with M. laciniatus locally adapted in low snowpack years, while an extremely high snowpack year contributed to average local maladaptation of M. guttatus.

adaptation
phenology
speciation
divergence
environmental variation
evolution
natural selection
temporal variation
==== Body
pmcINTRODUCTION

Closely related species often occupy distinct habitats within areas of geographic overlap (Coyne & Orr, 2004; Rundle & Nosil, 2005; Schluter, 2009). Species boundaries in sympatry can be maintained by divergent adaptation to differing environments that generate pre- and/or post-zygotic isolation via selection against immigrants (e.g. habitat isolation) and intermediate hybrids (e.g. extrinsic post-zygotic isolation; Rundle & Nosil, 2005; Coyne & Orr 2004). However, identifying the individual traits that contribute to reproductive isolation is difficult because as speciation progresses, multiple forms of reproductive isolation accumulate between divergent lineages which hinders the measurement of individual barriers (Coyne & Orr, 2004). Once species are partially isolated, neutral processes also cause the accumulation of phenotypic divergence. Therefore, experimental manipulation is necessary to determine whether habitat isolation is a key barrier and which traits are involved in adaptive divergence (Ramsey et al., 2003; Hall & Willis, 2006; Lowry et al., 2008). Measuring phenotypic selection on closely related sympatric species in a reciprocal transplant between native environments can identify whether divergent adaptation is important for reproductive isolation and speciation (Lowry et al., 2008; Anderson et al., 2015). If environmental adaptation is important for reproductive isolation and maintenance of species boundaries, then we expect each species to have higher fitness in its habitat and to detect divergent phenotypic selection in the direction of species’ differences (Anderson et al., 2015), indicating habitat isolation and local adaptation to current conditions.

Despite these predictions, asymmetry in local adaptation, or maladaptation to native habitats, is often observed in reciprocal transplants (Hereford, 2009; Kooyers et al., 2019; Dittmar & Schemske, 2023). This could be due to the prevalence of temporal fluctuations in selection (Siepielski et al., 2009) which can influence the strength of reproductive isolation and species boundaries (Tataru et al., 2023). Additionally, current mean environmental conditions may be different than the conditions during original divergence (Coyne & Orr, 2004; Rundle & Nosil, 2005). If species are no longer adapted to their current niches due to climate change, then we expect to find local maladaptation (Crespi, 2000) and selection for the non-local species’ phenotype (Wilczek et al., 2014; Kooyers et al., 2019; Anderson & Wadgymar, 2020). In this scenario, spatially varying selection is acting contrary to reproductive isolation and could erode species boundaries by driving each towards the other’s phenotypic optimum.

Variation in selection in longitudinal studies is usually measured on a seasonal or inter-annual timescale, a short-term measurement of population evolution. This leads to the question of how short-term fluctuations affect long-term evolutionary change in species and populations (Wadgymar et al., 2017; Dittmar & Schemske, 2023). A synthesis of phenotypic selection studies by Kingsolver & Diamond (2011) found that while temporal variation in selection is common, it does not significantly inhibit long-term directional selection. Therefore, even with fluctuations in the direction of selection, there may be an overall directional or stabilizing trend over a long enough period of time. For example, Wadgymar et al. (2017) found that interannual fluctuations in the direction of selection on the perennial plant Boechera stricta gave rise to long-term patterns consistent with stabilizing selection and local adaptation. Due to the arduousness of temporally replicating reciprocal transplant experiments, it is not currently well understood how short versus longer-term patterns of natural selection compare and impact species’ divergence and reproductive isolation.

Another understudied agent of evolutionary change is episodic selection resulting from sudden and extreme environmental events. Long-term studies which capture the evolutionary responses to these dramatic ecological changes are rare, yet critical for understanding the maintenance of species boundaries (Coltman et al., 1999; Pemberton et al., 2022). Through their long-term study of the Galapagos finches, Grant & Grant (1993) were able to illustrate how an extreme El Niño year altered patterns of selection on beak size and rates of interspecific gene flow. Strong episodes of selection such as this can impact biodiversity by affecting the long-term evolutionary trajectory of species adaptation and reproductive isolation (Grant & Grant, 1993). Extreme climatic events are predicted to become more frequent with continuing anthropogenic climate change (Bailey & van de Pol, 2016).

The Mimulus guttatus species complex, a group of closely related wildflowers, is well-suited for studying the effect of natural selection on species boundaries and reproductive isolation (Wu et al., 2008). Mimulus guttatus and M. laciniatus (syn. Erythranthe guttata and E. laciniata) are locally adapted to contrasting microhabitats throughout the central and southern Sierra Nevada, CA and are incompletely reproductively isolated (Vickery, 1964; Tataru et al., 2024). Mimulus guttatus occurs in moist meadows and seeps while M. laciniatus occurs on nearby rocky granite outcrops (Ferris et al., 2014; Ferris & Willis, 2018). Rocky habitats are harsh with shallow soils, high light intensity, extreme temperatures, and a short growing season (Ferris et al., 2014). Mimulus laciniatus has adapted to its granite habitat with smaller plant size, lobed leaves, self-fertilizing mating system with reduced floral size, and early rapid flowering (DeMarche et al., 2013; Ferris & Willis, 2018). Typical M. guttatus populations are later-flowering, predominantly outcrossing, round-leaved, and have larger plant sizes and flowers (Awadalla & Ritland, 1997; Ferris et al., 2014). Early flowering time is a drought escape strategy and acts as a key pre-mating barrier and contributor to local adaptation in the Mimulus guttatus species complex (Hall & Willis, 2006; Lowry et al., 2008; Ferris et al., 2017; Mantel & Sweigart, 2019). Plant taxa adapted to marginal habitats often have self-fertilizing mating systems which can offer reproductive assurance in a harsh environment where pollinators may be scarce (Wright et al., 2013). Additionally self-fertilizing plants often have small flowers meaning less water is lost through evapotranspiration in floral tissue (Galen, 2000). Lobed leaves may be able to better regulate temperature and water loss through a reduced boundary layer (Nicotra et al., 2011) and are unique to Mimulus species found in rocky habitat (Ferris et al., 2014; Ferris & Willis, 2018; Ferris, 2019).

Here, we investigate spatial and temporal variation in natural selection on interspecific divergence between Mimulus guttatus and M. laciniatus with data from five years of repeated reciprocal transplant experiments under dramatically varying climatic conditions (Figure 1). Previously, hybrids and parents were transplanted into each species’ habitat in Yosemite National Park, CA during years of historically low (Ferris & Willis, 2018) and high snowpack (Tataru et al. 2023). Weaker divergent selection on a key isolating barrier, flowering time, was found during the high snowpack year which could erode species differentiation. We repeated these transplants during similarly low snowpack years (2021, 2022) and a year of exceptionally high snowpack (2023; Figure 1a,b). In this study, we address the following questions: (1) Does the strength and direction of selection on adaptive and reproductively isolating phenotypes fluctuate? (2) How do short versus long-term patterns of divergence compare? (3) Are M. guttatus and M. laciniatus locally adapted? (4) How does episodic selection affect the trajectory of species divergence? By examining both interannual variation and overarching signatures of selection over a 10-year period, we were able to elucidate the effects of a changing climate on the maintenance of biodiversity in this system.

MATERIALS AND METHODS

Repeated reciprocal transplant design

The reciprocal transplants used hybrids created from inbred M. guttatus (YVO 18; 37.723366, −119.746433; 5,395 ft elevation) and M. laciniatus lines (WLF 47; 37.841533, −119.59385; 7993 ft elevation; Figure 1c). Advanced generation hybrids were used because recombination over multiple generations breaks up parental haplotypes and enables us to make inferences about selection on individual phenotypes (Nagy, 1997; Hall & Willis, 2006; Groh & Coop, 2024). Parents were crossed to generate F1 hybrids (WLF 47 maternal × YVO 18 paternal), which were self-fertilized to generate the F2 seeds used in 2021. 500 F2 hybrids were self-fertilized and propagated with single seed descent to produce the F3 seeds used in 2022 and to develop the 300 F6 lines used in 2023. In 2021 and 2022, we pooled seeds from 200 hybrid maternal families for use in the transplant experiment. To prevent inbreeding depression in parental genotypes in 2023, we intercrossed inbred lines within each parental population (WLF and YVO).

Our repeated reciprocal transplants took place in 2021, 2022, and 2023 in Yosemite National Park, CA, USA. We replicated experiments performed in 2013 and 2019 using the same four study sites (Figure 1c): two granite outcrops (Granite 1 [37.810700, −119.485200; 8,436 ft elevation], Granite 2 [37.843702, −119.573120; 8,117 ft elevation]) and two meadows (Meadow 1 [37.767781, −119.772720; 6,605 ft elevation], Meadow 2 [37.755968, −119.803031; 6,163 ft elevation]. In 2021 and 2022, each site comprised of 50 randomized blocks of 18 and 24 individuals, respectively, in a 4×6 formation with one inch spacing between individuals. In 2021, blocks contained 6 F2 hybrids and 6 of each parent and Granite 2 comprised 75 instead of 50 blocks. The 2021 experiment totaled 600 individuals of each genotype (M. guttatus, M. laciniatus, F2 hybrids) in meadow habitat and 750 of each genotype in granite habitat with 4,050 total plants. In 2022, each randomized block contained 16 F3 hybrids with 4 of each parent, totaling 1,600 F3 hybrids and 400 parents per habitat type and 4,800 overall experimental plants. In 2023, local parental genotypes originating from each site were planted only into their origin site, to compare the fitness of local M. guttatus and M. laciniatus parents with the parents of experimental hybrids originating from nearby sites (WLF and YVO). We used 2–3 intrapopulation parental crosses and original inbred lines, distributed evenly across blocks and sites, which did not differ in mean seed count and were combined for downstream analyses (Table S1). Each site had 60 blocks of 35 individuals in a 5×7 formation, composed of 3–4 of parental genotypes (YVO, WLF, and local), and 16–23 F6 hybrids for 200 of each parent and 1,500 F6 hybrids (5 replicates of 300 F6 lines) at each site. There was a total of 1,200 of each parent and 3,000 F6 hybrids in each habitat type with 8,400 overall experimental plants.

We cold-stratified parental M. laciniatus and hybrid seeds for 10 days and parental M. guttatus seeds for 5 days at 4°C to break dormancy and synchronize germination (Sunshine Mix #4 Professional Growing Mix, Sun Gro) for one week in greenhouses at the University of California (UC) Davis (2021) and UC Merced (2022–2023). Following germination, cotyledons were planted into block designs described above. Mortality within the first week was attributed to transplant shock and dead germinates were replaced.

Data collection

To gain an understanding of microhabitat variation between and within sites, environmental measurements were taken weekly from experimental blocks. We measured soil moisture (%) using the SM150T Soil Moisture Sensor (Dynamax), soil surface temperature (°F) using a laser thermometer, and light intensity (μmol m-2 s-1) using a MQ-200X Sunlight Quantum Meter (Apogee Instruments). Further, we recorded survival, herbivory (presence or absence), and adaptive phenotypes previously identified to be under divergent selection (Ferris et al., 2014; Ferris & Willis, 2018; Tataru et al., 2023). At first flower emergence, we recorded the date and measured plant height (to apical meristem; mm), stigma and anther lengths (mm), and corolla width (mm). Flowering time was determined as the number of days from transplant to the first open flower. Herkogamy (i.e. stigma-anther separation) was calculated by subtracting stigma length from longest anther length, as a metric of outcrossing ability. To assess leaf shape, the first true leaf was collected one week after flowering. Leaves were digitally scanned and analyzed to calculate leaf lobing index, a measure of degree of lobing, using ImageJ (Schneider et al., 2012; Ferris et al., 2015).

Assessing fitness metrics

We collected fruits and counted fruit and seed number and calculated their correlation using Kendall’s correlation analysis. We found that selection gradients β values derived from fruit number (2013) were similar to those from seed number (all other years) in both strength and direction (Table S5 and S6) likely due to a strong positive correlation between seed and fruit number (Table S7). Therefore we compare fruit and seed number derived selection gradients across years. Plants that produced flowers but did not produce fruit or seed were indicated with a value of zero. To assess local adaptation, we used three metrics of fitness for parental genotypes: survival to flowering (‘survival’), mean seed number per reproductive individual (‘fecundity’), and mean seed number per planted seed (‘total fitness’). To assess long term patterns of fitness and local adaptation, we included previously published data from two additional years (Ferris & Willis, 2018; Tataru et al., 2023) and averaged survival, fecundity, and total fitness across all years (2013, 2019, 2021, 2022, 2023) for each habitat and species. Data from 2013 was omitted from averaged fecundity and total fitness calculations because only fruit number was counted.

Environmental data analysis

To find associations between environmental variables and plant survival in each habitat (granite or meadow), we used linear mixed effects models with survival as the dependent variable and soil moisture, time, soil temperature and light levels as fixed independent variables, and block as a random effect. Survival in this analysis was calculated as the proportion surviving in each block during each environmental survey. We tested for interactions between independent variables and performed model selection for the best-fit model using AIC selection criteria using R packages nlme v3.1 (Pinheiro & Bates, 2000; Pinheiro et al., 2023) and MuMIn v1.47.5 (Bartoń, 2018). To test whether foliar herbivory and habitat are independent, we used chi-squared tests to test for correlation and created contingency tables. To understand the effect of herbivory on plant fitness, we used ANOVA models using seed count as the dependent variable, and herbivory presence/absence, habitat, site, and block as independent variables and interactions between them, and performed model selection using the R packages stats v3.6.2 (R Core Team, 2023) and MuMIn v1.47.5 (Bartoń, 2018).

Phenotypic selection analyses

To investigate plasticity between habitats in parental genotypes, we compared trait means using t-tests and calculated the coefficient of variation (CV; standard deviation divided by mean; (Pélabon et al., 2020) using the R package stats v3.6.2 (R Core Team, 2023). Although hybrid composition may differ genetically between habitats, we compared trait means of hybrids (aggregated) for a putative indication of plasticity. To investigate linkage between traits, we created a phenotypic trait correlation matrix using hybrid individuals, grouped separately by habitat and year, using R packages Hmisc v5.1–1 (Harrell Jr, 2023) and corrplot v0.92 (Simko, 2021). We found generally weak correlations (r < 0.50) between traits, except between leaf area and plant height or flower width in 2022 and 2023 (Figure S2, S3, S4). Therefore, leaf area was removed from following phenotypic selection analyses.

To understand the strength and direction of selection on traits, we used multivariate linear and quadratic phenotypic selection analyses on hybrid datasets. Quantitative traits (leaf lobing, plant height, outcrossing, flower width, flowering time) were standardized to a mean of 0 and s.d. of 1. We used seed count as the dependent variable with plant height, days to flowering, flower width, herkogamy, leaf lobing, and their interactions as fixed independent variables, and block nested in site as a random effect. We used two types of models to account for overdispersion of zeros in seed count using the R package glmmTMB v1.1.8 (Brooks et al., 2017; Magnusson et al., 2017). Multivariate negative binomial models were used to assess whether a trait is associated with the production of any seeds. Additionally, we used zero-truncated Poisson analysis to assess whether traits are associated with the number of seeds produced (i.e. omitting zero values). We performed model selection for the best-fit model using the R package MuMIn v1.47.5 (Bartoń, 2018). Phenotypes and interactions included in the top model were considered to be under selection. If a trait did not appear in the model, we consider there to have been no selection on that phenotype. To extend the comparison to the 2013 experiment which used only fruit number, we repeated all analyses using fruit number as the dependent variable. Lastly, we repeated phenotypic selection analyses on a combined dataset of hybrid genotypes from all years (2013–2023) using fruit number as the dependent variable with plant height, days to flowering, flower width, leaf lobing, and their interactions as fixed independent variables, and site and year as random effects. To visualize the fitness-phenotype relationships we plotted selection gradients using the R package visreg (Breheny & Burchett, 2017).

RESULTS

Environmental variation causes spatial and temporal variation in survival

Fine-scale environmental variation was significantly tied to survival trends in both habitats across years and fluctuated both spatially between habitats as well as temporally with snowpack (Figure 2, Table S2). Soil moisture availability throughout the season varies significantly between M. laciniatus’s granite and M. guttatus’s meadow habitat in each year with moisture in granite outcrops declining rapidly after an early season plateau, while in meadows there is a gradual, linear decline throughout the growing season. The amount of soil moisture available and the exact shape of the soil moisture curve fluctuated from year to year with higher snowpack years (2019, 2023) showing shallower declines in soil moisture over time in each habitat (Figure 2). Light intensity and temperature are higher on average in M. laciniatus’s outcrop environment except in the extreme high snowpack year of 2023 where temperatures were cooler than usual and similar between habitats throughout the growing season (Figure 2).

Across all years and habitats soil moisture and time (day) were significantly associated with survival, with decreasing soil moisture over time associated with a corresponding decrease in proportion of plants surviving (p < 0.001). In Mimulus guttatus’s meadows, light intensity (p < 0.001) also significantly affected survival in 2021 and 2022, while soil temperature was impactful in 2022 and 2023 (p < 0.001; Table S2). Higher temperatures and light intensities were associated with lower survival. In M. laciniatus’s granite habitat time, soil moisture, light and temperature all impacted survival in the same directions as described for M. guttatus’s meadow habitat in all three recent years (p < 0.0001; Table S2). In 2023, the exceptionally high snowpack year, temperature but not light significantly affected survival (p = 0.5672). Across habitats in 2023, light was not associated with survival. In summary, we found that temperature was significantly associated with survival in meadow habitat in high, but not low snowpack years. In granite, light significantly affected survival in low, but not high snowpack years. Further, 2021–2023 differed from 2019 where soil moisture was the only significant environmental factor (Tataru et al., 2023). Light and temperature were not measured in 2013 (Ferris & Willis, 2018). Similar to 2013 and 2019, herbivory pressure differed significantly between habitats with greater herbivory in M. guttatus’s meadows in 2021 (meadow 22%, granite 1.97%; χ2=552.3, d.f. = 1, p < 0.0001), 2022 (meadow 37.75%, granite 2.58%; χ2=921.78, d.f. = 1, p < 0.0001), and 2023 (meadow 23.17%, granite 8.74%; χ2=325.91, d.f. = 1, p < 0.0001). However, only in 2021 did herbivory and the interaction of herbivory x site significantly affect total fitness in the best-fit ANOVA model (Table S3). Overall, herbivore pressure varied spatially, but its effect on fitness varied temporally.

An extreme climatic event dramatically shifted the growing season

The extreme high snowpack year (2023) differed from other years in several aspects. The growing season was delayed by six weeks or more in both habitats due to the snowpack, with surveys not beginning in granite until previous years’ surveys had ended (Figure 2a,b). Soil moisture declined more gradually than usual in both habitats, while light and surface temperatures were lower (Figure 2c,d). Light intensity was notably lower in both meadow habitats (Figure 2e) and Granite 2 (Figure 2f). In granite in 2023, there was no characteristic temperature spike usually coinciding with mortality (Figure 2h). Unlike in previous years plants were still alive in meadows at the end of the growing season: 41.7% (Meadow 1) and 50.7% (Meadow 2) alive at the last survey before snowfall (Figure 2a). Of the individuals still alive in Meadow 1, many had not flowered indicating a significant life-history delay due to this extreme environmental event (0% of local genotype, 31.0% of hybrids, 44.1% of M. guttatus, 70.7% of M. laciniatus). Meadow 2 was similar (0% of local genotype, 47.3% of hybrids, 40.4% of M. guttatus) except most individuals of M. laciniatus flowered (92.0%). The 2023 growing season was 36–70 days longer than previous years in granite habitats (Granite 2013–2022: 33–87 days; 2023: 112–113 days) but comparable to other years in meadows (Meadow 2013–2022: 74–160 days; 2023: 119–146 days).

Interannually fluctuation versus long-term directional selection

We found temporal and spatial fluctuations in the strength and direction of selection on each quantitative trait in both our negative binomial (any versus zero seeds) and zero-truncated Poisson (seed number excluding zero values) selection analyses (Figure 3). Based on observed trait differences between species, we predicted selection in M. laciniatus’s granite habitat for early flowering time, smaller flowers, increased leaf lobing, shorter plants, and decreased stigma-anther distance (herkogamy), with the inverse direction of selection in M. guttatus’s meadows (Figure 3a–e). While there was divergent selection between habitats in line with the direction of species differences for all traits except flower size in the historic drought of 2013, this was not the case in the other years (Figure 3). Instead, phenotypic selection was frequently in the opposite direction of expectations. The exceptions to this were flowering time in our negative binomial analysis (Figure 3f) and herkogamy (stigma-anther distance) in our zero-truncated Poisson models (Figure 3o). As predicted, there was stronger selection for early flowering in M. laciniatus’s granite habitat in every year but one (Figure 3f, Table 1) and divergent habitat selection on herkogamy with selection for reduced stigma-anther separation in granite in most years (Figure 3o) which matches the native M. lacinaitus’s self-fertilizing mating system. Generally, we found that selection gradients (β) based on seed number (zero-truncated Poisson) were stronger than those derived in our negative binomial analyses (selection on qualitative seed production; Figure 3; Tables 1, 2). For many traits (leaf shape, plant height, herkogamy) selection was more temporally stable in M. guttatus’s meadow habitat than in M. laciniatus’s harsh rocky outcrops (Figure 3m,n,o). Additionally, we found the strongest signatures of phenotypic selection in the extreme high snowpack year (2023) relative to other years which confirms our predictions that this year represents an episodic selection event (e.g. flower width, leaf lobing, plant height, herkogamy; Figure 3; Tables 1, 2).

We saw that temporal fluctuation in selection was associated with variation in snowpack for some traits. For example, in the high snowpack years (2019 and 2023), selection on flowering time was reversed from our predictions with later flowering increasing fitness in granite whereas earlier flowering was advantageous in meadows (Figure 3k, Table 2, S4). For other traits such as flower width, we found that selection was more variable and unpredictable based on snowpack level. There was generally selection for larger flowers with the strongest selection for larger flowers in granite habitat in 2023, contrary to predictions (Figures 3g, l). In fact, the high snowpack and delayed growing season of 2023 (Figure 2) seemed to cause not only atypically strong selection but also selection in unexpected directions, particularly in M. laciniatus’s granite habitat where there was weak selection for later flowering (Figure 3k), strong selection against leaf lobing (Figure 3m), for taller plants (Figure 3n), and bigger flowers (Figure 3l) in that year.

Although we found interannual fluctuations in the strength and direction of selection on all traits, when we combined data across years we saw long term patterns of directional selection (Tables 1, 2; Figure 4, S5). In both zero-truncated Poisson and negative binomial analyses, we found stronger long-term selection for earlier flowering time and taller plants in granite, selection for smaller flowers in both habitats, and long-term selection for more lobed leaves in granite with negligible selection on leaf shape in meadows. Herkogamy was not included in the combined dataset due to its absence in 2013. Therefore, while inter-annual fluctuations in the strength and direction of selection are common in this system, we still find longer-term divergent selection in the direction of species’ differences in both flowering time and leaf shape.

Phenotypes are adaptively plastic between habitats

Parental genotypes exhibited plasticity in phenotypic traits between habitats, with the non-native species exhibiting phenotypic plasticity in the direction of the local species’ phenotype. In all three years, M. laciniatus grew more M. guttatus-like when growing in its non-native meadow habitat with significantly longer flowering time, taller plant height, larger leaf area (Tables S8, S9, S10). These results are similar to those from 2013 with plasticity in every trait except M. laciniatus’s flower size, in the direction of the local species’ phenotype (Ferris & Willis 2018). In its nonnative granite habitat, M. guttatus had significantly shorter mean flowering time and smaller plant height (2021 and 2023) and significantly smaller leaf area and increased leaf lobing like M. laciniatus (2022 and 2023; Tables S8, S9, S10). Similarly, mean hybrid phenotypes were different between habitats with shorter flowering time, smaller plant height, smaller leaf area in hybrids grown in the granite habitat, in the direction of the M. laciniatus (Tables S8, S9, S10). In 2019, hybrids in granite differed by having larger flowers and larger leaf area relative to meadow (Tataru et al. 2023). Overall, an earlier mean flowering time in granite is common throughout years (2013–2023) across hybrids and parents.

Local adaptation fluctuates temporally with long-term maladaptation of M. guttatus

We found temporal fluctuation in the signature of local adaptation between M. laciniatus and M. guttatus parental genotypes. In 2013 we found a clear signature of local adaptation in total fitness, but this was not true for other years in our study (Figure 5). There was also often a difference in local adaptation between the individual components of fitness: survival and fecundity. Mimulus laciniatus had higher survival to flowering than M. guttatus in both habitats in all years, except in 2019 when there was low overall survival of parental genotypes (Figure 5a–e). With all five years combined, M. laciniatus had higher average survival in both habitats (Figure 5f). Patterns of fecundity differed between years with parental species having higher fecundity in their native habitats in 2021 indicating local adaptation and selection against immigrants (Figure 5g–k). However, in 2022, M. laciniatus had higher fecundity in both habitats (Figure 5j), and in 2023, we found reciprocal local maladaptation where M. guttatus had higher fecundity than M. laciniatus in granite and vice versa in meadows (Figure 5k). With all years combined, we found overall evidence of local maladaptation in mean fecundity (Figure 5l); however, if the extreme snowpack of 2023 is excluded then local adaptation in fecundity is seen (Figure S1). For total fitness, M. laciniatus and M. guttatus did equally well in granite in high snowpack years (2019 and 2023; Figure 5n, q), but in low snowpack years M. laciniatus had higher fitness and home-site advantage (2013, 2021, 2022; Figure 5m, o, d). Total fitness in meadow was not as predictable based on snowpack with M. laciniatus having higher total fitness across recent years (2021–2023) but lower fitness in 2013 and 2019. Overall, M. laciniatus had higher total fitness in both habitats (Figure 5r, S1).

DISCUSSION

In this study, we investigated spatial and temporal variation in selection on two sympatric Monkeyflower species, M. guttatus and M. laciniatus. We compared multiple large-scale replicated reciprocal transplants conducted in low (2013, 2021, 2022) and high (2019, 2023) snowpack years including one extreme year resulting in episodic selection (Figure 1a, b). Using experimental hybrids, we found fluctuating selection on phenotypic traits involved in reproductive isolation and environmental adaptation. Within each year, phenotypic selection was often contrary to our predicted directions, but we saw a longer-term trend of divergent selection in the direction of species’ differences in both flowering time, a key reproductive isolating trait, and leaf shape. Our results suggest that local adaptation fluctuates across years but on average, M. guttatus is maladapted to its meadows and that habitat isolation may therefore be asymmetric. Finally, we find that episodic selection, caused by an extreme snowfall event, may disrupt local adaptation and impact species boundaries.

Inter-annual fluctuation in natural selection does not erase signatures of long-term divergence

Temporal and spatial variation in the strength and direction of selection has been well characterized as a mechanism that maintains genetic variation within species (Hereford, 2009; Siepielski et al., 2009, 2017; Bergland et al. 2014; Bitter et al. 2023). However, few studies have empirically tested how fluctuating selection impacts the strength of reproductive isolation between sympatric species (but see Campbell & Powers, 2015). We uncovered annually fluctuating selection on traits between low and high snowpack years, with more consistent patterns of selection on plant height, leaf lobing, flowering time and plant height. There was selection for earlier flowering time in granite across most years (Figure 3). However, in high snowpack years, later flowering was associated with higher seed number in M. laciniatus’s granite habitat which is in the opposite direction of species’ differences (Figure 3). This reflects a relaxation of divergent selection in high snowpack years. Mimulus guttatus that flower later in the season tend to produce higher seed (Mojica & Kelly, 2010; Mojica et al., 2012), particularly during a longer, high snowpack growing season with lowered desiccation risk. Selection on leaf lobing also fluctuated across years but typically with selection patterns opposite predictions (Figure 3).

Annual selection fluctuated unpredictably for flower width with combined selection for smaller flowers (Figure 3; Tables 1, 2), counter to species differences where M. guttatus has larger flowers than M. laciniatus. Flower size affects reproductive success and can contribute to population differentiation and reproductive isolation (Bradshaw et al., 1995; Brunet, 2009; Schiestl & Schluter, 2009; Venail et al., 2010; Krizek & Anderson, 2013). Flower size is often correlated with plant size indicating that early flowering, small stature, and small flowers are linked (Troth et al., 2018). An antagonistic pleiotropy has been demonstrated in M. guttatus where later flowering alleles also increase plant and flower size (Mojica et al., 2012; Monnahan & Kelly, 2015). Although plants with larger flowers have increased fecundity, they also lower viability, resulting in overall selection for smaller flowers (Mojica & Kelly, 2010). Although we saw selection for early flowering and small flowers, we find that selection consistently favored larger plants in both habitats across years suggesting that in our species’ these traits are genetically distinct (Tables 1, 2). Although we initially predicted selection for smaller plant sizes in granite based on M. laciniatus’s short stature, our findings suggest that larger plants always have higher fecundity. Size at reproduction is a key determinant of fecundity across taxa (Aarseen & Taylor, 1992; Aarssen & Jordan, 2001; Kingsolver & Pfennig, 2004; Troth et al., 2018). Despite selection for larger sizes, average plant height was still smaller in granite than meadow habitats (Tables S8, S9, S10), likely constrained by limited water and nutrients in rocky outcrops. We found weak selection on herkogamy perhaps indicating that selection on stigma-anther distance is no longer as strong as during the earlier stages of speciation (Harder & Johnson, 2009).

Even with interannual fluctuations, we found longer-term divergent selection on flowering time and leaf shape between habitats when we combined data across years (Tables 1, 2; Figures 4, S5). In the combined analysis, there was stronger selection for early flowering in granite than meadow habitat which aligns with our predictions based on species’ differences. We also found a distinct signature of divergent selection on leaf lobing in predicted directions. Lobed leaf shape is thought to be involved in adaptation to dry, exposed habitats like M. laciniatus’s rocky outcrop environment (Nicotra et al., 2011; Ferris, 2019) and therefore may also be involved in habitat adaptation and isolation (Ferris & Willis, 2018; Tataru et al., 2023). A signature of overall divergent selection on these key traits in across transplants suggests ongoing adaptive differentiation and the potential for long-term maintenance of species boundaries (Chapurlat et al., 2020).

M. guttatus is locally maladapted

Organisms are expected to have a fitness advantage in their home environment because spatial heterogeneity should favor the evolution of local adaptation (Hastings, 1983). However, temporal changes in environmental conditions are theorized to constrain adaptation if there are opposing selection pressures over time (Stearns, 1992; Kawecki & Ebert 2004) which may negatively impact reproductive isolation and species boundaries. We predicted that if habitat isolation is an important reproductive isolating barrier between our sympatric Mimulus, then there would be a fitness advantage of the local species in each habitat. Instead, we found that local adaptation fluctuated temporally (Figure 5m–q) with evidence of average local maladaptation of M. guttatus (Figure 5r). Annual fluctuations in snowpack was predictive of relative fitness of the parent species in granite with M. laciniatus and M. guttatus performing equally well in high snowpack years (5n,q), but in low snowpack years the drought adapted-M. laciniatus had the home-site advantage (Figure 5m,o,p).

Although local adaptation is thought to be ubiquitous in natural populations (Hereford, 2009), fitness trade-offs are not always found in reciprocal transplants (Bennett & Lenski, 2007; Lowry et al., 2009). For example, in his classic metanalysis Hereford (2009) found that 29% of reciprocal transplant studies failed to detect local adaptation. It is possible that our reciprocal transplant did not detect existing trade-offs due to logistical constraints or that there are weak to no fitness trade-offs associated with adaptation for M. laciniatus across environments (Fry, 1996; Leimu & Fischer, 2008). Alternatively, an increase in dramatic interannual environmental fluctuations due to anthropogenic climate change may mean that M. guttatus is no longer at its adaptive optima. Local adaptation is predicted to evolve in an annual outcrossing plant with large effective population sizes like M. guttatus, but may be constrained by fluctuating selection, gene flow, genetic drift, low genetic diversity, or the genetic architecture of underlying traits (Kawecki & Ebert, 2004; Anderson et al., 2011). Our results could indicate an adaptive lag where, due to recent rapid climate change, the adaptive landscape of M. guttatus’s habitat has changed to be more M. laciniatus-like and M. guttatus’s evolutionary response has not kept pace (Lane et al., 2012; Mills et al., 2013; Kooyers et al., 2019). This pattern of adaptive lag and local maladaptation has been found in other populations of M. guttatus (Kooyers et al., 2019). Over time this trend of adaptive lag could erode species boundaries by allowing colonization of M. guttatus’s meadow by M. laciniatus potentially increasing hybridization or lead to population extinction.

Asymmetry in local adaptation, or a lack of trade-offs, has been found in other transplant experiments (Hereford, 2009; Gosden et al., 2015; Latreille & Pichot, 2017; Toll & Willis, 2018). In a repeated reciprocal transplant between serpentine and sandstone adapted populations of Leptosiphon parviflorus Dittmar & Schemske (2023) found that while serpentine populations had a temporally consistent local fitness advantage, sandstone populations only had a local advantage in two out of four years. While stable spatial heterogeneity should favor the evolution of local adaptation, temporal fluctuation in the environment should select against local specialization and for phenotypic plasticity (Kawecki & Ebert 2004). The average fitness advantage of M. laciniatus across environments could be due to adaptive plasticity (Ghalambor et al., 2015). We found that both M. guttatus and M. laciniatus exhibit plasticity in the direction of the local species’ phenotype in each habitat (Tables S8, S9, S10). Mimulus laciniatus’s occupation of a harsh, highly variable environment (Figures 2 & 3) may have driven the evolution of adaptive plasticity (Via & Lande, 1985) allowing it to thrive in both habitats. However, M. guttatus was more plastic than M. laciniatus in all phenotypes except fitness indicating that plasticity, at least in the traits we measured, is unlikely to account for M. laciniatus’s performance advantage across habitats and the lack of fitness trade-offs.

Episodic selection weakens species divergence

The high Sierra Nevada snowpack of 2023 resulted in the delayed onset of spring in mid to high elevations by approximately 4–6 weeks and dramatically impacted patterns of soil moisture, light intensity, and surface temperature (Figure 2). In shifting the growing season, the extreme snowpack contributed to local maladaptation between parental species (Figure 5j), as evidenced by a signature of local adaptation in fecundity reappearing when 2023 was excluded (Figure S1). Our results indicate that in M. guttatus’s meadow habitat, this delay was detrimental to hybrids and M. guttatus by delaying flowering phenology which negatively impacted fitness. Additionally, none of the experimental local meadow genotypes flowered before snowfall, indicating a negative effect on native populations. However, many meadow plants remained alive at the first snowfall indicating potential facultative perenniality. Shorter day lengths were experienced in the late shifted season and critical photoperiod thresholds may not have been met to initiate the transition flowering (Friedman & Willis, 2013; Fishman et al., 2014; Kenney & Sweigart, 2016). End-of-season drought stress may also typically act as a flowering cue in these populations (Kooyers et al., 2015; Mantel & Sweigart, 2019), however in 2023 soil moisture in the meadows remained higher throughout the season than in previous years (Figure 2). The seasonal delay did not appear to negatively affect the fitness of M. laciniatus, likely because M. laciniatus flowers rapidly under shorter daylengths (Friedman & Willis, 2013; Ferris & Willis, 2018).

In granite, the delayed season in 2023 fostered a more favorable environment for the non-native M. guttatus with high water availability, lower light intensity and surface temperatures, and a longer growing season than is typical in the harsh habitat (Figure 2). This led to higher M. guttatus survival in granite than usual. Although M. guttatus still had lower survival than M. laciniatus, it had remarkably high fecundity (Figure 5e, k) resulting in similar total fitness between the two species in granite which was unprecedented in previous years (Figure 5q). In contrast to the meadow habitat, all plants had senesced and many had set seed by the first snowfall in granite habitat in 2023. This suggests that drought stress is more important than critical photoperiod as a flowering cue in M. lacinaitus’s habitat. It has also been shown that high-elevation populations of M. laciniatus have increased plasticity in critical photoperiod (Love & Ferris, 2024).

Our results demonstrate how an extreme climatic event may disrupt species boundaries. This atypical year produced unusually strong selection opposite of predicted directions on key life-history and morphology traits involved in adaptation and reproductive isolation: flowering time, flower size and leaf lobing (Figure 3). In the long term, episodic selection could expand or contract M. guttatus or M. laciniatus’ ranges, encourage the emergence of hybrid swarms or a new hybrid species (Grant & Grant, 1993), or shift their adaptive potential (Campbell-Staton et al., 2017). The rarity and unpredictability of extreme climatic events complicates predictions about the evolutionary response to climate change. Surveys of plant populations after extreme climatic events can monitor changes in the selective landscape and responses (Anderson, 2016; Bailey & van de Pol, 2016). In species with limited ranges, like the Sierra Nevada endemic M. laciniatus, will these events influence extinction risk? More longitudinal field studies are needed to understand long-term evolutionary consequences of extreme climatic events and changing climatic patterns on native plant species. For instance, do these events drive evolutionary change and does this outweigh selection acting during ‘normal’ periods? For species with a limited range, such as the Sierra Nevada endemic M. laciniatus, will these events influence extinction risk of populations?

Conclusions

Our decade-spanning repeated reciprocal transplants provide greater insight into the impacts of fluctuating selection on species’ boundaries. The foreign advantage of M. laciniatus in the non-native meadow habitat across years suggests either a low cost of adaptation to the granite habitat or recent change in the shape of the adaptive landscape of the meadow habitat leading to current maladaptation of M. guttatus. Furthermore, longer-term divergent selection on flowering time found in our combined phenotypic selection analysis should reinforce temporal and habitat isolation between these species. However, the lack of consistent divergent selection on other traits involved in reproductive isolation and local adaptation could lead to eventual species fusion. Finally, a bout of episodic natural selection due to extreme snowfall shifted population dynamics in a way that may have long term impacts on species’ range, abundance, and gene flow. Our study illustrates the strength of integrating manipulative field experiments across multiple spatial and temporal scales to gain a deeper understanding of the maintenance of biodiversity in the face of environmental change (Wadgymar et al., 2017; Dittmar & Schemske, 2023; Oakley et al., 2023).

Supplementary Material

1

ACKNOWLEDGEMENTS

This work was financially supported by the National Institute of General Medicinal Sciences of the National Institute of Health (NIH) under award number R35GM138224 to KGF. We thank the Yosemite National Park Service for permit support (YOSE-2021-SCI-0033, YOSE-2022-SCI-0051, YOSE-2023-SCI-0011), and Breeanne Jackson at Yosemite Field Station (doi:10.21973/n3v36c) for logistical support and providing accommodation. We thank Johanna Schmitt at the University of California (UC) Davis, and Danielle Edwards and Jason Sexton at UC Merced for providing cold-stratification and greenhouse resources. We also thank members of the Ferris Lab who assisted the project, Natalie Gonzalez and Jessica Scales. Lastly, we acknowledge the original inhabitants of the unceded land on which our research was conducted, the Southern Sierra Miwuk Nation, Bishop Paiute Tribe, Bridgeport Indian Colony, Mono Lake Kutzadikaa, North Fork Rancheria of Mono Indians of California, Picayune Rancheria of the Chukchansi Indians and the Tuolumne Band of Me-Wuk Indians.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study will be openly available in [repository name] at http://doi.org/[doi], reference number [reference number], upon acceptance of the manuscript.

Figure 1. Temporal and spatial scale of the experimental design. Annual changes in a) snow water content and b) accumulated precipitation during experimental years of low snowpack (2013, 2021, 2022) and comparatively higher snowpack (2019 and 2023) with the c) spatial layout of experimental sites: meadow 1 and 2 (filled circles) and granite 1 and 2 (open circles). Nearby parental populations are shown: M. guttatus YVO (filled square) and M. laciniatus WLF (open triangle). Representative meadow and granite habitats are shown in the insets. Water data are averaged from the California Department of Water Resources data stations GIN, WHW, and TUM, located nearby experimental sites along Tioga Road in Yosemite National Park (cdec.water.ca.gov).

Figure 2. Seasonal changes in (a,b) survival (%), (c,d) soil moisture (%), (e,f) light intensity (μmol m-2 s-1), and (g,h) soil surface temperature (°F) in each block per site in meadow and granite habitats in low snowpack years (2013, 2021, 2022) and high snowpack years (2019, 2023). Data from 2013 and 2019 are published in Ferris and Willis 2018 and Tataru et al. 2023 and re-analyzed here.

Figure 3. Visual representation of selection gradients (β) of phenotypic traits measured for hybrids during experimental years in granite (open circles) and meadow (filled circles) habitats. Shown are (a-e) predicted values in each habitat based on native species trait differences and estimated values from (f-j) negative binomial models and (k-o) zero-truncated Poisson models. Significance level is indicated with asterisks and negative graph areas are shaded grey for clarity. The absence of a value during an experimental year indicates that the trait was not present in the best-fit model. All β values are based on seed number, except for values from 2013 which are based on fruit number, as seed number was not recorded. Values from 2013 and 2019 are from Ferris and Willis 2018 and Tataru et al. 2023.

Figure 4. Patterns of cumulative selection with all years combined in meadow (left; filled circles) and granite (right; open circles) habitat in a zero-truncated Poisson analysis on (a,b) flowering time (days from planting), (c,d) width of first flower (mm), (e,f) leaf lobing index, and (g,h) plant height at flowering (mm to apical meristem). Y-axes show adjusted relative fitness (fruit number), statistically corrected for other variables included in the models, and the x-axes show unstandardized trait values for biological relevance. Data from 2013 and 2019 are from Ferris and Willis 2018 and Tataru et al. 2023 and re-analyzed here.

Figure 5. Survival to flowering (a-f), fecundity (mean seed number per reproductive individual; g-l), and total fitness (mean seed number per planted seed; m-r) of Mimulus guttatus (filled square, solid line) and M. laciniatus (open triangle, dashed line) in granite and meadow habitats, from each experimental year and averaged across years (f, l, r). Calculations of fecundity and total fitness for 2013 are based on fruit number because seed number was not recorded; 2013 was omitted from averaged fecundity and total fitness. Values from 2013 and 2019 are from Ferris and Willis 2018 and Tataru et al. 2023.

Table 1. Negative binomial analysis on whether or not a plant produced seeds on individual experimental years 2021–2023 and on a combined dataset of 2013, 2019, 2021, 2022, and 2023.

Trait	Meadow	Granite	
2021	2022	2023	Combined	2021	2023	Combined	
Flowering time (FT)	−0.900***	-	−0.35761**	−0.42795***	−0.6318*	−0.7304**	−0.80732***	
Flower width (FW)	−0.2393	0.551***	0.21311	−0.06862*	−0.2882	0.8929***	−0.09141**	
Leaf lobing (LL)	−0.143	−0.028	−0.06197	-	-	-	-	
Herkogamy (H)	-	-	0.06031	-	-	-	-	
Plant height (PH)	0.903***	0.137	0.66220***	0.60214***	1.461**	2.5706**	0.85300***	
FW × PH	−0.291*	-	−0.25121*	-	−1.057	-	0.11827**	
FT × H	-	-	0.23698*	-		-	-	
LL × PH	−0.232	0.096*	-	-		-	-	
FT × PH	-	-	-	−0.14995***	-	-	−0.56967***	
FW × LL	-	-	−0.48341**	-	-	-	-	
H × PH	-	-	-	-	-	-	-	
FT × FW	-	-	-	−0.09303*	-	-	-	
n =	200	258	775	2014	166	214	992	
Values indicate the selection gradient (β), or strength of selection. Asterisks indicate significance of selection in a trait: ***p < 0.001; **p < 0.01; *p < 0.05. Missing values indicate traits which were not included in the best-fit model. Herkogamy was excluded to from the combined analysis due to its absence from the 2013 dataset. Data from 2013 and 2019 are published in Ferris and Willis 2018 and Tataru et al. 2023 and re-analyzed here.

Table 2. Zero-truncated Poisson analysis, using number of seeds produced if any, on individual experimental years 2021–2023 and on a combined dataset of 2013, 2019, 2021, 2022, and 2023.

Trait	Meadow	Granite	
2021	2022	2023	Combined	2021	2023	Combined	
Flowering time (FT)	−1.827***	−0.289***	−0.19540***	−1.08156***	−0.432	0.1733	−1.12177***	
Flower width (FW)	0.957***	0.178*	−0.09655*	−0.23893***	0.415***	1.2539***	−0.19218*	
Leaf lobing (LL)	−0.374***	0.094**	−0.21639***	0.05601	−0.794***	−4.9779***	0.19775	
Herkogamy (H)	−0.436***	0.167***	0.19504***	-	−1.610***	−0.8008***	-	
Plant height (PH)	0.615***	0.563***	0.57737***	0.96239***	2.15***	−4.0555***	1.05603***	
FT × FW	1.964***	0.542***	0.31703***	−0.21770***	-	−0.8541***	−0.16277	
FT × H	−0.956***	−0.314***	0.31146***	-	−2.022***	-	-	
FT × LL	-	-	-	-	-	-	−0.30463**	
FT × PH	−0.651***	−0.420***	−0.09749*	−0.11394	1.110*	-	−0.34790**	
FW × LL	0.715***	-	−0.42251***	−0.05836	-	1.6669***	0.09612	
FW × PH	−0.601***	−0.154***	−0.12595***	0.07574	-	−1.6953***	0.19489*	
LL × PH	−0.186***	-	0.24088***	0.03987	−1.517***	−8.0169***	0.47229	
H × PH	0.401***	0.072**	−0.07153***	-	0.769***	−1.9776***	-	
FT × FW × LL	-	-	-	-	-	-	0.27679**	
FT × LL × PH	-	-	-	-	-	-	−0.22341*	
FT × FW × PH	-	-	-	0.10875*	-	-	-	
FW × LL × PH	-	-	-	−0.16163*	-	-	−0.30395***	
n =	128	134	202	998	87	51	547	
Values indicate the selection gradient (ß), or strength of selection. Asterisks indicate significance of selection in a trait: ***p < 0.001; **p < 0.01; *p < 0.05. Missing values indicate traits were not included in the best-fit model. Herkogamy was excluded to from the combined analysis due to its absence from the 2013 dataset. Data from 2013 and 2019 are published in Ferris and Willis 2018 and Tataru et al. 2023 and re-analyzed here.

COMPETING INTERESTS

None declared.
==== Refs
REFERENCES

Aarseen LW , Taylor DR . 1992. Fecundity allocation in herbaceous plants. Oikos 65 : 225–232.
Aarssen LW , Jordan CY . 2001. Between-species patterns of covariation in plant size, seed size, and fecundity in monocarpic herbs. Ecoscience 8 : 471–477.
Anderson JT . 2016. Plant fitness in a rapidly changing world. New Phytologist 210 : 81–87.26445400
Anderson JT , Eckhart VM , Geber MA . 2015. Experimental studies of adaptation in Clarkia xantiana: III: Phenotypic selection across a subspecies border. Evolution 69 : 2249–2261.26257193
Anderson JT , Wadgymar SM . 2020. Climate change disrupts local adaptation and favours upslope migration. Ecology Letters 23 : 181–192.31729141
Anderson JT , Willis JH , Mitchell-Olds T . 2011. Evolutionary genetics of plant adaptation. Trends in Genetics 27 : 258–266.21550682
Awadalla P , Ritland K . 1997. Microsatellite variation and evolution in the Mimulus guttatus species complex with contrasting mating systems. Molecular Biology and Evolution 14 : 1023–1034.9335142
Bailey LD , van de Pol M . 2016. Tackling extremes: Challenges for ecological and evolutionary research on extreme climatic events. Journal of Animal Ecology 85 : 85–96.26433114
Bartoń K. 2018. MuMIn: Multi-Model Inference. R package version 1.42.1. https://CRAN.R-project.org/package=MuMIn.
Bennett AF , Lenski RE . 2007. An experimental test of evolutionary trade-offs during temperature adaptation. Proceedings of the National Academy of Sciences of the United States of America 104 : 8649–8654.17494741
Bradshaw HD , Wilbert SM , Otto KG , Schemske DW . 1995. Genetic mapping of floral traits associated with reproductive isolation in monkeyflowers (Mimulus). Nature 376 : 762–765.
Breheny P , Burchett W . 2017. Visualization of Regression Models Using visreg. The R Journal 9 : 56–71.
Brooks ME , Kristensen K , Benthem KJ van , Magnusson A , Berg CW , Nielsen A , Skaug HJ , Maechler M , Bolker BM . 2017. Modeling Zero-Inflated Count Data With glmmTMB. bioRxiv preprint. bioRxiv.
Brunet J. 2009. Pollinators of the Rocky Mountain columbine: Temporal variation, functional groups and associations with floral traits. Annals of Botany 103 : 1567–1578.19414518
Campbell-Staton SC , Cheviron ZA , Rochette N , Catchen J , Losos JB , Edwards S V . 2017. Winter storms drive rapid phenotypic, regulatory, and genomic shifts in the green anole lizard. Science 357 : 495–498.28774927
Campbell DR , Powers JM . 2015. Natural selection on floral morphology can be influenced by climate. Proceedings of the Royal Society B: Biological Sciences 282 .
Chapurlat E , Le Roncé I , Ågren J , Sletvold N . 2020. Divergent selection on flowering phenology but not on floral morphology between two closely related orchids. Ecology and Evolution 10 : 5737–5747.32607187
Coltman DW , Smith JA , Bancroft DR , Pilkington J , MacColl ADC , Clutton-Brock TH , Pemberton JM . 1999. Density-dependent variation in lifetime breeding success and natural and sexual selection in Soay rams. American Naturalist 154 : 730–746.
Coyne JA , Orr HA . 2004. Speciation. Sunderland, Massachusetts: Sinauer Associates.
Crespi BJ . 2000. The evolution of maladaptation. Heredity 84 : 623–629.10886377
DeMarche ML , Rice KJ , Sexton JP . 2013. Niche partitioning between close relatives suggests trade-offs between adaptation to local environments and competition. Ecology and Evolution 3 : 512–522.23531923
Dittmar EL , Schemske DW . 2023. Temporal Variation in Selection Influences Microgeographic Local Adaptation. American Naturalist 202 : 471–485.
Ferris KG . 2019. Endless forms most functional: uncovering the role of natural selection in the evolution of leaf shape. American Journal of Botany 106 : 1532–1535.31797365
Ferris KG , Barnett LL , Blackman BK , Willis JH . 2017. The genetic architecture of local adaptation and reproductive isolation in sympatry within the Mimulus guttatus species complex. Molecular Ecology 26 : 208–224.27439150
Ferris KG , Rushton T , Greenlee AB , Toll K , Blackman BK , Willis JH . 2015. Leaf shape evolution has a similar genetic architecture in three edaphic specialists within the Mimulus guttatus species complex. Annals of Botany 116 : 213–223.26070644
Ferris KG , Sexton JP , Willis JH . 2014. Speciation on a local geographic scale: The evolution of a rare rock outcrop specialist in Mimulus. Philosophical Transactions of the Royal Society B: Biological Sciences 369 : 27–29.
Ferris KG , Willis JH . 2018. Differential adaptation to a harsh granite outcrop habitat between sympatric Mimulus species. Evolution 72 : 1225–1241.29603731
Fishman L , Sweigart AL , Kenney AM , Campbell S . 2014. Major quantitative trait loci control divergence in critical photoperiod for flowering between selfing and outcrossing species of monkeyflower (Mimulus). New Phytologist 201 : 1498–1507.24304557
Friedman J , Willis JH . 2013. Major QTLs for critical photoperiod and vernalization underlie extensive variation in flowering in the Mimulus guttatus species complex. New Phytologist 199 : 571–583.23600522
Fry JD . 1996. The Evolution of Host Specialization: Are Trade-Offs Overrated? The American Naturalist 148 : S84–S107.
Galen C. 2000. High and dry: Drought stress, sex-allocation trade-offs, and selection on flower size in the alpine wildflower Polemonium viscosum (Polemoniaceae). American Naturalist 156 : 72–83.
Ghalambor CK , Hoke KL , Ruell EW , Fischer EK , Reznick DN , Hughes KA . 2015. Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature. Nature 525 : 372–375.26331546
Gosden TP , Waller JT , Svensson EI . 2015. Asymmetric isolating barriers between different microclimatic environments caused by low immigrant survival. Proceedings of the Royal Society B: Biological Sciences 282 .
Grant BR , Grant PR . 1993. Evolution of Darwin’s Finches Caused by a Rare Climatic Event. Proceedings of the Royal Society B: Biological Sciences 251 : 111–117.
Groh JS , Coop G . 2024. The temporal and genomic scale of selection following hybridization. Proceedings of the National Academy of Sciences of the United States of America 121 : 1–11.
Hall MC , Willis JH . 2006. Divergent Selection on Flowering Time Contributes To Local Adaptation in Mimulus Guttatus Populations. Evolution 60 : 2466.17263109
Harder LD , Johnson SD . 2009. Darwin’s beautiful contrivances: evolutionary and functional evidence for floral adaptation. New Phytologist 183 : 530–545.19552694
Harrell F. Jr 2023. Hmisc: Harrell Miscellaneous. R package version 5.1–2.
Hastings A. 1983. Can spatial variation alone lead to selection for dispersal? Theoretical Population Biology 24 : 244–251.
Hereford J. 2009. A quantitative survey of local adaptation and fitness trade-offs. American Naturalist 173 : 579–588.
Kawecki TJ , Ebert D . 2004. Conceptual issues in local adaptation. Ecology Letters 7 : 1225–1241.
Kenney AM , Sweigart AL . 2016. Reproductive isolation and introgression between sympatric Mimulus species. Molecular ecology 25 : 2499–2517.27038381
Kingsolver JG , Diamond SE . 2011. Phenotypic selection in natural populations: What limits directional selection? American Naturalist 177 : 346–357.
Kingsolver JG , Pfennig DW . 2004. Individual-level selection as a cause of cope’s rule of phyletic size increase. Evolution 58 : 1608–1612.15341162
Kooyers NJ , Colicchio JM , Greenlee AB , Patterson E , Handloser NT , Blackman BK . 2019. Lagging adaptation to climate supersedes local adaptation to herbivory in an annual monkeyflower. American Naturalist 194 : 541–557.
Kooyers NJ , Greenlee AB , Colicchio JM , Oh M , Blackman BK . 2015. Replicate altitudinal clines reveal that evolutionary flexibility underlies adaptation to drought stress in annual Mimulus guttatus. New Phytologist 206 : 152–165.25407964
Krizek BA , Anderson JT . 2013. Control of flower size. Journal of Experimental Botany 64 : 1427–1437.23404902
Lane JE , Kruuk LEB , Charmantier A , Murie JO , Dobson FS . 2012. Delayed phenology and reduced fitness associated with climate change in a wild hibernator. Nature 489 : 554–557.22878721
Latreille AC , Pichot C . 2017. Local-scale diversity and adaptation along elevational gradients assessed by reciprocal transplant experiments: lack of local adaptation in silver fir populations. Annals of Forest Science 74 .
Leimu R , Fischer M . 2008. A meta-analysis of local adaptation in plants. PLoS ONE 3 : 1–8.
Love JM , Ferris KG . 2024. Local adaptation to an altitudinal gradient: The interplay between mean phenotypic trait variation and phenotypic plasticity in Mimulus laciniatus. Perspectives in Plant Ecology, Evolution and Systematics 63 : 125795.
Lowry DB , Hall MC , Salt DE , Willis JH . 2009. Genetic and physiological basis of adaptive salt tolerance divergence between coastal and inland Mimulus guttatus. New Phytologist 183 : 776–788.19549130
Lowry DB , Rockwood RC , Willis JH . 2008. Ecological reproductive isolation of coast and inland races of Mimulus guttatus. Evolution 62 : 2196–2214.18637837
Magnusson A , Skaug HJ , Nielsen A , Berg CW , Kristensen K , Maechler M , Bentham KJ van , Bolker BM , Brooks ME . 2017. glmmTMB: Generalized Linear Mixed Models using Template Model Builder. R package version 0.1.3.
Mantel SJ , Sweigart AL . 2019. Divergence in drought-response traits between sympatric species of Mimulus. Ecology and Evolution 9 : 10291–10304.31632643
Mills LS , Zimova M , Oyler J , Running S , Abatzoglou JT , Lukacs PM . 2013. Camouflage mismatch in seasonal coat color due to decreased snow duration. Proceedings of the National Academy of Sciences of the United States of America 110 : 7360–7365.23589881
Mojica JP , Kelly JK . 2010. Viability selection prior to trait expression is an essential component of natural selection. Proceedings of the Royal Society B: Biological Sciences 277 : 2945–2950.
Mojica JP , Lee YW , Willis JH , Kelly JK . 2012. Spatially and temporally varying selection on intrapopulation quantitative trait loci for a life history trade-off in Mimulus guttatus. Molecular Ecology 21 : 3718–3728.22686425
Monnahan PJ , Kelly JK . 2015. Naturally segregating loci exhibit epistasis for fitness. Biology Letters 11 .
Nagy ES . 1997. Selection for native characters in hybrids between two locally adapted plant subspecies. Evolution 51 : 1469–1480.28568619
Nicotra AB , Leigh A , Boyce CK , Jones CS , Niklas KJ , Royer DL , Tsukaya H . 2011. The evolution and functional significance of leaf shape in the angiosperms. Functional Plant Biology 38 : 535–552.32480907
Oakley CG , Schemske DW , McKay JK , Ågren J . 2023. Ecological genetics of local adaptation in Arabidopsis: An 8-year field experiment. Molecular Ecology 32 : 4570–4583.37317048
Pélabon C , Hilde CH , Einum S , Gamelon M . 2020. On the use of the coefficient of variation to quantify and compare trait variation. Evolution Letters 4 : 180–188.32547779
Pemberton JM , Kruuk LEB , Clutton-Brock T . 2022. The Unusual Value of Long-Term Studies of Individuals: The Example of the Isle of Rum Red Deer Project. Annual Review of Ecology, Evolution, and Systematics 53 :327–351.
Pinheiro J , Bates D . 2000. Mixed-Effects Models in S and S-PLUS. ew York: Springer.
Pinheiro J , Bates D , R Core Team. 2023. nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1–163.
Ramsey J , Bradshaw HD , Schemske DW . 2003. Components of reproductive isolation between the monkeyflowers Mimulus lewisii and M. cardinalis (Phrymaceae). Evolution 57 : 1520–1534.12940357
Rundle HD , Nosil P . 2005. Ecological speciation. Ecology Letters 8 : 336–352.
Schiestl FP , Schluter PM . 2009. Floral isolation, specialized pollination, and pollinator behavior in orchids. Annual Review of Entomology 54 : 425–446.
Schluter D. 2009. Evidence for ecological speciation and its alternative. Science 323 : 737–741.19197053
Schneider CA , Rasband WS , Eliceiri KW . 2012. NIH Image to ImageJ: 25 years of Image Analysis. Nature Methods 9 : 671–675.22930834
Siepielski AM , Dibattista JD , Carlson SM . 2009. It’s about time: The temporal dynamics of phenotypic selection in the wild. Ecology Letters 12 : 1261–1276.19740111
Siepielski AM , Morrissey MB , Buoro M , Carlson SM , Caruso CM , Clegg SM , Coulson T , Dibattista J , Gotanda KM , Francis CD , 2017. Precipitation drives global variation in natural selection. Science 355 : 959–962.28254943
Simko WT . 2021. R package ‘corrplot’: Visualization of a Correlation Matrix. (Version 0.92).
Stearns S. 1992. The evolution of life histories. Oxford, UK: Oxford University Press.
Tataru D , De Leon M , Dutton S , Machado-Perez F , Rendahl A , Ferris KG . 2024. Fluctuating selection in a Monkeyflower hybrid zone. bioRxiv: 1–33.
Tataru D , Wheeler EC , Ferris KG . 2023. Spatially and temporally varying selection influence species boundaries in two sympatric Mimulus. Proceedings of the Royal Society B: Biological Sciences 290 .
Toll K , Willis JH . 2018. Hybrid inviability and differential submergence tolerance drive habitat segregation between two congeneric monkeyflowers. Ecology 99 : 2776–2786.30365165
Troth A , Puzey JR , Kim RS , Willis JH , Kelly JK . 2018. Selective trade-offs maintain alleles underpinning complex trait variation in plants. Science 361 : 475–478.30072534
Venail J , Dell’Olivo A , Kuhlemeier C . 2010. Speciation genes in the genus Petunia. Philosophical Transactions of the Royal Society B: Biological Sciences 365 : 461–468.
Via S , Lande R . 1985. Genotype-environment interaction and the evolution of phenotypic plasticity. Evolution 39 : 505–522.28561964
Vickery RK . 1964. Barriers to gene exchange between members of the Mimulus guttatus complex (Scrophulariaceae). Evolution 18 : 52.
Wadgymar SM , Daws SC , Anderson JT . 2017. Integrating viability and fecundity selection to illuminate the adaptive nature of genetic clines. Evolution Letters 1 : 26–39.30283636
Wilczek AM , Cooper MD , Korves TM , Schmitt J . 2014. Lagging adaptation to warming climate in Arabidopsis thaliana. Proceedings of the National Academy of Sciences of the United States of America 111 : 7906–7913.24843140
Wright SI , Kalisz S , Slotte T . 2013. Evolutionary consequences of self-fertilization in plants. Proceedings of the Royal Society B: Biological Sciences 280 .
Wu CA , Lowry DB , Cooley AM , Wright KM , Lee YW , Willis JH . 2008. Mimulus is an emerging model system for the integration of ecological and genomic studies. Heredity 100 : 220–230.17551519
