Pollinator‐mediated plant coexistence requires high levels of pollinator specialization

Abstract In pollen‐limited plant communities, the foraging behavior of pollinators might mediate coexistence and competitive exclusion of plant species by determining which plants receive conspecific pollen. A key question is whether realistic pollinator foraging behavior promotes coexistence or exclusion of plant species. We use a simulation model to understand how pollinator foraging behavior impacts the coexistence dynamics of pollen‐limited plants. To determine whether pollinators are likely to provide a biologically important coexistence mechanism, we compare our results to bee foraging data from the literature and from a novel experimental analysis. Model results indicate that strong specialization at the level of individual foraging paths is required to promote coexistence. However, few empirical studies have robustly quantified within‐bout specialization. Species‐level data suggest that foraging behavior is sufficient to permit pollinator‐mediated coexistence in species‐poor plant communities and possibly in diverse communities where congeneric plants co‐occur. Our experiments using bumblebees show that individual‐level specialization does exist, but not at levels sufficient to substantially impact coexistence dynamics. The literature on specialization within natural foraging paths suffers from key limitations, but overall suggests that pollinator‐mediated coexistence should be rare in diverse plant communities.


| INTRODUC TI ON
Two properties of pollinator biology pose challenges to plant species coexistence, particularly the persistence of rare plant species in species-rich communities. First, pollen is transferred quickly, so successful plant reproduction depends on sequential or nearly sequential visits to conspecific plants (Karron et al., 2009;Richards et al., 2009), that is, flower constancy (Waser, 1986). The sequence of plants visited will depend on the spatial arrangement of plants and pollinator preferences for the plant species in the local environment Saleh & Chittka, 2007). In the context of dietary generalist pollinators, the frequency of pollination events to a particular plant species does not scale linearly with plant relative abundance. Instead, if pollinator visits are independent, then the probability that a pair of visits will be sequential to the same plant species is the square of the plant species' relative abundance.
Thus, pollination by generalist pollinators should tend to exclude rare plants from the community via a positive feedback loop.
Second, foraging pollinators (i.e., bees) tend to preferentially visit common plant species (Augspurger, 1980;Bruninga-Socolar et al., 2016;Bruninga-Socolar & Branam, 2022;Dauber et al., 2010;Ghazoul, 2006;Kunin, 1997a;Levin & Anderson, 1970;Moeller, 2004;Podolsky, 1992). This tendency coupled with the dependence of plant reproduction on sequential visits by bees suggests that successful pollination of rare plant species by generalist pollinators should be infrequent, and pollen limitation should exclude rare plant species from the community (Ferrière et al., 2004). Yet diverse communities of bee-pollinated plants are ubiquitous in nature, and some of these communities are pollen-limited at least some of the time (Ashman et al., 2004;Knight et al., 2005). How do rare plants find interaction partners and persist in such communities?
For rare species in pollen-limited communities, adequate pollination requires some form of preferential visitation within foraging paths (Benadi & Gegear, 2018). This preferential visitation might take the form of specialized interactions with pollinator species (Benadi, 2015;Valdovinos et al., 2013), spatial clumping among conspecific plants , or specialization within generalist pollinator species at the level of individuals or even of single foraging paths (Brosi & Briggs, 2013). The potential contribution of individual-level foraging specialization within generalist pollinator species has received little attention.
Generalist pollinator species may have an innate, or average, preference for some plant species over others (e.g., due to morphological adaptations such as tongue length (Harder, 1985) or balancing nutritional requirements (Williams & Tepedino, 2003)), but will nevertheless visit additional species depending on the foraging context . Generalist pollinator species may exhibit adaptive foraging in response to their environment, despite a general preference for common plant species, and therefore a segment of the pollinator community might visit rare plant species enough to reliably pollinate them (Benadi, 2015;Benadi et al., 2012;Bruninga-Socolar et al., 2016Goulson, 1994;Kunin & Iwasa, 1996;Revilla & Křivan, 2016;Valdovinos et al., 2013;Waser, 1978). Recent theoretical work on the role of pollination in plant species coexistence explores variation in pollinator species preferences and finds that the ability of generalist pollinator species to switch plant interaction partners and to specialize on different plant species promotes plant species coexistence (Benadi, 2015;Revilla & Křivan, 2016;Valdovinos et al., 2013).
For generalist pollinators, sequential visits to rare plant species may also result from variation between conspecific individuals in short-term foraging behavior, that is, within individual foraging paths. Existing models of plant-pollinator interactions and their effects on plant species coexistence assume that all pollinators or all pollinators of the same species forage identically (e.g., Benadi, 2015;Benadi et al., 2012Benadi et al., , 2013Essenberg, 2012;Levin & Anderson, 1970).
Yet Heinrich (1976) qualitatively demonstrated that individual-level specialization exists in bumblebees and recent work has shown small-scale differences among individuals in the sequence of plants visited in foraging paths in the same experimental array (Saleh & Chittka, 2007). Individual-level variation in foraging behavior should be especially important for species coexistence when a small number of species dominates the pollinator community (i.e., the species abundance distribution ubiquitous in ecological communities;McGill et al., 2007, Song & Feldman, 2014, Winfree et al., 2015. However, to our knowledge, very few studies have quantified variation among individual bees at the level of their foraging paths, where sequential visits to conspecific individuals of a rare plant species may occur. Song and Feldman (2014) use coupled plant-pollinator population dynamics models to show that adaptive foraging among individuals of a single pollinator species promotes plant species coexistence.
However, plant competition for resources other than pollinators has a stronger effect on plant coexistence than competition for pollinators in that study. Therefore, these results might fail to capture the dynamics of pollen-limited systems where pollen limitation is the mechanism limiting plant population sizes. Further, while previous work shows interesting avenues for plant-pollinator interactions to influence plant species coexistence, no study has modeled variation in the stepwise foraging decisions of individual bees (i.e., intraspecific variation among foraging paths), a key puzzle piece in understanding the dynamics of plant coexistence in pollen-limited systems because of the requirement that pollen be delivered in sequential or nearly sequential visits for successful pollen transfer.
In this study, we focus on the potential contribution of individuallevel specialization to coexistence dynamics in a pollen-limited twoplant system. First, we use observed foraging paths of Bombus terrestris in experimental flower arrays to estimate the level of individual-level specialization in a model generalist pollinator species. We then develop a simple, computationally efficient simulation model for the coexistence dynamics of two pollen-limited plants competing for shared pollinators. We use this simulation to explore coexistence dynamics across a region of parameter space that spans variation in the average preferences of two pollinator species, variation in the magnitude of individual-level variance around those averages, and variation in the starting abundances of the two plants. We ask whether and when realistic levels of individual variation in path-level foraging preferences meaningfully alter coexistence dynamics. Finally, we review the sparse existing literature on the degree of individual specialization in nature and assess whether observed bee preferences are sufficient to support pollinator-mediated plant coexistence.

| Experimental data on individual bee foraging
Because individual variation in foraging paths is rarely measured empirically, we present experimental data that allows us to quantify variation in individual bee foraging behavior. We observed foraging bumblebees (B. terrestris) in a laboratory setting, using a repeated measures design to test for statistical differences in individual foraging behavior. Naïve, marked B. terrestris workers foraged in an artificial meadow with 42 artificial flowers consisting of blue discs with two types of floral guides: 21 flowers with yellow floral guides containing pollen rewards and 21 flowers with orange floral guides containing nectar rewards (Appendix 1: Figure A1). These two artificial flower types simulate two plant species from which bees get different rewards. We tested 41 B. terrestris workers in 199 individual trials across 20 combinations of pollen and nectar quantity and quality (Appendix 1: Table A1). Each worker was not tested for all combinations of rewards. Each worker completed a mean of 4.85 ± 0.72 trials, and a maximum of 16 trials. For each trial, we recorded the number of flowers visited of each type. To be counted as a visit, the bee had to directly contact and collect the floral reward. After each trial, we calculated the amount of reward collected by measuring the sugar solution left in the nectar flowers and weighing the pollen collected by the bee. See Konzmann and Lunau (2014) for additional details.

| Data analysis
To assess whether individual bees differed in their foraging preferences, we modeled individual bee foraging decisions (probability of visiting pollen vs. nectar flowers) as a function of the quality and quantity of pollen and nectar rewards as well as a random effect of individual bee. If all bees are identical in preference, then the random effect variance should be indistinguishable from zero and including the random effect should not substantially improve the model fit. Therefore, we sought direct inference on the random effect variance to examine whether the data rule out near-zero variance, where random effect variance indistinguishable from zero would indicate no differences between individual bees. Because frequentist model fits may underestimate the uncertainty in variance parameters (Kéry, 2010), we fit the model under a Bayesian mode of inference using vague priors and five chains, each with 1000 iterations of burn-in and 50,000 iterations of sampling. The Gelman-Rubin diagnostic was 1.0 for all margins, including the multivariate Rhat ("coda" package; Plummer et al., 2006). We compared prior and posterior distributions for the standard deviation to assess whether the data constrain the standard deviation away from zero.
To confirm whether the random effect of individual variation should be in the best fit model, we compared the models with and without individual variation using both a frequentist likelihood ratio test and a Bayesian indicator variable analysis (Appendix 1). All analyses were done in the program R (R Core Team, 2017) using JAGS (Plummer, 2003) via the R package "rjags" (Plummer, 2016) and are publicly available .

| Simulation model
We wish to investigate the dynamics of a plant-pollinator community with two pollen-limited annual plant species that are obligate on two shared pollinator species. To do so, we use a stochastic difference Our modeling framework simulates the number of plants of each species in each successive generation based on the number of successful pollen transfers and the probability that a fertilized ovule recruits to the adult stage in the next generation. For this study, we treat the latter probability as fixed because we are most interested in exploring parameters related to bee foraging. The model allows the number of successful pollen transfers to depend on bee foraging preferences, plant abundances, and the decay of pollen transfer after bees visit intervening flowers (Richards et al., 2009). We assume that pollinator populations are not limited by plant populations (e.g., they are nest site limited), such that the number of flower visits per year is constant. This assumption allows us to model the effects of a particular pollination regime on plant coexistence without modeling changes in the bee community over time. The extent to which bees are limited by non-food resources such as nest sites is understudied (Harmon-Threatt, 2020; Roulston & Goodell, 2011). However, a few studies document a positive relationship between nest site availability and bee abundance and species richness across nesting types (Potts et al., 2003(Potts et al., , 2005 and demonstrate that neither floral nor nesting variables alone can explain bee population and community dynamics (Potts et al., 2005). For analytical simplicity and inferential clarity about pollinator-mediated plant coexistence, we conceptualize our system as a nest site-limited bee community. The model does not include a spatial component or other mechanisms of plant coexistence, such as competition for non-pollinator resources.
We further assume that the identity of the plant in the ith visit delivered by the jth bee is independent of the plant identity on the previous visit and is conditioned on the bee's foraging preferences, described below. Thus, our simulation computes the probability that any single bee visit will deliver conspecific pollen, samples binomially from this probability and the total number of visits delivered by a given bee, and then sums across all bees in the community. Our full model is given by where P nt is the population size of plant species n in year t, S nt is the number of successful pollen transfers to plant species n in year t, and ρ n is the probability (treated as constant) that a fertilized ovule of plant species n recruits to an adult in the next generation (Table 1).
S nt is obtained by summing V n , the per-bee number of successful pollen transfers to plant n, across all bee individuals j of all bee species i. To calculate V nij for a given bee individual, we draw from a binomial distribution with number of trials v ij equal to the total number of floral visits performed by the individual bee (treated as constant; Table 1). The per-visit probability of successful pollen transfer is Ranges of parameter values in the simulation model.
Bee abundance, A i 100 per bee species (Appendix 1) Both elevated: 150 bees (Appendix 1: Figure  Pollen transfer probabilities, τ k Sequentially for each visit in a bee foraging path: 1, 1, 1, 1, 0.7, 0.5, 0.2 (Richards et al., 2009) A vector whose 1st value is the probability of conspecific pollen transfer from immediately sequential visits, the 2nd value is the probability given one intervening visit to another plant species, and so on.
Both elevated: 30,000 visits ( Figure A12) Both lowered: 20,000 visits ( Figure A13) Note: Bee species 2 is constrained to always have a positive mean preference because the system is symmetric with respect to bee species mean preferences; that is, the scenario where the bee species have opposite plant preferences (one has a positive mean and one has a negative mean) is captured by allowing bee species 1 to vary to negative values, and the situation where both species have negative mean preferences is the same as the situation where they both have positive mean preferences, except the identity of the mutually preferred plant species is switched. Parameter values are also provided for 12 model runs assessing the sensitivity of our results to the values related to the bee and plant populations.
equal to the probability that sequential visits result in conspecific pollen transfer τ k multiplied by the probability that visits are sequential. Note that τ k is a vector whose first value is the probability of pollen transfer from immediately sequential visits, the second value is the probability of pollen transfer given one intervening visit to another plant species, and so on ( Table 1). The values of τ k are taken from Richards et al. (2009) who model variation in pollen deposition by individual pollinators. We assume that a visit that simultaneously delivers pollen from multiple previous visits to conspecific plants results in more total reproduction than a bee visit delivering pollen from fewer conspecifics.
The probability of sequential visits to plant n is a function of the foraging preference α ij of the individual bee (constrained to take values greater than 0), the density of the plant species in year t, and a parameter β that controls whether the bee preferentially forages on common (β > 1) or rare (β < 1) plants. Generalist bees may forage preferentially on rare plant species to balance nutritional requirements, for example (Cook et al., 2003, Hendriksma & Shafir, 2016 see Section 4). α ij and β thus modify the true density of each plant species to an effective density determined by each bee's species-and individual-level preferences (α and β, respectively), and the probability that a given visit-pair involves plant species n on both visits is the square of this effective density. Plant species 1 is preferred when α ij is greater than 1, and plant species 2 is preferred when α ij is less than 1. Each bee's α is sampled from a log-normal distribution specified by its species mean preference and standard deviation, which represents the intraspecific variation in preference around the mean (Table 1): where μ i is the mean foraging preference of the ith bee species and σ i is the standard deviation of the ith bee species. The mean and standard deviation are chosen a priori to span a wide range of possible values (Table 1). Values of β are systematically varied in the model to explore effects of preference for rarity on bee foraging and plant coexistence (Table 1).
All simulations were run in program R (R Core Team, 2017) for 100 plant generations and the R script is publicly available . In the infinite time limit, the only possible outcome of our simulation is eventual extinction of both plant species, but we found that all simulations rapidly settled into a welldefined basin of attraction in many fewer than 100 iterations, and moreover that in all cases (for all choices of parameter values and starting conditions), 100 replicate simulations from identical starting conditions universally found the same basin of attraction and remained there after 100 iterations. We ran simulations to explore the range of parameter space of the bee species means, standard deviations, and preference for rarity to determine how bee preferences affect plant coexistence (Table 1). Parameters unrelated to bee foraging were held constant in most simulations (

| Coexistence in our simulation context
Traditionally, ecological coexistence is analyzed based on the invasibility criterion: Can a population increase when rare (Chesson, 2000)?
In pollinator-mediated coexistence, this criterion is difficult to meet because sequential visitation rates to a rare plant species should approach zero in the limit of low relative abundance. Instead, we hypothesized that coexistence will often be stable within a limited basin of attraction at intermediate relative abundance. Therefore, we measure the strength of coexistence as the fraction of communities with both plant species persisting after 100 generations. For each combination of parameter values, we run 100 simulations to calculate this fraction. Thus, our criterion for coexistence is not invasibility, but rather the requirement that neither plant is likely to leave the basin of attraction and go extinct over our 100-generation time window (Caswell, 1978;Socolar & Washburne, 2015;Valdovinos et al., 2013).

| Experimental data on individual bee foraging
The posterior distribution of the standard deviation of the random effect was constrained away from zero (Figure 1), indicating that variation in the random effect, attributable to variation among individual bees, is included in the best model. The standard deviation of the random effect is mathematically equivalent to the standard deviation of bee species preferences in our simulation model. In our data analysis, the standard deviation (sigma or σ) takes a value close to 1. In our simulation model, increasing the standard deviation from 0 to 1 in agreement with our experimental result had a very minor effect on plant species coexistence when there was no preference for rare plants (Figure 2c), but improved coexistence slightly when there was moderate preference for rare plants ( Figure 2b). When there is no interspecific variation (bee species' means are equal), intraspecific variation only improves coexistence when the standard deviation of at least one species is 3, a level much higher than our experimental result ( Figure 2). The likelihood ratio test and Bayesian indicator variable analysis confirmed that individual variation is included in the best fit model (Appendix 1).

| Effects of variation in bee foraging preferences on plant coexistence
Our simulation model of individual bee foragers requires strong bee specialization for plant coexistence. Bee specialization in the model is a result of density-dependent preference for rare plants and/or strong variation in bee preference, such that different individuals strongly prefer different plant species. In our model, this variation could arise from differential species-specific mean preferences (μ i ), high levels of individual variation around the means (σ i ), or differences among individual bees in their preference for common versus rare plants (β).  Figure 2). Strongly opposing differences in mean preference always promote plant coexistence (Figure 2, left half of all plots), even when density-dependent preference for rare plants is weak (Figure 2b) or non-existent (Figure 2c,d), that is, pollinators prefer to visit common plants. When initial plant population sizes vary such that one plant species starts out as rare, stronger bee specialization (greater difference between μ 1 and μ 2 ) is required for coexistence in our model ( Figure 3). With a strongly skewed ratio of initial plant population sizes (e.g., 625:1) extremely high levels of specialization are required for coexistence (Figure 3). When the initial plant population sizes are equal, smaller differences between the bee species' means maintain coexistence ( Figure 3).
In our model, high standard deviation around the mean bee spe- bee's preference is drawn from a log-normal distribution with a mean of 1.3 or −1.3, depending on its species identity, and no standard deviation (Equation (7)). The resulting value from this distribution (3.67 from a mean of 1.3; 0.27 from a mean of −1.3) is α in Equation (3). For the purposes of this example, we assume that β = 1 (no preference for rarity) and the effective density of plant species 1 thus becomes (α*P 1t )/((α*P 1t ) + P 2t ) (Equation (3)). For a bee species mean of 1.3, the effective density of plant species 1 is 0.79, an increase from its "true" density of 0. individual pollinator specialization within foraging paths for plant coexistence; that is, the vast majority of visits within a single foraging trip must be to a single plant species.
For the case where the initial population sizes of the plant species are unequal, for example with a ratio of 25:1, the bee species means must differ in magnitude by 4.4 for coexistence to occur, that is, the mean preference of bee species 1 is 2.2 and the mean preference of bee species 2 is −2.2 ( Figure 3). Repeating our calculation above assuming that β = 1 and there is no intraspecific variation in bee foraging behavior, the α values for bee species 1 and bee species 2 are 9.03 and 0.11, respectively (Equation (7)). Using Equation (3) . When the bee species prefer different plant species, as indicated by one positive mean and one negative mean, plant species coexistence is supported (left side of all plots). When the bee species' means are equal, there is no interspecific difference in preference. High standard deviation (high within-species individual variation) around both bee species means increase coexistence. The greatest coexistence is obtained when there are high levels of variation around both bee species' means (high standard deviation; upper right plot of each panel has the greatest extent of orange compared to other plots within the same panel). When only one bee species' standard deviation is high, coexistence occurs less than in the previous case but more than when variation is low for both species (compare top row or right column of plots to bottom left plot of all panels). When bees prefer common plants ( Figure 3), which corresponds to delivering over 90% of visits to one plant when the plants are in equal mixture.
Our experimental results using bumblebees suggest that individual-level specialization exists in foraging pollinators in the lab (i.e., intraspecific variation; Figure 1). However, the level of observed variation (SD = 1) does not appreciably affect plant coexistence in our model when there is no preference for rare plants (Figure 2c), but did improve coexistence slightly when there was moderate preference for rare plants (Figure 2b). Other studies have analyzed floral choice among bumblebee and honeybee workers, but did not calculate means of individual behavior for quantitative comparison (Grüter et al., 2011;Heinrich, 1976Heinrich, , 1979 (Thomson, 1981;Waser, 1986).
Very little previous work has investigated effects of specialization of individual foragers of a generalist pollinator species on plant species coexistence. Song and Feldman (2014) found that some pollinators specialize on the plant species that is rarer in a two-species community only when its local density is high, which they suggest is more likely when total plant density in the community is high. We add that spatial clumping of conspecifics within a plant community could have a similar effect as bees forage optimally on plants that are close together. A recent empirical study suggests that bee responses to conspecific spatial clumping of plants in a multi-species community may drive their stepwise foraging choices within foraging paths, affecting the delivery of sequential visits to the same plant species and resulting in persistence of rare species in the community (Bruninga- . Most studies of the effects of pollinator visitation on plant species coexistence are conducted in a network context (e.g., Valdovinos et al., 2013) or using differential equations that model populations over time without modeling foraging paths (e.g., Revilla & Křivan, 2016). Both approaches elide the short-term dynamics of effects of bee foraging behavior on pollen delivery to plants.
A key question is whether the levels of pollinator specialization that drive pollinator-mediated coexistence in our model are widespread in nature. We consider four mechanisms that might deliver sufficiently specialized within-bout foraging dynamics: specieslevel specialization, individual-or bout-level specialization, densitydependent preference for rare plants, and spatial or temporal clumping of flowers.

| Species-level specialization
Most bees are not dietary specialists, and very few bee species are monolectic (specialized on the pollen of a single plant species) (Michener, 2000). However, oligolecty (specialization on a plant genus or family) is not uncommon (Michener, 2000). For example, 43% of bee species in the tribe Anthidiini are oligolectic at the level of plant family, subfamily, or tribe (Müller, 1996), and 30% of bee species in a region of subtropical Brazil are oligolectic (Schlindwein, 1998). Thus, in species-poor systems with a single plant per family, oligolectic bees might provide a powerful mechanism for plant species coexistence if they behave as monolectic.

F I G U R E 3
With no preference for rare versus common plants (β = 1) and no intraspecific bee variation (σ = 0), the two plant species coexist (orange regions of graph) with large differences between the mean bee species' preferences. Coexistence is more difficult as the ratio of the initial plant population sizes increases, that is, if one plant species starts out as rare compared to the second plant species. When the initial plant population sizes are equal, coexistence is obtained at smaller differences between the mean bee species' preferences.
Note that competition between two species of congeneric plants for oligolectic pollinators is precisely analogous to the two-plant scenario for generalist pollinators in our model.
Documented cases of strongly monolectic pollinators are rare and often involve unusual examples of coevolution, such as the classic example of long-tongued moths and long-spurred orchids (e.g., Netz & Renner, 2017). In such specialized cases where bees and plants have a strong reciprocal preference, pollinators might easily mediate persistence of the plant in arbitrarily species-rich systems, and the selective forces that guided the evolution of such elaborate signaling seem likely to involve pollen limitation (Kiester et al., 1984).
However, we expect that these extraordinary cases account for the persistence of only a small minority of plant species.

| Individual-or bout-level specialization
Evidence of high individual bee specialization in nature is limited.
While multiple studies report the frequency of visits to a given plant species within a foraging bout, these studies generally are not accompanied by data on the relative abundance of that plant within the community. Thus, it is possible (and in our view likely) that reported cases of high apparent specialization simply reflect preferential visitation to common plant species (or the commonest plant among the genus or family preferred by an oligolectic bee). Nevertheless, we note that many reported foraging paths are entirely restricted to a single plant species, which suggests that pollinator-mediated coexistence might be possible. In an alpine system, 77% of individual bumblebees visited only one plant species within a foraging bout (Brosi & Briggs, 2013), and in an Australian garden, 88% of foraging trips of a stingless bee consisted of visits to only one plant species (White et al., 2001). However, only 35% of bumblebee foraging paths in a German meadow visited exclusively one plant species, and 37% of foraging paths included visits to at least three plant species (Raine & Chittka, 2007). In an agricultural system in Uruguay, approximately 80% of individuals of two bumblebee species collected only one type of pollen in corbicular pollen loads (Rossi et al., 2015).
However, only 60% of individuals carried one type of pollen in nectar expressed from the abdomen, including most of the bees whose corbicular pollen loads were also tested (Rossi et al., 2015). Thus, studies that only examine corbicular pollen loads as a test of constancy may underestimate the diversity of plants visited within foraging paths. Among solitary bees, the species Ceratina australensis shows a high degree of specialization among individual bees as measured by comparing the composition of individual larval pollen provisions to mean population-level pollen composition (Gaiarsa et al., 2022).
Three species in the genus Osmia collected only one plant family in 44%-58% of pollen loads (each pollen load corresponds to one foraging path, subject to the caveat above; Eckhardt et al., 2014).
However, two additional Osmia species showed no specialization within foraging paths at all (Cane, 2011). It is worth noting that nectar is a refillable resource in many plant species whereas pollen is not; therefore, individual bees foraging for pollen may be more likely to forage on rare plant species if the pollen resources offered by their preferred or common plant species have been exhausted.
One way to circumvent the need for data on the relative abundance of plants is to use Thomson's interview method of assessing bee preference, in which individual foragers are experimentally confronted with a choice of flowers of different plant species offered to the bees by the researcher (Thomson, 1981). Two studies use this method to document high specialization of bumblebee foragers on either of two congeneric plant species (Raine & Chittka, 2005;Wilson & Stine, 1996). In a mixed field of white and red clover (Trifolium repens L. and T. pratense L., respectively), Wilson and Stine (1996) show that 68% of Bombus vagans workers chose white clover when interviewed if the previous flower they had visited was also white clover.
88% of workers chose red clover when interviewed if the previous flower they had visited was also red clover (Wilson & Stine, 1996). Raine and Chittka (2005) calculate a bee species-specific constancy index that compares the number of visits to the same plant species as the previous flower visited to the number of visits to a different plant species than the previous flower visited. They calculated constancy indices for three bumblebee species and found that the constancy indices for these species ranged from partial constancy to complete constancy (Raine & Chittka, 2005). These results suggest that bout-level specialization in bumblebees might be sufficient to promote coexistence even of congeneric plants. The 79% sequential visitation rate required by our model sits between the 88% and 68% sequential visitation rates to red and white clover, respectively, documented by Wilson and Stine (1996).

| Density-dependent preference for rare plants
Density-dependent preference for rare plants is under-explored and further empirical work in natural systems is needed. We expect that bee preference for rare plant species may occur in nature, for example, when generalist bees require a certain resource that only a specific plant species can provide, and that plant species happens to be rare in the community (Williams & Tepedino, 2003). Visitation to rare plants may provide pollinators a competitive advantage against other pollinator species (Valdovinos et al., 2013). Rare plants may also offer respite to generalist pollinators from harmful secondary compounds present in common plants, in some systems (Bukovinsky et al., 2017;Eckhardt et al., 2014), or allow generalist pollinators to balance collection of multiple necessary nutrients, for example, essential amino acids (Cook et al., 2003;Hendriksma & Shafir, 2016).

| Spatial or temporal clumping
Because of the high levels of individual bee specialization required by our simulation model for plant coexistence, we predict that in speciesrich, pollen-limited plant communities, rare plant species should occur in clumps of high local abundance. Optimally foraging animals maximize energy gained per unit of time or effort (Schoener, 1971), or per nutrition type (e.g., amino versus fatty acids (Ruedenauer et al., 2021)).
In the case of pollinators, transit between flowers is a key component known to be pollen-limited (e.g., Crone & Lesica, 2006) and recent work highlights that heterogeneity in floral resource availability affects pollinator visitation and plant reproduction (Labonté et al., 2023). More broadly, in many systems, certain plants or plant communities flower during short periods of the year, for example, spring ephemerals (Kudo et al., 2008).

| CON CLUS IONS
The role of pollinators in mediating plant coexistence is of major interest both as a potentially important aspect of modern coexistence theory and for its basic and applied implications for pollen limitation: as a rule, diverse plant communities should not be severely pollenlimited unless their pollinators tend to promote coexistence rather than exclusion. Previous work suggests that realistic levels of individual specialization can play a role in rescuing plant species coexistence when differences in mean species preferences-or overall densitydependent preference for rarity-exists but is weak. Outside this regime, our results suggest that only extremely high levels of individual pollinator specialization are capable of maintaining coexistence in pollen-limited plant communities. Species-level specialization among bees is variable in nature, but may rarely be sufficient to promote coexistence among congeneric co-flowering plants. Individual-and bout-level foraging specialization, coupled with spatial clumping of rare plants, might be sufficient to provide a more general coexistence mechanism. We provide empirical evidence of intraspecific variation among bumblebee workers, but at a level insufficient to contribute strongly to plant coexistence in our simulation model. In the existing literature on bee foraging preference, we find that few studies permit rigorous quantification of bee foraging specialization, which requires quantitative data on the relative abundance of local plant species.
However, a handful of studies using Thomson's interview method suggest that bumblebees might conceivably promote coexistence even among pollen-limited congeners. We conclude that pollinator specialization should be included in models of plant coexistence and propose that future empirical work in pollen-limited plant communities investigate the role of pollinator specialization in the persistence of those communities, particularly the persistence of rare plant species.

ACK N OWLED G M ENTS
We thank Rachael Winfree, Elizabeth Crone, and Tina Harrison for their thoughtful comments on this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R code used in this paper were deposited in the A PPE N D I X 1

Simulation model parameter choice
We selected biologically reasonable parameters based on the study system with which we are most familiar (Bruninga-Socolar et al., 2016, where we estimate that roughly 200 bees may pollinate a patch of flowers and roughly 500,000 flowers may populate a patch. This field system differs from our simulation model system because the most abundant plant species is not annual. We determined the probability that a fertilized ovule recruits to the adult stage in the next generation based on the values we chose for the parameters in Table 1.

ANALYSIS
When comparing a random effects model with one without the random effect, standard AIC-based approaches to model selection are difficult to implement because estimation of the number of effective degrees of freedom associated with the random effect is challenging (Kéry, 2010). Instead, we applied two alternative approaches to assess the strength of evidence for the inclusion of the random effect. First, we fit models with and without the random effect in the R package lme4, and we compared their fits using a likelihood ratio test. The likelihood ratio test between the models with and without the random effect representing variation among individual bees was significant (χ 2 = 125.67, p < .0001).
Second, we performed fully parametric Bayesian model selection using indicator variables (Hooten & Hobbs, 2015)  Note: For example, in each of five tests, bees were presented with a different value for nectar quantity while nectar quality, pollen quantity, and pollen quality were held constant. In the next series of five tests, bees were presented with a different value for nectar quality while nectar quantity, pollen quantity, and pollen quality were held constant, and so on for all 20 tests. See Konzmann and Lunau (2014) for additional details.

TA B L E A 1
Combinations of rewards used in 20 tests of Bombus terrestris foraging behavior.

F I G U R E A 2
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for bee abundance of each species are both elevated to 150 bees compared to abundances of 100 for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 3
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for bee abundance of each species are both lowered to 50 bees compared to abundances of 100 for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 4
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for bee abundance are asymmetric with species 1 set to 100 individuals and species 2 set to 150 individuals compared to abundances of 100 each for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 5
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for bee abundance are asymmetric with species 1 set to 100 individuals and species 2 set to 50 individuals compared to abundances of 100 each for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 6
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for the probability of conspecific pollen transfer are all lowered (0.7, 0.7, 0.7, 0.7, 0.49, 0.35, 0.14) compared to the set for our main model results (1, 1, 1, 1, 0.7, 0.5, 0.2; Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 7
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the values for the probability of conspecific pollen transfer decay faster (1, 0.9, 0.7, 0.5, 0.3, 0.2, 0.1) than the set for our main model results (1, 1, 1, 1, 0.7, 0.5, 0.2; Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 8
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the probabilities of pollen becoming an adult plant of each species are both elevated to 0.2 compared to the value 0.04 used for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 9
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the probabilities of pollen becoming an adult plant of each species are both lowered to 0.01 compared to the value 0.04 used for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 1 0
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the probabilities of pollen becoming an adult plant of each species are asymmetric with species 1 set to 0.04 and species 2 set to 0.05 compared to 0.04 for both species for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 11
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the probabilities of pollen becoming an adult plant of each species are asymmetric with species 1 set to 0.04 and species 2 set to 0.03 compared to 0.04 for both species for our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.

F I G U R E A 1 2
We ran the model 12 times to test the sensitivity of our results to the parameters listed in Table 1. Here, the lifetime visits per bee individual are set to 30,000 for each bee species compared to 25,000 used in our main model results ( Figure 2). All other parameters in Table 1 stayed the same. All parameters given in Table 1 stayed the same.