How early does the selfing syndrome arise? Associations between selfing ability and flower size within populations of the mixed-mater Collinsia verna.

PREMISE: Widespread associations between selfing rate and floral size within and among taxa suggest that these traits may evolve in concert. Does this association develop immediately because of shared genetic and/or developmental control, or stepwise with selection shaping the evolution of one trait following the other? If the former, then association ought to appear within and across selfing populations. We explore this fundamental question in three populations of the mixed-mater Collinsia verna where autonomous selfing (AS) ability has been shown to be under selection by the pollination environment.
METHODS: We grew clonal replicates of C. verna in a controlled environment to characterize broad-sense genetic correlations among traits within populations and to assess whether divergence in mating system and floral traits among these populations is consistent with their previously observed selection pressures.
RESULTS: As predicted by their respective pollination environments, we demonstrate significant genetic divergence among populations in AS ability. However, patterns of divergence in floral traits (petal, stamen, and style size, stigmatic receptivity, and stigma-anther distance) were not as expected. Within populations, genetic variation in AS appeared largely independent from floral traits, except for a single weak negative association in one population between flower size and AS rate.
CONCLUSIONS: Together, these results suggest that associations between selfing rate and floral traits across Collinsia species are not reflected at microevolutionary scales. If C. verna were to continue evolving toward the selfing syndrome, floral trait evolution would likely follow stepwise from mating system evolution.

The shift from outcrossing to selfing is one of the most common evolutionary transitions in flowering plants (Barrett et al., 1996; Takebayashi and Morrell, 2001; Igić et al., 2008). Because of the widespread positive association between flower size and outcrossing rate seen among taxa, it is typically assumed that the evolution of the mating system fundamentally coincides with the evolution of floral size traits. Indeed, numerous studies document reduced petal size, decreased stigma‐anther distance, early stigmatic receptivity, and decreased pollen‐to‐ovule ratios, among other traits in more highly selfing species compared to more outcrossing congeners (Lloyd, 1965; Wyatt, 1984; e.g., Ritland and Ritland, 1989; Parker et al., 1995; Goodwillie, 1999; following Foxe et al., 2009; Kalisz et al., 2012), or even among populations within species (Fishman and Stratton, 2004; Whitehead et al., 2018). Although this trait constellation—termed the “selfing syndrome” (Müller, 1883; Ornduff, 1969; Richards, 1997; Sicard and Lenhard, 2011)—has long been recognized (Darwin, 1876), the evolutionary processes shaping the phenotypes associated with the transition to a highly selfing mating system remain unclear.

Two mechanistic explanations for the evolution of the selfing syndrome are constraints, including genetic, developmental, and/or functional, and correlated selection, but these hypotheses are not mutually exclusive. The constraints hypothesis is predicated on a set of shared genes underlying floral traits that synchronizes development and may correlate floral trait values (Krizek and Fletcher, 2005). For example, if the distance between stigmas and anthers (stigma‐anther distance or ‘SAD’) decreases along with flower size because of shared genetic control of styles, filaments, and petals, and if autonomous selfing ability increases with decreased SAD, then selection favoring high autonomous selfing (AS) rates should directly result in decreased flower size. A common scenario that promotes AS is when plants are limited by pollinator service and/or mate availability (the “reproductive assurance hypothesis”: Stebbins, 1957; Jain, 1976; Lloyd, 1979, 1980, 1992; Schoen et al., 1996; Morgan and Wilson, 2005). In this situation, reduced SAD and flower size serves as evidence of an evolutionary response to selection favoring AS ability (e.g., Barrett and Shore, 1987; Dole, 1992; Holtsford and Ellstrand, 1992; Brunet and Eckert, 1998; Chang and Rausher, 1998; van Kleunen and Ritland, 2004; Chen et al., 2009; Brys and Jacquemyn, 2011). However, some work has shown that AS ability and SAD are functionally independent (Carr and Fenster, 1994; Weber and Goodwillie, 2009) and that petal size and SAD are, in part, under independent genetic control (Carr and Fenster, 1994; Fenster and Ritland, 1994a; Lin and Ritland, 1997; Fishman et al., 2002).

Alternatively, selection may directly favor the combination of AS for reproductive assurance and selfing syndrome traits, which could also produce trait associations. For example, because selfing individuals need not rely on pollinators, selfing individuals that redirect resource investments away from attractive floral traits toward seed development can have greater fitness (reviewed in Cutter, 2019). Evidence suggests, however, that the conditions where floral construction or maintenance costs exceed the benefits of attractive, long‐lived flowers may be exclusive to completely selfing species (reviewed in Galen, 1999; and in Obeso, 2002; but see Kelly et al., 2008). For species experiencing unpredictable or even only sporadic opportunities for outcross pollination, a combination of high AS rates and small flower sizes may be less beneficial and could exact high fitness costs through ovule and/or pollen discounting (Holsinger et al., 1984; Lloyd, 1992; Harder and Wilson, 1998). Therefore, species that reproduce via a mixture of selfing and outcrossing (“mixed mating”; e.g., Vogler and Kalisz, 2001) and thus gain a significant proportion of their fitness through pollen export may experience competing selective pressures to maintain both high selfing and outcrossing abilities. In a similar vein, mixed‐mating species would not be expected to experience relaxed selection on floral traits that could also contribute to the selfing syndrome (Cutter, 2019). In either selective scenario (positive or relaxed), we would not expect strong associations between floral traits and AS ability until populations are highly selfing.

The genus Collinsia (Plantaginaceae), presents an exemplary system for exploring the intersection of mating system and floral evolution, because all Collinsia species are self‐compatible and potentially mixed mating (Kalisz et al., 2012). Nevertheless, Collinsia species vary dramatically in flower size and selfing ability, with small‐flowered species generally exhibiting higher population‐level selfing rates. The consistent divergence of floral sizes and associated selfing phenotypes documented for at least six sister species pairs strongly suggests that the selfing syndrome has evolved repeatedly in this genus (Baldwin et al., 2011; Kalisz et al., 2012).

Here, we investigate associations among floral traits and AS ability at the population and genotypic levels in Collinsia verna to evaluate the extent to which such associations occur at microevolutionary scales, even before species and strong mating system divergence. Collinsia verna is a large‐flowered species that diverged from its smaller‐flowered sister C. violaceae < 3.6 mya (Baldwin et al., 2011; Kalisz et al., 2012). Previous work demonstrated that AS is favored and that field selfing rates are higher and positively correlated with the degree of pollinator failure across C. verna populations (Kalisz et al., 2004). Currently unknown is the degree to which the observed pattern of selection among populations has shaped the phenotypic and genetic variation among floral traits and mating systems. To fill this gap, specifically, we ask four questions. (1) Has selfing ability diverged across C. verna populations as expected based on previously observed variable selection for reproductive assurance (Kalisz et al., 2004)? We predict that populations with the most unreliable pollination services and highest in situ selfing rates will have the highest AS rates under common controlled conditions. (2) Do we see evidence of similar among‐population divergence for floral size traits? We expect smaller average flower size in populations with greater AS rates if these traits have evolved in concert. Within populations, (3) what is the extent of standing genetic variation in AS ability and floral traits, and the potential for evolutionary change in these traits, as measured by broad‐sense heritability and evolvability (sensu Houle, 1992; Hansen et al., 2003; Hansen et al., 2011); and (4) is there evidence of broad‐sense genetic correlations between AS ability and floral traits that suggest incipient selfing syndrome evolution?

MATERIALS AND METHODS

Study species

Collinsia verna Nutt., blue‐eyed Mary (Plantaginaceae), is a winter annual found in undisturbed moist deciduous forests of eastern North America (Ohio, Indiana, Kentucky, Tennessee, Illinois, Missouri, and Pennsylvania, USA and Ontaria, Canada). Ranging from mixed mating to highly outcrossing across populations, t  = 0.59–0.93 (Kalisz et al., 1999; Knapczyk, 2007), the outcrossing rate also fluctuates across years within each population (Kalisz et al., 2004; Kalisz et al., 2012). Thus, it is not known to what extent variation in selfing across populations is due to genetic divergence in autonomous selfing (AS) ability. As is characteristic for the genus, flowers occur in whorls in the upper nodes of the stem and are bilaterally symmetric, consisting of two upper and three lower petals. The middle lower petal is folded, forming a keel that encloses the style and four stamens (Baldwin et al., 2011). The keel must be pushed apart by bees to achieve outcrossing. Spatial separation of anthers and stigma is temporally dynamic (termed “movement herkogamy”) because of continued elongation of styles and filaments over the floral lifespan. Filament elongation places sequentially dehiscing anthers at the front of the keel to form a pollen zone (Armbruster et al., 2002), while the style elongates into this zone and becomes receptive (Kalisz, 1989). The timing of stigma‐anther contact, as well as onset of stigmatic receptivity during floral lifespan, varies across Collinsia species, and both occur earlier in highly selfing species (Kalisz et al., 2012). However, these traits can vary within species as well (Lankinen et al., 2007; Randle et al., 2009; Spigler and Kalisz, 2013).

Common environment studies of selfing ability and floral traits

We used a common garden experiment (Clausen et al., 1940) to evaluate the magnitude of genetic divergence among C. verna populations. Genotypes for this study originated from three large (>1 × 106) populations in southwestern Pennsylvania where mates are not limiting. We focused on these three populations because they have been studied extensively with respect to their pollination biology in the field and inbreeding depression (Kalisz and Vogler, 2003; Kalisz et al., 2004). Previous work shows the populations at Braddock Trail (BT) and Enlow Fork (EF) experienced higher pollination failure and exhibited higher selfing rates than at Ten Mile Creek (TMC) (Kalisz et al., 2004). In 1999, the time of previous pollination failure studies, we collected seeds from 75 randomly selected maternal plants in each population. We planted a single seed per family (hereafter “genotype”) in plant growth facilities and germinated seedlings under warm conditions (23°C/12°C day/night with 10 h days). We grew seedlings under cool conditions (13–15°C/12°C day/night with 10 h days) to keep them in a vegetative, nonflowering state for 1–2 months. Once rosettes had 3–4 pairs of true leaves, we generated six clonal replicates (hereafter “clones”) per genotype. To generate the clones, we removed the apical meristem from each genotype to induce the formation of lateral branches. We then excised branches when they were ~1 cm long and planted them individually in trays containing vermiculite that were kept continually moist. After about one month, once roots had formed, we potted each clone into a 4 in plastic pot containing Fafard #4 (Sun Gro Horticulture, Agawam, Massachusetts, USA) potting soil, and three clones per genotype were placed in a randomized block design in each Conviron (model PCG Flex, Pembina, North Dakota, USA) growth chamber (15°C day/10°C night, 12‐hour day) where they remained until fruit set. Additional propagation details can be found in Appendix S1.

To measure AS ability, we isolated three clones per genotype in a growth chamber where they were undisturbed except for watering. The folding of the keel petal to surround the stigma and anthers prevents outcross pollination without the intrusion of insects, particularly bees. There were no insects in our growth chambers, and the controlled environmental conditions (i.e., identical pots, soil, temperature, and watering conditions), suggest there were no differences in resources among individuals that would influence our results. Furthermore, C. verna flowers produce fruits when hand‐pollinated in growth chambers, under the same conditions as used in the current experiment. Thus, any fruits that formed must have occurred via autonomous self‐pollination. Once fruits matured, we counted flower and fruit number on whorls 2–5 and estimated AS rate as the number of fruits divided by the total number of flowers produced (following Charlesworth and Mayer, 1995).

For floral morphological metrics, we placed the remaining three clones per genotype in a pollinator‐free greenhouse. Using the sequential dehiscence of stamens (1–4) to define floral stages within a flower (Kalisz et al., 1999), we collected up to five replicate flowers in each floral stage from whorls 2–5 on all clones and stored them in 70% ethanol until image analysis. We captured digital images of each flower with one upper, one lower, and one side of the keel petal removed to expose all floral parts, and we used ImageJ (National Institutes of Health, Bethesda, Maryland, USA; Rasband, 1997–2018) to measure keel area and the lengths of each stamen, the style, and the upper, lower, and keel petals. From these measurements, we calculated the size and location of the “pollen zone”, which is the floral region bounded by the inner and outermost dehisced anthers (Armbruster et al., 2002). For each stigma, we measured either the shortest distance to the pollen zone or, when within this zone, the mean distance to each dehisced anther to represent SAD.

Using the same three clones per genotype as for floral morphology, we measured the timing of stigmatic receptivity. For each clone, we tested receptivity on three fresh flowers in each of stages 2–4 (C. verna stigmas are strictly nonreceptive at stage 1). We placed the excised style and stigma of each flower between two glass slides, added a drop of hydrogen peroxide, and scored the intensity of bubbling viewed under a dissecting microscope on a scale of 1 (no bubbles for nonreceptivity) to 5 (vigorous bubbling for high receptivity). Occasionally we encountered plants that scored no receptivity for all three stages. To ensure these plants were capable of producing receptive stigmas, we also checked a single stage 5 flower, i.e., postcorolla dehiscence.

Data analysis

The hierarchical structure of our data required a structured analysis (Appendix S1c,d). We averaged replicate measures of floral traits made on each clone prior to analyses. To reduce the number of highly correlated variables, we extracted one principal component (flower size) that represented 99.4% of the variation in all petal traits, and in another analysis, we extracted another principal component (pollen zone) that represented 99.7% of the variation in the size and location of the pollen zone (Appendix S2). In both analyses, we used the basic ‘prcomp’ function without scaling or centering, and without additional rotations. All analyses were conducted in R (R Core Team, Vienna, Austria; R Core Team, 2020).

For analysis of trait divergence among populations, we pooled clone‐level estimates and analyzed genotype mean values because floral traits and AS were measured on different clones. Using the ‘manova’ R function, we performed an analysis of variation (ANOVA) testing the multivariate response variable including flower size, style length, pollen zone, SAD, stigmatic receptivity, and AS for differences among populations. After finding a significant difference among populations in the multivariate trait centroids, we performed univariate ANOVA (‘aov’ R function) for each trait independently followed by post hoc pairwise comparisons using Tukey's honestly significant difference test (Tukey's HSD ‘R’ function). Univariate and multivariate quantile‐quantile plots indicated that all variables were reasonably close to a normal distribution, and collectively they appeared to follow a multivariate normal distribution. Variance was highly similar among populations for each variable (Appendix S3).

To estimate standing genetic variation in AS ability and floral traits, we first partitioned the variation among clones in mean floral trait values (stage 4 only) within and among genotypes to obtain a genetic variance‐covariance matrix (V ). Using clones of the same genetic individual did not allow us to partition genetic variance further into its component parts: additive, dominance, and interlocus interactive variance. Because AS ability and floral traits were measured on different clones, we were also unable to include AS directly in estimating V . Instead, we estimated genetic variation for AS separate from other traits and substituted the genotypic level covariance between AS and all floral traits as a proxy for component‐derived covariances. For a complete genetic variance‐covariance matrix, V , we added a row and column to V , placing AS genetic variance in the diagonal element and the covariances among genotypic trait means and AS in the appropriate rows and columns.

Statistically, this required a multivariate linear mixed‐effects model (‘MCMCglmm’; Hadfield, 2010) to partition the genetic variance‐covariance matrix from the full phenotypic variance‐covariance matrix. The simplest form of this statistical model—called the “animal model”—can be written as . Here, trait for individual is a function of the population mean trait value , the breeding value , and residual error, (Lynch and Walsh, 1998; Wilson et al., 2010). The breeding value estimates V above and is a random effect in the model that is structured by relatedness, which in our case was the genotype. We fit separate models for each population, and for all floral traits and AS together as a multivariate response. To validate model convergence, we graphically analyzed Markov Chain Monte Carlo (MCMC) trace and density plots and verified that all Gelman‐Rubin diagnostic values were ≤1.01. We ran 130,000 iterations, discarded the first 10,000, and saved every 100th sample to reduce autocorrelation among MCMC samples. The effective sample size was at least 1000 for each trait.

We calculated two metrics of broad‐sense evolutionary potential: heritability and evolvability. Whereas broad‐sense heritability, H  = V /(V  + V ), describes genetic determination of phenotypic variation, evolvability, I = V /z 2, can be interpreted as the percent change in population mean expected per unit of directional selection on trait z (Houle, 1992; Hansen et al., 2003). Rather than strict additive effects, our measures of H and I also include dominance and epistatic effects. Previous work in C. verna has shown that maternal effects, which could further inflate estimates of genetic variation, decline rapidly and dramatically through its life cycle (Thiede, 1998).

Finally, we looked for evidence of both phenotypic and genetic correlations between AS ability and floral traits. However, because of the dynamic changes in flower size, pollen zone, and stigmatic receptivity occurring across floral lifespan in Collinsia, we first characterized phenotypic correlations within each flower stage. Floral development in C. verna can be described as movement herkogamy. During development, the short style lengthens, large SAD declines to near zero, and nonreceptive stigmatic tissues become more receptive (Appendix S2) all while four stamens sequentially dehisce. Rather than testing a hypothesis of phenotypic correlations, we focused simply on presenting the patterns of variation, including potential changes in correlations across floral lifespan, to illustrate the dynamics with floral age. To avoid autocorrelation among stages within a flower when testing significance of pairwise genetic correlations involving AS and floral traits, we only used data from fully mature, stage 4 flowers. We estimated Pearson correlation coefficients based on genotypic means to evaluate phenotypic correlations. Although tempting to consider correlations among genotypic means as genotypic correlations, Roff (1997) pointed out scenarios where this may not be true. We computed mean genetic correlations from the posterior modes of the variance‐covariance matrix using the definition of a correlation as a standardized covariance (Via, 1984; Roff, 1997), and we also computed sample genetic correlations for each MCMC sample. To evaluate significance, we calculated the 95% highest posterior density of the marginal distribution for each variance, covariance, and correlation estimate, and we highlighted those intervals that do not include 1, ‐1, or 0.

RESULTS

Population divergence in floral traits and AS ability

Populations differed significantly in floral and mating system trait space (Pillai2,155 = 0.231, F 12,302 = 3.29, P < 0.001; Figure 1). Flower size was not significantly different among populations (MSP = 14.6, F 2,155 = 2.09, P = 0.128, Figure 1A), but there was a significant difference for pollen zone (as a composite trait) among populations (MSP = 4.91, F 2,155 = 5.05, P = 0.008, Figure 1B), which indicates the pollen zone was further from the base of the flower in BT than TMC, while EF was intermediate. Neither style length (MSP = 1.28, F 2,155 = 2.33, P = 0.101) nor SAD (MSP = 0.117, F 2,155 = 0.654, P = 0.521) differed among populations (Figure 1C, D). Stigmatic receptivity differed significantly among populations (MSP = 4.96, F 2,155 = 7.04, P = 0.001), with BT and EF having 17% and 13% greater receptivity scores than TMC (Figure 1E). For the mating system, we document significant genetic differentiation in AS ability among populations (MSP = 0.251, F 2,155 = 7.42, P < 0.001). Populations BT and EF had 38% and 58% greater AS rates relative to population TMC, respectively (Figure 1F). Although flower size was not significantly different among populations as a composite trait, we did find that upper petals (MSP = 13.1, F 2,155 = 14.2, P < 0.001) and lower petals (MSP = 6.36, F 2,155 = 7.76, P < 0.001) were larger in BT and EF than TMC (Appendix S4). Thus, the two populations with higher AS ability also had flowers with longer petals.

Figure 1

Population‐level differences in mean floral traits for C. verna were nonsignificant for (A) flower size, (C) style length, and (D) SAD, according to ANOVA, and significant for (B) pollen zone, (E) stigmatic receptivity score, and (F) the proportion of selfed fruits, i.e., AS. Each sample in these analyses is a genotypic mean, which was first calculated as the mean of all flowers on each clone, then averaged over all clones in a genotype. Means and standard errors are from univariate analysis of variation (ANOVA) models. Different letters above errors bars within a plot indicate significant pairwise differences between populations according to Tukey's Honestly Significant Difference (HSD) test

Genetic variation within populations and evolutionary potential

We found significant evidence of genetic variation within all three populations for all floral and mating system traits (diagonal elements in Table 1), but genetic covariances were mostly not distinguishable from zero. Only pollen zone had statistically nonzero genetic covariance with flower size, style length, and SAD, and only in population TMC. In proportion to the total phenotypic variation, genetic variation ranged from moderate to low for most traits, indicated by the broad‐sense heritability and suggesting relatively low genetic determination in floral and mating trait variation. In BT and TMC, heritability was highest for AS ability (0.38 and 0.43, respectively) followed by SAD (0.33 and 0.32) compared with other floral traits (Table 1). In population EF, heritability was highest for SAD (0.70) followed by pollen zone and style length (0.43 and 0.41). Flower size heritability was notably low, especially in EF. In terms of evolvability, SAD was highest across populations (14–65) and was three orders of magnitude above the next highest trait, AS ability (0.09–0.030; Table 1). Thus, although many traits have comparable levels of broad‐sense heritability, greater evolvability of AS and SAD suggests that we would expect greater proportional responses in those traits relative to flower size traits under a given unit of directional selection.

Table 1

Broad‐sense heritability, evolvability, and genetic variance‐covariance (lower diagonal in bold) and correlations (upper diagonal) among traits in flowers with all four anthers dehisced

Population† Trait‡	H 2	I	Flower size	Style length	Stamen length	Stigma‐anther dist. (SAD)	Stigma recep. (SR)	Auto. selfing (AS) ability
BT
Flower	0.170	0.003	2.064*	0.029	0.253	0.216	0.024	–0.200
Style	0.254	0.004	0.023	0.290*	0.290	–0.469*	0.321	–0.157
Stamen	0.181	0.002	0.197	0.085	0.296*	0.245	0.379	–0.069
SAD	0.321	14.066	0.108	–0.088	0.047	0.122*	–0.107	0.005
SR	0.144	0.014	0.015	0.077	0.092	–0.017	0.200*	0.340
AS	0.383	0.173	–0.114	–0.023	–0.012	0.000	0.048	0.018*
EF
Flower	0.041	0.001	0.546*	0.252	0.101	0.021	–0.109	–0.010
Style	0.411	0.005	0.118	0.398*	0.359*	–0.185	0.398	0.200
Stamen	0.434	0.005	0.068	0.205	0.819*	0.468*	0.379	0.313
SAD	0.699	64.909	0.008	–0.059	0.213	0.252*	0.123	0.134
SR	0.270	0.029	–0.048	0.150	0.205	0.037	0.355*	–0.205
AS	0.339	0.089	–0.006	0.039	0.105	0.016	–0.041	0.012*
TMC
Flower	0.249	0.005	3.246*	0.410	0.788*	0.193	–0.035	–0.082
Style	0.251	0.003	0.359	0.235*	0.373*	–0.533*	0.118	0.030
Stamen	0.289	0.003	0.948*	0.121*	0.446*	0.484*	–0.057	0.029
SAD	0.314	37.076	0.110	–0.082	0.102*	0.100*	–0.181	–0.072
SR	0.253	0.042	–0.042	0.038	–0.025	–0.038	0.439*	0.074
AS	0.433	0.296	–0.038	0.004	0.005	–0.005	0.012	0.017*

Note: Genetic variances and covariances are displayed in the lower triangle in boldface, and genetic correlations are displayed in the upper triangle. Using the 95% credible intervals on the marginal posterior distribution, we evaluated significance of variance‐covariance values as those intervals not including zero and marked them with an asterisk. For genetic correlations, we also tested whether the limits included one. Broad‐sense heritabilities (H ) are listed in the first column. Evolvability (I) is listed in the second column. All correlations were significantly less than 1.0.

Populations = Braddock Trail (BT), Enlow Fork (EF), and Ten Mile Creek (TMC).

Abbreviations for traits in Column 1 refer to trait names in the header.

Floral trait correlations

Flower size was strongly positively correlated with style length and pollen zone location at all flower stages and in all populations but showed little correlation with SAD and stigmatic receptivity (Figure 2). In population EF only, flower size appeared weakly correlated with stigmatic receptivity for stages 2 and 3, and in BT it was weakly negatively correlated with AS (Figure 2). Style length was positively correlated with pollen zone in all stages of floral development and negatively correlated with SAD; SAD was only weakly correlated with pollen zone and only in stage 4. No traits exhibited appreciable correlations with AS ability in any stage of development (Figure 2). Although correlation patterns among stages of floral development were largely similar among populations, there were some minor differences (Figure 2).

Figure 2

Correlation matrix showing patterns of phenotypic variation across floral lifespans. We show Pearson's correlation coefficients among genotypic means for flower size (FS), style length (SL), stamen size (ST), stigma‐anther distance (SAD), stigmatic receptivity (SR), and autonomous selfing (AS) ability in C. verna. Subcolumns contain correlations for each trait pairing within each stage of floral maturity, and subrows contain correlations within each population: Braddock Trail (BT), Enlow Fork (EF), and Ten Mile Creek (TMC). Correlation strength is indicated by color gradient and by printed values for coefficients greater than ±0.3

Across floral stages, the strength of the negative correlation declined between style length and SAD (Figure 2), which we interpret to be—at least in part—as a result of style growth (Figure 3). Shorter styles early in floral development are below the pollen zone and hence have high SAD. The growth of style length and pollen zone width together resulted in SAD near zero in stage 4 flowers (Figure 3). Interestingly, SAD in EF lost its negative correlation with style length in stage 4 flowers and gained a positive correlation with the pollen zone (Figure 2). Although the concurrence of style lengthening with increasing receptivity are normal developmental processes that likely share strong links, these traits were only mildly correlated in a few floral stages, particularly for EF (Figure 2).

Figure 3

Illustration of movement herkogamy showing covariation between style length and stigma‐anther distance (SAD, upper) and SAD versus stigmatic receptivity (SR) score (lower) for each of the three populations of C. verna: Braddock Trail (BT), Enlow Fork (EF), and Ten Mile Creek (TMC). Points represent genotypic means for each four stages of floral maturity, colored from stage 1 in light gray to stage 4 in dark gray. Dashed lines demarcate zero SAD

Genetic correlations

The overall pattern of genetic correlations underlying the floral trait correlations we illustrate was largely similar across populations. While particular pairwise genetic correlations were not universally significant in all populations, the direction of each correlation was consistent among populations. Of the significant correlations we found, the positive correlation between flower size and pollen zone in TMC was the greatest (Table 1). The pollen zone was also strongly correlated with style length in EF and TMC, and with SAD in TMC. Style length was negatively correlated with SAD in BT and TMC. Stigmatic receptivity and AS were both only weakly correlated with other traits (all nonsignificant). Notable among the weak correlations with AS was a negative correlation with flower size in BT that was coupled with a positive correlation between AS and stigmatic receptivity. We might expect higher receptivity in smaller flowers to correlate with AS ability if genetic correlations were guiding evolution toward the selfing syndrome. However, we also found a similar level of evidence to the contrary. The AS in EF had weak positive correlations with style length and pollen zone coupled with a negative correlation with receptivity. None of these weak genetic correlations were statistically significant (Table 1).

DISCUSSION

We show evidence of genetic divergence among Collinsia verna populations in selfing ability and floral traits. For selfing ability, the direction of divergence matches expectations based on prior estimates of pollinator visitation. However, the direction of divergence for floral traits did not match our expectations and could be interpreted as contrary to the selfing syndrome with longer petals coinciding with higher AS ability. Furthermore, the genetic correlations between floral traits and AS ability were very weak and inconclusive. Low correlations among phenotypic traits suggest that floral size, floral developmental rates, and AS can evolve independently in C. verna, and as such does not support a hypothesis of constraints shaping the selfing syndrome. In addition, the lack of an association between flower size and selfing ability could reflect the competing selective pressure for outcrossing exerted by pollinators and selfing for reproductive assurance when pollinators fail to deliver outcrossed pollen in C. verna. Nevertheless, if the microevolutionary patterns seen in this population are applicable to the macroevolutionary patterns seen across the genus Collinsia, then our data suggest that evolution of the mating system precedes that of floral size traits and the associations typical of the selfing syndrome may not yet be apparent while populations retain potentially high outcrossing rates.

An evolutionary response to selection for reproductive assurance in Collinsia verna

Reproductive assurance is considered the leading factor driving the evolution of selfing in plants (Busch and Delph, 2012). Yet, we know of only three studies that have characterized selection on reproductive assurance traits in natural populations (Mimulus: Fenster and Ritland 1994a; Aquilegia: Herlihy and Eckert, 2002; Clarkia: Moeller and Geber, 2005). Previous work in C. verna demonstrated the importance of reproductive assurance via AS by linking increased selfing rates in the field to the increased degree of pollinator failure, which varied among populations and years (Kalisz et al., 2004). Building on this prior work, our common environment study revealed meaningful differentiation in AS ability among these same wild populations, with plants from populations where field‐estimated pollinator failure rates were highest also exhibited higher inherent AS ability. In addition, although inbreeding depression can vary across populations (e.g., Byers and Waller, 1999; Spigler et al., 2017), it is consistently low in these three populations (Kalisz et al., 2004), suggesting that selection by pollinators may be countering the evolution of higher rates of selfing (Spigler and Kalisz, 2017). Linking our work with prior knowledge of the selective environment presents compelling evidence that the mating system has evolved in accordance with the reproductive assurance hypothesis and provides a rare view of an evolutionary response to the mating environment in the wild.

If we were to assume that AS rates were equal among our study populations at some point in the past, then our findings demonstrate an increase in mean AS ability as great as 58%. This magnitude of change would easily be achieved given the amount of genetic variation and evolvability we measured for AS. Even if our estimates of total genetic variation overestimate heritable (additive) genetic variation (Bartkowska and Johnston, 2009) and if our evolvability estimates are somewhat inflated, they predict more change than needed to match the change we observed. Other experiments that have selected for AS report comparable levels (10–32%) of evolvability per unit of directional selection (Bixby and Levin, 1996; Bodbyl Roels and Kelly, 2011). Bixby and Levin (1996) found that just two cycles of selection on the mating system of Phlox cuspidata increased selfing ability 2–10 times that of the original populations. Moreover, our high evolvability estimates suggest the potential for continued mating system divergence among our study populations in future generations if selection pressure on selfing ability continues.

Although heritability for SAD was high relative to other floral traits, we interpret the extremely high evolvability of SAD primarily as an artifact of the small mean value for this trait. Opedal et al. (2017) pointed out that dividing by the mean value is meaningless when calculating a composite trait such as SAD, the difference between pistil and stamen length. We followed their advice to also consider the absolute SAD, but this did not appreciably change our evolvability estimates.

We did not find evidence of incipient associations between selfing ability and floral traits pointing toward the selfing syndrome, finding instead that the evolutionary response of AS ability in these wild populations appears to be independent of floral traits. If floral size traits had evolved in conjunction with the mating system in our study populations, either because of genetic trait correlations or selective trade‐offs between floral size and AS ability, then we would expect populations BT and EF to have smaller flowers and reduced SAD compared to population TMC. Instead, we not only found similar SAD and composite flower size among populations, but when we analyzed petals independently, we found significantly larger upper and lower petals in populations BT and EF, the opposite of the predicted pattern for the selfing syndrome. Notably, the divergence in upper petal size (~5–6%) was relatively small compared to AS (~38–58%), but we also emphasize that the differences we report in AS are not likely due to random genetic processes given that our study populations are extremely large with high estimated levels of genetic polymorphism (Dunn et al., 2006; Knapczyk, 2007) and low levels of divergence among populations within this species in general, F ST = 0.09 (Knapczyk, 2007). Altogether, the disconnect between mating system and floral traits within C. verna suggest that they can evolve independently, supporting similar evidence from artificial selection studies in both Phlox (Bixby and Levin, 1996) and Mimulus (Bodbyl Roels and Kelly, 2011).

A path toward the selfing syndrome

In their natural setting, high outcrossing rates, long floral life spans, and the delayed selfing mechanism in C. verna (Kalisz et al., 1999; Kalisz and Vogler, 2003; Kalisz et al., 2004) indicate that individuals can acquire significant fitness gains via pollen export if and when pollinators are available. Considering also that reduced flower size as a consequence of adaptation to selfing contributes to pollen discounting (reviewed in Busch and Delph, 2012), C. verna may be experiencing selection pressure to maintain, or even increase, flower size. Interestingly, genetic variation for both AS and flower size were highest in TMC, the population with the smallest flowers, lowest AS ability, and least pollinator failure. If natural selection erodes standing genetic variation, then our results could feasibly reflect the results of high‐ pollinator failure in BT and EF having caused selection for both larger, more attractive flowers, and improved AS ability for reproductive assurance. While this idea would require field testing for validation, pollinator preference for larger flowers has been demonstrated for C. verna's congener C. parviflora (Elle and Carney, 2003). A scenario of balancing selection on pollen export (via flower size) in years of high pollinator availability when opportunities for outcrossing are great, and on reproductive assurance (selfing ability) in years of high pollinator failure, might dilute or preclude trade‐offs between selfing ability, inbreeding depression, and flower size in these populations. Moreover, it could result in both high AS ability and large flowers, and more broadly it would contribute to the maintenance of mixed mating and low inbreeding depression in this species.

The key denominator underlying the functional link between AS and flower size is often considered to be SAD. Despite predicting that SAD would be positively and negatively correlated to flower size and AS, respectively, we found no relationships either phenotypically or genetically in any study population. Importantly, this independence was not simply due to a lack of significant genetic variation in SAD. Studies in Leptosiphon (Goodwillie et al., 2006) and Mimulus (Fishman et al., 2002; Kelly and Mojica, 2011) documented independent quantitative trait loci controlling petal size and SAD, which could explain the absence of a correlation between these traits in our results. Furthermore, aspects of C. verna's floral biology could dissolve the functional link between SAD and flower size, including temporal changes in herkogamy and variation in the timing of stigmatic receptivity over floral lifespan (Kalisz et al., 1999; Kalisz et al., 2012). The possibility of stigma‐anther contact at different stages of floral development could facilitate selfing at different stages in different flowers, which would in turn obscure any true relationships among SAD and floral traits when focused on a single floral stage. Nevertheless, SAD did show positive and negative genetic correlations with stamens and styles, respectively, albeit with some inconsistency among populations, but we found no evidence for a genetic correlation between stigmatic receptivity and other floral traits.

In summary, our results suggest that the strong divergence in mating system and floral traits seen across Collinsia species, and in fact across angiosperms at large, are not mirrored at the microevolutionary scale for C. verna. Thus, the association of traits comprising the selfing syndrome may be more likely to occur via selection on the constellation of traits or perhaps even relaxation of selection on floral traits only once selfing reaches fairly high levels. Of the many studies examining the transition to self‐compatibility from self‐incompatible taxa reviewed in Igić et al. (2008) and Goldberg et al. (2010), the pathway involving self‐incompatibility genes was shown to be largely unidirectional. The mechanics of this transition appear to be especially well understood in the Solanaceae (e.g., Markova et al., 2017; Raduski, 2018), where incompatibility, once lost, is rarely regained in a lineage. This may not be the case in mixed mating taxa with reproductive assurance mechanisms as the control point for selfing ability. More focus on mating system evolution in self‐compatible species will shed new light on differences in the paths, processes, and pace of evolution of high selfing relative to self‐incompatible taxa.

AUTHOR CONTRIBUTIONS

D.V. and S.K. designed and carried out the experiment. R.S. and S.K. conducted initial statistical analyses and R.M. completed all statistical analyses and revised the final draft. All authors contributed to writing the manuscript.

Appendix S1. Experimental design.

Click here for additional data file.

Appendix S2. Floral traits and variable reduction.

Click here for additional data file.

Appendix S3. Multivariate analysis of variance.

Click here for additional data file.

Appendix S4. Petal and stamen divergence.

Click here for additional data file.