Phenotypic Selection on Flower Color and Floral Display Size by Three Bee Species

Plants exhibit a wide array of floral forms and pollinators can act as agent of selection on floral traits. Two trends have emerged from recent reviews of pollinator-mediated selection in plants. First, pollinator-mediated selection on plant-level attractants such as floral display size is stronger than on flower-level attractant such as flower color. Second, when comparing plant species, distinct pollinators can exert different selection patterns on floral traits. In addition, many plant species are visited by a diverse array of pollinators but very few studies have examined selection by distinct pollinators. In the current study, we examined phenotypic selection on flower color and floral display size by three distinct bee species, the European honey bee, Apis mellifera, the common eastern bumble bee, Bombus impatiens, and the alfalfa leafcutting bee, Megachile rotundata, foraging on Medicago sativa. To estimate phenotypic selection by each bee species and for all bees combined simultaneously and on the same group of plants, we introduce a new method that combines pollinator visitation data to seed set and floral trait measurements data typical of phenotypic selection study. When comparing floral traits, all bee species selected on the number of racemes per stem and the number of stems per plant, two components of floral display size. However, only leafcutting bees selected on hue or flower color and only bumble bees selected on chroma or darkness of flowers. Selection on chroma occurred via correlational selection between chroma and number of open flowers per raceme and we examine how correlational selection may facilitate the evolution of flower color in plant populations. When comparing bee species, the three bee species exerted similar selection pattern on some floral traits but different patterns on other floral traits and differences in selection patterns were observed between flower-level and plant-level attractants. The trends detected were consistent with previous studies and we advocate the approach introduced here for future studies examining the impact of distinct pollinators on floral trait evolution.


INTRODUCTION
Plants exhibit a high level of floral trait diversity. Flower size, flower color, flower shape, and various aspects of floral display size can vary among plants in a population, among populations of a plant species or among plant species (Brunet, 2009;Dart et al., 2012). The role of pollinators in shaping such floral diversity has been of great interest to evolutionary biologists (Galen, 1996;Fishman and Willis, 2008;Harder and Johnson, 2009;Sletvold et al., 2017). In the last three decades, the attention has focused on identifying the role of pollinators, as opposed to other biotic or to abiotic factors, as agent of selection on floral traits (Strauss and Whittall, 2006;Parachnowitsch and Kessler, 2010;Caruso et al., 2019). The two literature reviews of phenotypic selection in plants have indicated that selection on floral traits by pollinators tend to be greater than by herbivores (Parachnowitsch and Kessler, 2010;Caruso et al., 2019) but can be of similar strength as selection by abiotic factors (Caruso et al., 2019).
To isolate the impact of pollinators on selection of floral traits, it has been suggested to measure phenotypic selection in two groups of plants, one group of hand-pollinated and one group of open-pollinated plants (Fishman and Willis, 2008;Sandring and Ågren, 2009;Parachnowitsch and Kessler, 2010;Sletvold et al., 2017). Selection gradients are estimated for each group and the difference in the selection gradients between the hand-pollinated and the open-pollinated treatments is attributed to pollinator-mediated selection on the floral traits of interest. When concentrating on directional selection for studies that compared hand-pollinated and open-pollinated treatments, two patterns emerged. First, pollinators differentially selected on distinct categories of floral traits. Selection was strongest on floral traits associated with pollinator efficiency such as the length of the corolla tube, followed by plant level traits associated with pollinator attraction such as floral display size and finally selection was weakest for flower level traits associated with pollinator attraction such as flower size and flower color (Caruso et al., 2019). Second, distinct pollinators had different impacts on the selection of floral traits and, among plant species, long-tongue flies or birds tended to exert the strongest selection on floral traits and Lepidoptera the weakest (Caruso et al., 2019).
Few studies have compared selection by distinct pollinators within a plant species (Sahli and Conner, 2011;Worley, 2012, 2013). Conflicting selection among pollinators was identified for some floral traits while for other traits distinct pollinators exerted similar patterns of selection (Sahli and Conner, 2011). For example, in Polemonium brandegeei, hummingbirds selected for stigmas exserted beyond the anthers and for longer and wider corolla tubes while hawkmoths selected for stigmas recessed below the anthers and for narrower corolla tubes Worley, 2012, 2013). These studies examined each pollinator separately and on different sets of plants (Kulbaba and Worley, 2013) or examined pollinators separately and combined in cages (Sahli and Conner, 2011). But plants in natural populations are differentially visited by distinct pollinators whose abundance and efficiency vary and it would be useful to quantify the impact of the major pollinators on the floral traits of interest simultaneously and on the same group of plants.
Pollinator visitation has been used as a proxy for reproductive success in some phenotypic selection studies (Campbell et al., 1997;Zhao and Wang, 2015). Here, we propose to combine pollinator observations with measurements of seed set and floral traits of plants to examine phenotypic selection on floral traits by distinct pollinators. Each plant in the population is expected to receive differential visits by the distinct pollinator species and a different proportion of its flowers will be visited by each species. Such proportions can be used to differentially attribute seeds set on a plant to the distinct pollinators. In addition, data on pollinator efficiency can be combined with flower visits data to proportionally attribute seeds to each pollinator species. Relative fitness (RF) and phenotypic selection can then be measured within each bee species. Phenotypic selection by all pollinators combined can be measured using total seed set per plant to calculate RF. This approach is developed here to illustrate how phenotypic selection by all pollinators combined can be differentially attributed to each pollinator species as we study phenotypic selection on flower color and floral display size by three bee species in Medicago sativa. We determine whether selection by pollinators is stronger for plant-level attractants like floral display size than for flower-level attractants such as flower color (Caruso et al., 2019). We also examine whether the three bee species exert similar or different patterns of selection on these floral traits and how this translates into the overall pattern of selection on the plants. Measuring selection on floral traits by different bee species on the same group of plants provides a more realistic depiction of pollinator-mediated selection on floral traits in plant populations.

Study Species
Medicago sativa L. is an open-pollinated perennial legume that requires bees for seed production. Flowers are clustered into racemes and plants exhibit variation in the number of open flowers per raceme, number of racemes per stem (inflorescence), and number of stems per plant and flower color can also vary, ranging from shades of purple, to white, to yellow (Bauer et al., 2017).
Medicago sativa flowers require tripping for pollination, where pollinators apply pressure to the keel of the flowers which releases its anthers and stigmas. Flowers remain open following tripping but there is little evidence of further pollen deposition by pollinators on already tripped flowers (J. Brunet, pers. obs.). The tripping rate, the proportion of visited flowers that are tripped by a pollinator, varies among bee species (Cane, 2009;Brunet and Stewart, 2010;Pitts-Singer and Cane, 2011). Typically, alfalfa leafcutting bees have the highest tripping rate, followed by bumble bees and finally honey bees (Pitts-Singer and Cane, 2011;Brunet et al., 2019). Honey bees (Apis mellifera) and alfalfa leafcutting bees (Megachile rotundata) are used as managed pollinators in alfalfa seed production fields. In addition, many wild bee species, including the common eastern bumble bee (Bombus impatiens), are known to visit and effectively pollinate alfalfa (Brookes et al., 1994;Brunet and Stewart, 2010).

Experimental Set Up
Five patches of M. sativa with 81 plants per patch, initially planted 0.3 m apart, were set up in a linear arrangement at the West Madison Agricultural Research Station in Madison, WI. One bumble bee hive was set up at the center edge, one honey bee hive 30 m away, and a leafcutting bee domicile was set up at the northwest corner facing southeast. About 1.2 lbs of leafcutting bees were released prior to M. sativa peak bloom. Floral trait and fitness measurements were obtained from all flowering plants in two center patches where each plant was numbered.

Floral Traits
The floral traits examined in this study included components of floral display size and flower color. For floral display size, we recorded the number of stems per plant, racemes per stem and open flowers per raceme. For each flowering plant, we counted the number of stems and the number of racemes per stem on ten randomly selected stems or on all stems if a plant had fewer than ten stems. The number of open flowers per raceme was recorded on ten randomly selected racemes per plant. The average number of racemes per stem and open flowers per raceme were tabulated for each plant.
Flower color was determined from spectral measurements of the banner petal for three flowers per plant using the USB 4000 spectrophotometer (Ocean Optics, Orlando, Fl., 350-1,000 nm). Reflectance data were analyzed using Spectra Suite v.10.7.1 software (Ocean Optics). Flowers of M. sativa do not reflect in the UV range (Bauer et al., 2017) and spectral measurements were taken in the visible light range (400-700 nm). We used equations from Endler (1990) as modified by Smith (2014) to calculate three components of flower color: chroma (darkness or saturation), hue (color), and reflectivity (brightness). Details of these calculations can be found in Bauer et al. (2017). A plant value represented the average of the three flowers.
An alternative would have been to use a hexagon color vision model, a method that considers bee photoreceptors when quantifying color (Chittka, 1992;Chittka and Kevan, 2005). We have used such models to examine how flower color affected the choice of plants by bees (Bauer et al., 2017). However, in this study, while bees may be doing the selection, they are selecting on the plant traits and not on their perception of those traits. We thus chose hue, chroma and receptivity to describe flower color. While the best method to quantify flower color when pollinators are selecting on the trait may deserve further attention, such discussion is beyond the scope of this study.

Female Reproductive Success
We used the total number of seeds produced per plant as a measure of female reproductive success. On each plant, on ten randomly selected stems or all stems if a plant had fewer than ten stems, we counted the number of pods per stem. A pod is a fruit developing from one flower on a raceme. We collected ten randomly selected fruiting racemes per plant and placed each one in an individually marked paper coin envelope. In the laboratory, the number of pods per raceme were recorded and pods were shredded to obtain the number of seeds per raceme. For each plant, we obtained the average number of mature seeds per pod per raceme and, using the 10 fruiting racemes per plant, we calculated the average number of mature seeds per pod on a plant. To obtain the total number of seeds set per plant, we multiplied the average number of pods per stem by the average number of seeds per pod and multiplied this value by the number of stems produced on a plant.

Proportion of Seeds Attributable to Each Bee Species
To estimate the proportion of seeds on each plant attributable to a given bee species, we used available data on the number of flower visits to a plant by each of the three bee species. Pollinator visitation data were collected on these plants during a two-week observation period at peak bloom for M. sativa the year of the study (Bauer et al., 2017). To determine the number of pollinator visits to a plant, we followed bees in a patch and two observers recorded each plant visited by a bee, the number of racemes visited on a plant and flowers visited per raceme on each plant until the bee left the patch or was lost to the observers. This provided floral visits by at least one of the three pollinators for most plants in the two patches (Supplementary Data). The pattern of visitation in the patches was typical for the major bee species visiting M. sativa throughout its flowering period. To attribute the number of seeds to a bee species based on the number of flower visits, for each plant, we multiplied the proportion of flowers visited by each of the three bee species by the number of seeds set on that plant. This approach links, for each plant, the pollinator visitation data to its seed set during that period, as seeds were collected about four weeks following pollinator observations, the period it takes for fruits and seeds to mature in this plant species.
The number of flowers visited by a bee species is a useful measure of pollinator visits, but to better link floral visits to seed set we also integrated the tripping rate of a bee species to the floral visitation data. In M. sativa, flowers must be tripped before they can produce seeds and tripping rate varies among bee species (Cane, 2009;Pitts-Singer and Cane, 2011). We obtained a second measure of pollinator visits which combined floral visits with the tripping rate of a bee species. Previous observations in the area indicated a tripping rate of 55% for bumble bee, 25% for honey bee (Brunet and Stewart, 2010) and 80% for leafcutting bee under warm temperatures typical of alfalfa seed-production fields (Brunet et al., 2019). For each plant, the number of flowers visited by a bee species on that plant was multiplied by the bee species specific tripping rate. We call this measure the number of flowers tripped by a bee species. For each plant, we calculated the number of tripped flowers by each bee species and the proportion of flowers tripped by each bee species. We multiplied these proportions by the number of seeds set on the plant to assign seeds to each of the three bee species based on the number of tripped flowers.

Plant Relative Fitness
Plant relative fitness (RF) was estimated by dividing the absolute fitness of a plant by the mean absolute fitness of the group of plants under consideration (Lande and Arnold, 1983). The absolute fitness of a plant was quantified as the number of seeds set on a plant. RF was obtained for all plants for which floral trait measurements were available (N = 153). We calculated RF of a plant over all bees, based on the total number of seeds it produced, and within each bee species. Within a bee species, RF was the number of seeds on a plant attributable to a given bee species, based either on the proportion of flowers visited or the proportion of flowers tripped by a specific bee species, divided by the mean for that bee species. Using this approach, the mean RF was 1.0 within each bee species and potential differences in seed production across pollinators were removed. We also calculated the opportunity for selection for overall RF and for RF by bee species based on proportion of visits or proportion of tripped flowers. Opportunity for selection was measured as the variance in RF.

Phenotypic Selection
To measure phenotypic selection, we examined the relationship between the trait value of a plant and its RF (Lande and Arnold, 1983). Each floral trait examined was scaled such that its mean was 0 and its variance was 1: (trait valuetrait mean)/trait standard deviation. We performed phenotypic selection analyses on the number of stems per plant, the number of racemes per stem, the number of open flowers per raceme, and hue, chroma and reflectivity. We first performed phenotypic selection analyses using RF calculated over the total seed set of a plant. We also examined phenotypic selection within each bee species, where RF was calculated as explained earlier, either based on proportion of flowers visited or proportion of flowers tripped by a bee species on each plant. RF was relativized and traits standardized within each bee species which eliminated any potential differences in traits or fitness across bee species. The number of plants was similar for overall fitness and within each bee species and represented plants with floral traits and seed set data (N = 153). The number of plants that received no visits by a specific bee species did vary, with a greater number of plants not visited by leafcutting bee (N = 107 plants), followed by honey bee (N = 46) and last bumble bee (N = 16).
We used regression analyses, examining linear and nonlinear regressions, to estimate various selection parameters, following the methods suggested by Lande and Arnold (1983). Untransformed variables were used to obtain the values of the selection coefficients. To obtain the statistical significance of the selection coefficients, RF values were log transformed in order to improve the model's residuals. This procedure was followed because selection coefficients are not known to be affected by a poorly fit model while the probability values are (Lande and Arnold, 1983;Mitchell-Olds and Shaw, 1987;Brodie and Janzen, 1996). In addition, due to the large number of zeros, the model's residuals for leafcutting bee still indicated a poor fit to the data after transforming RF. We therefore used bootstrapping to estimate the 95% confidence intervals around the selection coefficients and determine whether they were statistically significant (Davison and Hinkley, 1997). We used bootstrapping for all cases for comparison purpose. We performed 1,000 bootstraps using the bootstrap function in the package "boot" (Canty and Ripley, 2020) in R (version 3.6.1).
For directional selection, we estimated the selection differential (S i ), which represents the change in the population mean of trait (i) after selection (Arnold and Wade, 1984). The selection differential can be obtained from the slope of a linear regression between the standardized value of a trait and the corresponding plant RF. This coefficient includes both direct and indirect selection and multiple regression analyses were performed to isolate direct selection. The partial regression coefficient for a trait represents the selection gradient (βi) for that trait (i) and illustrates direct selection on a trait after removing indirect selection from all other traits present in the analysis. When traits are correlated, a trait that appears to respond to selection may simply be correlated to the trait under selection, hence the need to isolate direct selection. The coefficients S or β both represent directional selection and a positive value indicates that the phenotypic mean of a trait (i) increases under selection while it decreases when values of Si or βi are negative.
Because selection can also be non-linear and work on the shape of the trait distribution, we first added a quadratic term to the single regression and obtained the non-linear (quadratic) selection differential C ii (Table 1), where C 22 illustrates the non-linear (quadratic) term of the single regression. We then performed multiple regressions with linear, quadratic and cross product terms to obtain the non-linear or quadratic selection gradient γ ii , represented by the partial regression coefficient for the quadratic term, and to detect correlational selection γ ij using the partial regression coefficient for the cross product terms ( Table 1; Brodie, 1992;Roff and Fairbairn, 2012). The quadratic coefficient gradients were estimated as double the quadratic regression coefficients (Stinchcombe et al., 2008;Sahli and Conner, 2011).
We graphically illustrated the statistically significant cross product terms representing correlational selection gradients using the function "persp" in R (R Core Team, 2019). To represent the non-linear selection for the statistically significant quadratic selection gradients we used generalized additive models (GAMs) using the "mgcv" package in R (Wood, 2011). These models automatically fit a spline regression (Wood, 2011). Results from the GAMs were plotted using "ggplot2" (Wickham, 2016) and "gridExtra" (Auguie, 2017).
Besides using regression analyses, we also examined the distributional selection gradient on the floral traits (Henshaw and Zemel, 2017). This measures total selection on a trait and can be broken down into a directional component (dD) illustrating selection on a trait mean and a non-directional component (dN) that reflects selection on the shape of the trait distribution. This approach permits estimation of the general selection differential (S) and selection gradients (β). We used the R code available from 1 | Selection parameters obtained based on different regression analyses using plant relative fitness and standardized floral traits.

Model
Single regression Multiple regression Linear S Selection differential; slope; both direct and indirect selection β Selection gradient; partial regression coefficient; direct selection Non-linear (linear and quadratic terms) C ii Non-linear or quadratic selection differential is C 22 Non-linear (linear, quadratic and cross product terms) γ ii Non-linear or quadratic selection gradient; partial regression coefficient of quadratic term Non-linear (linear, quadratic and cross product terms) γ ij Correlational selection gradient; Partial regression coefficient of cross product term The subscript i indicates trait(i) and j indicates a separate trait. The methodology followed to obtain the different selection coefficients is summarized in Table 1. The letter P stands for probability and NS for not statistically significant.
Flr stands for open flowers per raceme and Chr for flower chroma. Github 1 to run distributional selection differential analyses on our data following Henshaw and Zemel, 2017. The opportunity for selection was 1.74 for overall fitness, calculated using total seed set per plant. When RF was based on the proportion of flower visits, the opportunity for selection was 1.24 for bumble bee, 2.53 for honey bee and 46.48 for leafcutting bee. The average seed set attributable to each bee species, based on proportion of flower visits, was 396.27 seeds per plant for bumble bee; 261.03 for honey bee; and 61.68 for leafcutting bee. When RF was based on the proportion of tripped flowers, the opportunity for selection was 1.15 for bumble bee, 3.16 for honey bee and 35.54 for leafcutting bee. The average seed set of a plant attributable to each bee species was 453.32 seeds for bumble bee, 187.24 seeds for honey bee, and 84.25 seeds for leafcutting bee.

All Bees Combined
Over all pollinators combined, we observed a positive directional selection differential S and selection gradient β on the number of stems per plant indicating selection to increase the number of stems per plant ( Table 2). For the number of racemes per stem, there was a statistically significant positive directional selection differential S and selection gradient β indicative of selection for an increase in the number of racemes. However, we also detected a statistically significant negative quadratic selection differential C 22 although the quadratic selection gradient γ ii was not statistically significant, suggesting indirect non-linear selection on the number of racemes per stem (Table 2)

Bumble Bee
For bumble bees, the positive directional selection differential S and gradient β were both statistically significant for the number of racemes per stem and for the number of stems per plant suggesting selection to increase both traits ( Table 3). In addition, we observed a statistically significant positive correlational selection gradient between the number of open flowers on a raceme and the darkness of a flower (chroma) (γ FlrChr ) for all cases except when RF was based on the number of tripped flowers and the statistical significance of selection coefficients were tested using bootstrapping (Table 3). Bumble bees favored plants with more open flowers per raceme and with darker flowers ( Figure  1B). There was also a positive correlational selection gradient between the number of open flowers and flower reflectivity (γ FlrRef ) but it was only statistically significant when RF was based on the proportion of flowers visited and when log transformed RF regression model was used to detect the significance of the selection coefficients (Table 3).

Honey Bee
For honey bee, there was a statistically significant negative quadratic selection differential (C 22 ) and quadratic selection  The statistical significance of the selection coefficients was determined using either regression models with log transformed relative fitness or using bootstrapping. The methodology followed to obtain the different selection coefficients is summarized in The statistical significance of the selection coefficients was determined using either regression models with log transformed relative fitness or using bootstrapping. The methodology used to obtain the different selection coefficients is summarized in Table 1. The letter P stands for probability and NS for not statistically significant. Rcps stands for racemes per stem and stpp stands for stems per plant. gradient (γ ii) for the number of racemes per stem indicative of non-linear selection ( Table 4). Results of the spline regression analysis indicates that honey bees exert some stabilizing selection on the number of racemes per stem (Figure 2A). For the number of stems per plant, both the positive directional selection differential S and selection gradients β were statistically significant but we also detected non-linear positive selection, suggestive of disruptive selection, with statistically significant quadratic selection differential (C 22 ) and gradient (γ ii ) but only when bootstrapping was used to determine the statistical significance of the selection coefficients ( Table 4). The coefficient of directional selection was much larger than the non-linear selection coefficient (Table 4) which translated into mostly directional selection for increased number of stems per plant as indicated by the spline regression analysis (Figure 2B). Patterns were similar whether RF was based on the proportion of visited or tripped flowers ( Table 4).

Leafcutting Bee
For leafcutting bee, only bootstrapping was used to determine the statistical significance of the selection coefficients. For the number of racemes per stem, we detected both directional and non-linear selection ( Table 5). There was a statistically significant positive selection differential S and gradient β but also a statistically significant negative quadratic selection differential C 22 and gradient γ ii at least when RF was based on the proportion of visited flowers ( Table 5). The spline analysis indicated that leafcutting bees exerted some stabilizing selection on the number of racemes per plant ( Figure 3A). For the number of stems per plant, we detected a positive directional selection differential S and gradient β favoring plants with more stems. Finally, we observed a statistically significant negative quadratic selection gradient γ ii for flower color or hue, indicating some stabilizing selection on hue by leafcutting bees (Table 5 and Figure 3B).

Distributional Selection Differential
When performing distributional selection differential (DSD) analyses, we detected positive directional selection for the number of racemes per stem and the number of stems per plant for all bees combined and for each bee species (Table 6).
We did not detect non-linear selection on any components of floral display size and did not detect selection on any components of flower color for either all bees combined or any of the bee species (Table 6). We present the DSD results to contrast with the results obtained using the Lande and Arnold (1983) approach. We will leave other studies to discuss discrepancies between approaches and below concentrate on the results obtained using the more traditional method originally proposed by Lande and Arnold (1983).

Selection on Flower Color Relative to Floral Display Size
The number of stems per plant and the number of racemes per stem were selected by all three bee species. In contrast,  (Sahli and Conner, 2011). The plants used in this study exhibited a high level of phenotypic variation in both flower color and floral display size and the variation in flower color was greater than typically occurs in wild M. sativa populations. We also observed strong opportunity for selection overall and within each pollinator species. The foraging behavior of pollinators may help explain the difference between selection on plant-level and flower-level attractants. Pollinators forage for rewards and their goal is to collect pollen and nectar to provide for their young and feed themselves. Components of floral display size such as the number of racemes per stem and the number of stems per plant are both indicative of the amount of resources available on a plant. Bumble bees can determine whether a flower offers pollen or not and can detect the number of pollen-producing flowers on a plant. They are attracted to inflorescences, a plant-level attractant, based on the number of pollen-producing flowers . Similarly, bumble bees may be able to detect the number of nectar-producing flowers on inflorescences (Makino and Sakai, 2007). Bumble bees, on the other hand, cannot distinguish between flowers presenting distinct amount of pollen, a flowerlevel attractant, unless it is linked to another trait such as flower size or flower color Thairu and Brunet, 2015). Similarly, while bees have innate preferences for flower color (Simonds and Plowright, 2004;Raine and Chittka, 2007), they learn to associate a flower color with a reward and can switch their preference of flower color for the color providing the most reward (Ings et al., 2009;Thairu and Brunet, 2015). The fact that plant-level attractants such as floral display size directly advertise resource availability to pollinators may help explain why, relative to flower-level attractants, they are more likely to be selected by pollinators within plant populations.
An association between a reward and flower color is more likely to occur among plant populations or plant species of distinct colors rather than within a population where the association between a color and a reward can be broken down by recombination . This may help explain why flower color polymorphisms are more common among than within plant populations (Narbona et al., 2018). Pollinators have been suggested as the selective agents responsible for flower color polymorphisms among populations (Streisfeld and Kohn, 2007) and in some cases the genes responsible for the change in flower color associated with each pollinator have been elucidated (Streisfeld et al., 2013). Similarly, the genetic basis of flower color differences has been elucidated for some plant species and shown to be responsible for the pollinator preference (Bradshaw and Schemske, 2003;Hoballah et al., 2007). However, it remains unclear whether the pollinator preference created the flower color diversification or whether the association between flower color and pollinator arose following the fixation of the flower color in the population or species due to a different factor.
Within plant populations, correlational selection between flower color and floral display size may facilitate the evolution of flower color via pollinators. Correlational selection, as was observed in the current study between number of open flowers per raceme and flower chroma, provides a mechanism to associate a flower-level attractant like flower color to a plantlevel attractant that advertise resource availability to a pollinator. Moreover, correlational selection leads to the development of genetic correlations between traits (Roff and Fairbairn, 2012) and The statistical significance of the selection coefficients was only determined using bootstrapping due to the high number of plants with no leafcutting bee visits. The methodology used to obtain the different selection coefficients is summarized in Table 1. The abbreviation NS stands for not statistically significant. Rcps stands for racemes per stem and hue for flower hue. correlational selection between a flower color and floral display size has been shown to increase the frequency of a color morph within a population even in the absence of differences between color morphs in seedling germination or survival (Gomez, 2000). Correlational selection by pollinators, between flower color and a plant-level attractant, may facilitate the maintenance of flower color polymorphism within plant populations. The role or correlational selection in the evolution of flower color in plant populations deserves more attention.

Distinct Pollinators and Selection of Floral Traits
Distinct pollinators can exert different or conflicting selection on floral traits (Galen, 1989;Sahli and Conner, 2011;Kulbaba and Worley, 2013) and in the current study, we found different patterns of selection on some floral traits by the distinct bee  Relative fitness was based on the proportion of floral visits. Total selection is measured by DSD while dD represents the directional and dN the non-directional component of selection. The general selection differential is S and β is the selection gradient. Statistically significant values are bolded.
traits. Overall, there is directional selection on the number of racemes per stem and stems per plant with indirect nonlinear selection on racemes per stem. There is also correlational selection between number of open flowers per raceme and flower chroma. Clearly bumble bees are solely responsible for the correlational selection while all three bee species exert directional selection on the number of stems per plant. While both honey bees and leafcutting bees exert some stabilizing selection on the number of racemes per stem, the overall selection is mostly directional. Bumble bees were the most abundant pollinators and better trippers than honey bees, the second most frequent visitors. The differential influence of the three bee species on floral traits indicates that the overall pattern of selection in a population will vary with the abundance and efficiency of its pollinators. We therefore, expect temporal or spatial variation in pollinators (Brunet, 2009;Narbona et al., 2018) to influence the temporal or spatial pattern of selection on floral traits (Kelly, 1992;Siepielski et al., 2009Siepielski et al., , 2013Narbona et al., 2018). However, environmental factors may also vary among populations or temporally within populations and affect floral trait evolution (Schemske and Bierzychudek, 2001;Strauss and Whittall, 2006;Caruso et al., 2019). Interestingly, yearly variation in abiotic factors can modify the pattern of correlational selection (Maad, 2000). Both pollinators and abiotic factors should be considered when examining phenotypic selection of floral traits over time or space (Narbona et al., 2018;Sletvold, 2019).

RF and Phenotypic Selection Within Bee Species
The methodology introduced in this study permits evaluation of phenotypic selection by distinct pollinators simultaneously, using the same set of plants. It more realistically describes the process of pollinator-mediated selection in natural populations. Sample sizes remain the same, over all bees and within each bee species, although the proportion of flowers not visited by a bee species may vary among bee species. Of interest is the fact that pattern of selection obtained over combined pollinators could be explained by the patterns observed for each bee species. Moreover, some selection patterns were only significant at the level of a bee species. For example, selection on flower color by leafcutting bees was not expressed at the whole plant level, likely because leafcutting bees were not as common in the study and were responsible for a lesser proportion of the seeds produced by the plants in the population.
We assigned the selection patterns observed in this study to pollinators rather than to another biotic or to abiotic factors. This approach was followed because the number of pollinator visits increase seed set in this plant species (Bauer et al., 2017); M. sativa plants set few seeds in the absence of pollinators (Bohart, 1957); plants were grown in a common environment minimizing variation in resource availability; and herbivory was not observed. Gathering pollination data in phenotypic selection studies will provide useful information on pollinator-mediated selection by distinct pollinators. We will further argue in a separate manuscript that comparing selection gradients between hand-pollinated and open-pollinated plants may not be the most efficient method to assign selection to pollinators (Brunet, in preparation).
The approach introduced in this study relies on good quality pollinator data and a link between visitation and seed set. The pollinator visitation data should be representative of the plant species under study over its flowering season. The plants used to collect pollinator data should represent the variation in floral display size that occurs spatially and temporally in the population. If pollinator types vary throughout the day or the flowering season, one should sample to reflect such variation. To link floral visits to seed set, it is best to sample seeds on visited plants after a period that reflects the time it takes for seeds to reach maturity. Finally, while applied to female reproductive success, the methodology could be extended to male reproductive success. In this case the proportion of floral visits to plant(i) is used for proportional visits by the distinct pollinators as it reflects the pollen leaving plant(i). The total seeds assignable to plant(i), on plant (i) if selfing occurs and on other plants in the population, represent the seed set for male function for plant(i). Results of this study illustrate how the approach proposed can attribute overall phenotypic selection patterns to individual pollinators and we therefore advocate the approach introduced here for future studies examining the impact of distinct pollinators on floral trait evolution.

CONCLUSION
The methodology introduced to isolate and combine the phenotypic selection patterns of distinct bee species on floral traits provides patterns of selection similar to what has been observed in previous studies. The selection patterns observed over all bees could be assigned to specific bee species. All three bee species selected for components of floral display sizes but not all bees favored components of flower color although the selection coefficients were strong. This difference between plantlevel and flower-level attractants could be explained by the fact that floral display size but not flower color directly advertises resource availability to pollinators. Spatial and temporal variation in the abundance of the distinct pollinators is expected to affect patterns of selection of flower traits, particularly for traits differentially selected by the distinct pollinators. Correlational selection between floral display size, a plant-level attractant, and flower color, a flower-level attractant, is expected to facilitate the evolution of flower color by pollinators within plant populations. Studies of pollinator-mediated selection would benefit from combining data on pollinator visitation rates together with seed set and measurements of floral traits when examining the impact of distinct pollinators on floral trait evolution.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JB conceived the study and wrote the manuscript. AB collected the data and samples from the field, processed the samples, and did preliminary data analyses. AF performed the data analyses for the manuscript, prepared the Figures, and provided comments. All authors contributed to the article and approved the submitted version.