Mimulus guttatus (syn Erythranthe guttata; 2n = 28; Phrymaceae) is a model for plant genetics and mating systems evolution (Wu et al., 2008). We founded experimental populations by randomly pairing 124 highly inbred lines (Troth et al., 2018). The inbred lines were made by single seed descent (6-12 generations) from a large, random selection of field plants from one natural population of M. guttatus on Iron Mountain in Oregon (44°24′03″N, 122°08′57″W (Willis, 1993; Arathi and Kelly, 2004)). Reciprocal crosses between paired plants produced 62 F₁ families. The F₁ plants were randomly paired and crossed producing F₂ seeds. F₂ seed from all families was combined into a common pool and this pool (our “Ancestral population”) was used to establish 20 experimental populations. We created replicate populations within three treatments: four populations for pure outcrossing, eight populations for mixed mating (10% outcross and 90% selfing) and eight populations for pure selfing.

In each generation, for a total of ten generations (generations 0 through 9 in Figure 1 ), plants were grown following the same protocol. Each of the 20 populations was founded by 30 mg of seed distributed uniformly over the soil surface of a tub. The tubs were contiguous growth arenas (30.5 x 38.1 cm in area) with soil depth of 5 cm. Each tub could draw water through two holes in the bottom. The soil was kept fully hydrated for the first 11 days with a combination of bottom water and top misting. Tubs were rotated every other day, and all populations were treated equivalently in terms of growth conditions throughout the generation. On day 14, fertilizer (250 ml of Jack’s blossom booster (10-30-20) in 25 gallons of water) was applied by bottom watering. There was no further watering until day 21, when all tubs received bottom water for one hour. Flowering typically began around days 21-23. A progressive drought regime was initiated after day 21. Each tub received 1000 ml of water on days 28 and 32, then 500 ml on days 36 and 39, and finally 250 ml on day 42. This progressive drying mimics field conditions where nearly all plants eventually die of desiccation.

In the outcrossing treatment, we enforced outcrossing using hand pollination. All plants flowering on days 24-25 within the outcrossing treatment populations were hand pollinated, randomly choosing plants as pollen donors or recipients. Approximately 60% of the plants had at least one flower by days 24-25 and only one flower was used per plant (an average of about 280 flowers per population). Any incidental self-pollination in outcrossing treatment populations should be quantitatively insignificant because M. guttatus plants from Iron Mountain produce very few seeds from autogamy in the first 24 hours of the floral lifespan even when there is no competing cross-pollen (see Figure 2 in Arathi et al., 2002). In the outcrossing treatment, receiving flowers had opened less than 8 hours before delivery of a saturating dose of cross-pollen via forceps (anthesis occurs between 4am and 8am and hand pollinations were conducted between 8am and noon on the day of anthesis). At the final harvest on days 45-47, only the hand pollinated fruits were collected and all seed was subsequently combined within each population.

In the mixed mating populations, 40 plants were randomly selected to receive hand pollination from another flowering plant on day 26-27. One flower was selected as either a pollen donor or recipient per plant. As described in the outcrossing treatment, the pollination was done in the morning hours between 8am and 12 pm on the select days. The hand pollinated fruits were marked to distinguish them from selfed fruits which were continuously produced until the end of the generation (routinely greater than 5 per plant). At the final harvest (days 45-47), we collected the hand pollinated fruits into one envelope and all others (those produced by selfing) into another. To make the founding seeds for each mixed mating population, we sampled 3mg from the hand-pollination collection and 27mg from the selfed collection (which enforces the 90% selfing condition for this treatment).

In the pure selfing populations, the plants were left undisturbed until the end of the generation. There was no opportunity for outcrossing in the pollinator free environment because the pollen is too large for air transport (and there is no wind anyway). All seeds were collected from the pure selfing populations (one pool per population). As described in the mixed mating populations, the number of selfed fruits per plant varied although it was routinely greater than 5 per plant across populations. At the final harvest (days 45-47) in each generation, all fruits were collected and the seeds combined.

In generation 9, we regrew 60 plants from ancestral seed pool and randomly paired and intercrossed these individuals to make 100 families. As described below, seeds from these families were grown alongside the descendant populations (Generations 10 and 11 in Figure 1 ). This “refresher generation” for the ancestral population eliminates any effects of seed age on phenotype (Franks et al., 2007).

The progeny from the last generation of selection (Generation 10) were grown simultaneously with the ancestral population to maturity without selection. We obtained phenotypic measurements from 1409 plants (217 from the ancestral population and about 60 from each experimental population). On the day that each plant opened its first flower, we recorded the day of flowering and the length of the widest leaf. On the first two flowers produced by each plant, we recorded corolla width and length, pistil length and anther length (Fishman et al., 2002). Stigma-anther separation was the difference between pistil and anther lengths. For analysis, we averaged measurements across the two flowers per plant. The 3^rd flower on each plant was hand pollinated using a saturating amount of pollen from IM767. IM767 is an inbred line derived from the Iron Mountain population with high pollen viability. We grew hundreds of IM767 individuals simultaneously with Generations 10 and 11 to provide a source of genetically identical pollen. Finally, we marked the 4^th flower produced by each plant and allowed it to set seed via autogamy. We eventually collected all seeds from these 3^rd and 4^th flowers for counting and/or weighing.

Variation in seed produced through autogamy involves two components. Many plants fail to produce any seed. Those that set seed exhibit approximately log-normal variation in counts. Considering fruits where we both counted and weighed seed, we find a strong allometric relationship between mass and count ( Supplementary Figure 1 ): Log₁₀(mass) = -0.992 + 0.775 Log₁₀(count). The slope (0.775) is positive but substantially lower than 1, indicating that individual seed mass declines as seed number per flower increases. Given these results, we decided to analyze autogamy as two components. The first is binary (0 or 1) based on whether a plant sets any seed at all (referred to as Prob of self-seed). Second, for plants that set self-seed, we log-transformed the count and denote this value as Log10(self-seed>0). With hand pollination, nearly all plants produced seed, again with approximately log-normal variation in both mass and counts. We calculated Log(mass+1) as reproductive capacity.

Alongside the plants used to take measurements in Generation 10, we grew 64 additional plants per population. These were used as pollen recipients (dams) for crossing from a randomly paired plant from the measured set within each population. The purpose of these crosses was to create outbred seed within each population for growth and measurement in Generation 11 (Kelly, 2005). The Generation 10 plants in the ancestral and outcrossing populations were already fully outbred, but plants within mixed mating and selfing populations were either partially or highly inbred. In Generation 11 (as shown in Figure 1 ), we grew plants from the ancestral population along with individuals from each descendant population (same genotypes grown in Generation 10). These descendant populations were used to generate outbred seed from each population. In addition, we obtained measurements on the same traits as in Generation 10 on a total of 1260 plants.

Generation 10 was initiated on April 4, 2023, generation 11 on August 15, 2023. Each generation was grown as two distinct “cohorts” with seeds distributed to soil 14 days apart (dates above refer to start of first cohort). Within generations, populations were balanced across cohorts. We also grew a collection of Generation 10 seeds from each population simultaneously with the outbred plants of Generation 11 in cohorts 3-4. The balance across cohorts within generations is important to statistically factor out any effect of cohort on trait means. Cohort is included as a causal factor in the ANOVA models (described below) because there are seasonal effects on plants grown in the greenhouse. The within generation analyses (e.g. Figures 2 , 3 below) use only the plants within that generation, while the integrated analysis ( Figure 4 below) considers all plants grown across all four cohorts.

First, to test whether experimental populations diverged from the ancestral population, we applied one-way ANOVAs with each trait as the response and population as the factor. We used Tukey post hoc tests to identify which population means differed significantly from each other. Next, we excluded the ancestors to test for differences among descendant populations. To test whether experimental populations diverged from each other, and whether divergences were predicted by treatment, we applied a nested ANOVA (population within treatment) to each trait measured in Generation 10. Cohort was included as a random factor. We fit a total of 9 Nested ANOVAs to consider all the traits, with population nested in treatment. The same analysis was done on Generation 11 plants. Finally, we analyzed the data from Generations 10 and 11 simultaneously. To distinguish changes due to allele frequency from the direct effects of inbreeding, we used regression with inbred/outbred status of each population included as a covariate. This model also included growth cohort as a categorical predictor. The dependent variable was each trait, independent variables included inbred/outbred status as a covariate. After inspecting the floral trait data, we applied principal component analysis (PCA) to corolla width, anther length, corolla length and pistil length using the correlation matrix. One of our response variables, whether plants produced any seed by selfing, is binary and not continuous. This may cause our ANOVA to be underpowered. For this variable, we also performed Generalized linear model fits (logistic regression) with a binary response and with mating system and cohort included as categorical predictors. All statistical analysis was done using Minitab version 20 (Alin, 2010).

The number of plants that progressed to flower varied among populations, and within population through time, but the average was appoximately 400. Total seed production was always much greater than the 30mg required, usually by more than tenfold (>300mg seed in total). For the plants grown and measured in generations 10 and 11, the means and standard deviations (SD) of each trait in each population are reported in Supplementary Table 1 . As expected, floral dimensions were strongly positively correlated in both generations ( Supplementary Table 2 ). For this reason, we applied principal component analysis to the four flower size measurements (corolla width and length, anther and pistil lengths) and found that PC1, which accounts for 77% of the variation, is a simple average of floral dimensions (loadings reported in Supplementary Table 3 ). PC2 is determined (mainly) by the exertion of both the stigma and anther from face of the flower, while PC3 measures stigma-anther separation. For subsequent tests of trait divergence, we analyzed PC1 as “flower size” and PC2 as “stigma-anther exertion.” For simplicity, we retain the difference (pistil length – anther length) as the trait “stigma-anther separation”. The stigma-anther separation value which had been calculated prior to the PCA exhibits a 0.96 correlation with PC3.

All populations diverged from the Ancestor in one or more traits, and for some traits (e.g., Leaf width), almost all populations differed from the Ancestor ( Supplementary Table 1 ). More importantly, there was a clear effect of the mating system treatment on the pattern of divergence ( Supplementary Table 4 ). The nature of this response differed among traits, with four qualitatively distinct outcomes. The first kind of response, illustrated by days to flower ( Figure 2 , left) and stigma-anther separation ( Figure 2 , right), is where mating systems show consistent differences in trait means that were maintained across both generations 10 and 11. Mean days to flower evolved to the lowest values in the outcrossing populations followed by the mixed mating and selfing populations. Importantly, creating outbred plants in generation 11 within the mixed mating and selfing populations did not alter the relative means of each treatment. This indicates a response driven by allele frequency change (adaptation) and not simply increased homozygosity in the inbreeding populations. Stigma-anther separation showed the opposite pattern to days to flower (the greatest change was in the selfing populations, the least in the outcrossing populations), but the pattern of differences was similar between generations 10 and 11. A previous study showed that both pistil length and anther length exhibit inbreeding depression in this population (Arathi and Kelly, 2003). However, each trait is reduced to a similar extent and so stigma-anther separation, which is the difference (pistil minus anther), has minimal inbreeding depression. In this situation, we predict minimal change from generations 10 to 11, as is observed in Figure 2 , right.

The second type of response is where differences among treatments evident in generation 10 evaporate in generation 11. This was observed for flower size ( Figure 2 , middle) and reproductive capacity ( Figure 3 , middle). In generation 10, flower size and reproductive capacity were significantly reduced in the inbreeding populations relative to the outbred, but the absence of differences in generation 11 implicates homozygosity and not evolution as the cause.

The third sort of response is non-significance in generation 10 but significance in generation 11. Both measures of autogamy (whether a plant sets any self-seed at all and the number of seed if it succeeds in selfing; Figure 3 ) exhibited this pattern, as did leaf width ( Supplementary Table 4 ). For these traits, the differences among mating systems only becomes clear when comparing outbred plants (generation 11) which implies that adaptation and inbreeding depression are having conflicting effects on trait means in generation 10. Plants in the mixed mating and selfing treatments did adapt to become more efficient at selfing, but this improvement is obscured by the general reduction in vigor (and thus in seed production) caused by higher homozygosity. The last category – no difference among mating systems in either generation – was observed only for stigma-anther exertion ( Supplementary Table 4 ). Finally, the variance among populations within treatments was always highly significant, usually greatest in Selfing treatment ( Figures 2 and 3 , Supplementary Table 5 ).

The treatment effects evident in Figures 2 and 3 are based on phenotypic measurements from plants grown simultaneously. We can analyze generations 10 and 11 simultaneously if genetic effects can be distinguished from seasonal differences between generations 10-11. To enable this discrimination of effects, we grew a collection of Generation 10 plants alongside the outbred plants in Generation 11. As a consequence, many genotypes, including the ancestral population, were measured in all four grow-ups (two cohorts within each generation). We use a regression model to distinguish inbreeding effects (present in the mixed mating and selfing populations, but not in the outcrossing or ancestral populations) from cohort-specific environmental effects. The full set of regression coefficients is reported as Supplementary Table 6 . Figure 4 illustrates this model fit for a collection of trait pairs. The change from ancestor (purple symbol in Figure 4 ) to the Generation 10 genotypic mean (red, blue, green circles in Figure 4 ) reflects the combined effect of adaptation and inbreeding depression. The phenotypic “release” from inbreeding depression is the transition from Generation 10 to Generation 11 (arrows to x’s in Figure 4 ). There is no release for the outbred populations because the inbreeding level does not change from generation 10 to 11 (circles and x’s are coincident for these populations).

For most traits, there was for a reduction in trait means from ancestral to Generation 10 populations. For example, both leaf width and reproductive capacity declined but to varying extents with mating system. Leaf width declined most severely in the outbred populations, reproductive capacity in the selfing treatment ( Figure 4A ). The transition to outbred plants largely reverses change in reproductive capacity (significant in Generation 10 but not 11) but exaggerates mating system differences in Leaf width (non-significant in Generation 10 but significant in 11). Flower size exhibits a trajectory similar to leaf width ( Figures 4B, D ). Days to flower ( Figure 4C ) and stigma-anther separation ( Figure 4B ) yield more complicated patterns with the inbreeding release apparently acting in different directions in the selfing and mixed mating treatments. However, inbreeding effects were small for both traits. The major exception to trait declines was autogamy ( Figure 4D ). Inbreeding depression causes declines of means in the treatments where selfing is most important to reproductive success. Adaptation is only evident after comparing outbred plants across treatments.

To understand the results, it was essential to experimentally disentangle adaptation from inbreeding depression ( Figure 4 ). Inbreeding depression is usually considered in relation to fitness and attributed to either deleterious recessive or partially recessive alleles or to overdominance (Charlesworth and Charlesworth, 1987; Dudash and Carr, 1998). Variation in quantitative traits like flower size or phenology may not be caused by rare partially recessive alleles or by overdominance, but that does not imply that they are unaffected by inbreeding. Any quantitative trait locus with dominance (partial or complete) will cause the mean of a trait to change with inbreeding (Lynch and Walsh, 1998). In fact, all else equal, intermediate frequency alleles generate much more inbreeding depression than rare alleles (Kelly, 1999a). In quantitative genetics, this effect is usually termed ‘directional dominance’ instead of inbreeding depression because trait changes are not always negative. Regardless, the trajectories of Figure 4 illustrate how adaptation and inbreeding can drive phenotypic changes in interestingly different ways. In some cases, such as the probability of setting seed by selfing, they are strongly antagonistic. Beyond the present experiment, the conflation of inbreeding depression with adaptation is a concern for any study comparing trait differences among populations that differ in outcrossing rate, either in space or through time.

The trait differences in generation 11 are easier to interpret genetically, they are not more biologically important than the differences evident in generation 10. From a demographic perspective, it is clearly important that pollinator loss increases homozygosity as well as spurring adaptation. Inbreeding depression can greatly reduce absolute individual fitness and thus population viability. Regarding evolution, while the crosses of generation 11 disentangled inbreeding from adaptation, these processes certainly interacted to determine trait changes over the course of the experiment. Selection acts on genetic variation in a different way in outcrossing populations, where allelic effects are an average across heterozygous and homozygous states, than in selfing populations where homozygous effects dominate. Inbreeding also accelerates the stochastic loss of genetic variation due to genetic drift and linked selection (Caballero and Hill, 1992; Caballero and Santiago, 1995; Pekkala et al., 2014; Busch et al., 2022). These processes can cause deleterious alleles, initially at low frequency, to reach fixation within a population. This “fixed load” (Willi et al., 2013) is a consequence of allele frequency change and is not revealed by crosses within a population (as in our contrast of plants from Generations 10 and 11). However, this process accelerates differentiation among replicate populations.

Given that autogamy was a key fitness component in this experiment, rapid evolution was expected. We had also anticipated that stigma-anther separation would decline in selfing and mixed mating treatments – a shorter distance between the anthers and stigmatic surface can increase the probability of pollen transfer (Carr et al., 1997; Carvallo and Medel, 2010; Opedal, 2018). It may also increase “delayed self-pollination” (Goodwillie and Weber, 2018), where self-pollination is deferred to the end of a flower’s lifetime. This is relevant given that Iron Mountain M. guttatus exhibits delayed selfing in a quantitative fashion (Arathi and Kelly, 2004). Across the genus, selfing species routinely have reduced stigma-anther separation relative to outcrossing sister species, e.g. M. nasutus (Fishman et al., 2002), M. micranthus (Carr and Fenster, 1994), M. platycalyx (Dole, 1992), and M. cupriphilus (Macnair, 1989).

The selection responses for other floral traits to pollinator loss/reduction were harder to anticipate. Flower size is typically reduced in selfng species. This may be driven by energetics: Why invest in expensive flowers when attracting pollinators is not required (Cruden and Lyon, 1985)? Alternatively, changes in corolla morphology may be essential to facilitate self-pollination. We observed reductions in flower size in the treatments with selfing, but differences in the means of each treatment group were mostly due to inbreeding depression ( Figure 2 ). The genetic correlation between flower size and days to flower was likely a key factor here. Genotypes that rapidly progress to flower tend to have smaller flowers than genotypes that develop more slowly (Kelly, 2008). Our experimental design (hand-pollination done only on Day 25 in the Outcrossing treatment) clearly imposed selection for rapid progression to flowering and there was a strong response ( Figure 2 ), which produced a correlated response to selection towards smaller flowers in the outcrossing populations. This correlated selection should have been weaker in the Mixed mating and mostly absent in the Selfing treatment populations. Interestingly, all experimental populations evolved smaller flowers relative to the Ancestors (dotted line in Figure 2 ), although perhaps under different pressures.

Studies of contemporary evolution, in the field or greenhouse, can address questions about the ‘order of events’ in the evolutionary transition from outcrossing to selfing. The selfing syndrome consists of traits that range from the direct act of reproduction (e.g. the number of seeds per flower when there is no visitation) to traits affecting that act (e.g. herkogamy and dichogamy) to traits affecting the likelihood of outcrossing as opposed to self-fertilization (e.g. attraction or reward of pollinators) to traits that affect the overall efficiency of reproduction under different mating systems (e.g. pollen/ovule ratios, flowering time). Should these various traits evolve simultaneously or sequentially? It is natural to think that the traits most closely linked to reproduction would evolve first followed by changes in ‘secondary’ characteristics that facilitate attraction of pollinators. Contrary to this view, flower size and floral rewards have shown immediate changes in field populations (Acoca-Pidolle et al., 2023; Bishop et al., 2023), simultaneous with changes in the selfing rate. In this study, we see substantial change not only in the primary traits (autogamy and herkogamy), but still also in the secondary traits like flower size and days to flower ( Figures 2 and 3 ). Current evidence thus suggests that pollinator loss will result in a rapid and multivariate response to selection.