A Tale of Two Morphs: Modeling Pollen Transfer, Magic Traits, and Reproductive Isolation in Parapatry

The evolution of the flower is commonly thought to have spurred angiosperm diversification. Similarly, particular floral traits might have promoted diversification within specific angiosperm clades. We hypothesize that traits promoting the precise positional transfer of pollen between flowers might promote diversification. In particular, precise pollen transfer might produce partial reproductive isolation that facilitates adaptive divergence between parapatric populations differing in their reproductive-organ positions. We investigate this hypothesis with an individual-based model of pollen transfer dynamics associated with heterostyly, a floral syndrome that depends on precise pollen transfer. Our model shows that precise pollen transfer can cause sexual selection leading to divergence in reproductive-organ positions between populations served by different pollinators, pleiotropically causing an increase in reproductive isolation through a “magic trait” mechanism. Furthermore, this increased reproductive isolation facilitates adaptive divergence between the populations in an unlinked, ecologically selected trait. In a different pollination scenario, however, precise pollen transfer causes a decrease in adaptive divergence by promoting asymmetric gene flow. Our results highlight the idea that magic traits are not “magic” in isolation; in particular, the effect size of magic traits in speciation depends on the external environment, and also on other traits that modify the strength of the magic trait's influence on non-random mating. Overall, we show that the evolutionary consequences of pollen transfer dynamics can depend strongly on the available pollinator fauna and on the morphological fit between flowers and pollinators. Furthermore, our results illustrate the potential importance of even weak reproductive isolating barriers in facilitating adaptive divergence.

divergence). While the S-locus governs which reproductive organ appears at the high position and which appears at the low position in a given flower, the actual heights of the high and low positions seem likely to be governed by many loci of small effect size unlinked with the S-locus.
Although the details of this are not known for any heterostylous species, Bissell and Diggle [6] did find that anther position and stigma position are genetically independent in species in Nicotiana; they argue that this independence provides useful evolvability to species in adapting to different pollinators. Other closely related heterostylous species are also known to exhibit differences in anther position, stigma position, and anther-stigma separation [e.g., 7], suggesting that the genetic independence found by Bissell and Diggle [6] may be general. We thus represent the heights of the reproductive organ positions as quantitative traits, x and y, with continuous values representing the additive effect of many alleles in an infinitesimal model [8,9]; while this architecture is not known to correspond to any particular heterostylous species, neither is it contradicted by any empirical findings of which we are aware, and it seems a simple and reasonable choice unlikely to bias our results in any important way. Environmental variance is not included in our model; in other words, the heritability of the traits is 1. Empirical evidence suggests that genetic variation and heritability for reproductive-organ positions is often quite high [10][11][12][13][14][15], so this seemed a reasonable choice for simplicity; evolution is likely more rapid in our model than it would be in real heterostylous systems as a result, but this does not affect our conclusions since we do not attempt to draw inferences regarding the absolute time required for adaptation. In any case, the stochasticity of the movement of pollen in our model (see Pollination phase) produces an effect very similar to the expected effect of environmental variance in reproductive-organ positions.
One consequence of our chosen genetic model is that the "polarity" of the S trait (whether S = 0 represents pin or thrum) is emergent: it is possible for the two populations to arrive at different polarities (almost always when the pollinator crossover probability was very low, since gene flow pushes both patches toward a shared polarity). This resulted in almost complete reproductive isolation between patches, because morphs anatomically well-suited to exchange pollen (a pin in one patch and a thrum in the other) were blocked from cross-fertilizing by pollen incompatibility due to their possession of the same S trait value. This phenomenon of polarity differences among closely related heterostylous populations has not been observed in nature, likely due to shared ancestry [16]. For this reason, only realizations in which the two patches arrived at the same polarity were used; this constituted the large majority of realizations, since gene flow tended to cause the polarities in the two patches to synchronize. Removing realizations with opposing polarities between the patches amounts to enforcing a shared genetic architecture for heterostyly among the plants of the two patches, a reasonable assumption given that they represent diverging populations within the same species. Not removing realizations with opposing polarities, on the other hand, would produce incorrect results; our model would show very high degrees of reproductive isolation in those cases, but that isolation would be due to a mechanism that is not observed empirically because the genetic architecture of heterostyly is fixed by common ancestry. By removing these realizations, we focus attention on the empirically justified case, in which reproductive isolation between populations due to precise pollen transfer is the result of smaller, quantitative differences in reproductive organ heights, as shown in Fig.   1B.
The reproductive state of individuals is tracked with several non-genetic (i.e., non-heritable) state variables. Each individual has a number of pollen grains, p, which can be taken up and transported by pollinators, a number of ovules, o, that can be fertilized in each year, and a style clogging index, s, representing the empirical fact that the style of a flower becomes progressively clogged by pollen tubes that decrease the probability of fertilization. The details of these mechanics will be discussed below.

Processes and scheduling
Overlapping yearly generations of individuals are modeled as consisting of three phases executed consecutively in each year: germination, mortality, and pollination. In the germination phase, seedling individuals are generated from all fertilized ovules, and then seedlings die with a probability inversely proportional to their relative fitness in their local patch until the population size of each patch is less than or equal to the carrying capacity of the patch; conceptually, this represents the germination of seeds and the subsequent trait-dependent natural selection of seedlings during maturation. In the mortality phase, individuals die randomly with a fixed probability; this may be taken to represent either truly random mortality or natural selection on traits not modeled. In the pollination phase, the surviving individuals are visited by pollinators that transport pollen between them, resulting in the fertilization of ovules that will germinate at the beginning of the following year.

Interactions
In the germination phase, seedlings interact through competition in the sense that the probability of mortality for each seedling depends upon the number of other seedlings alive; in other words, natural selection during this phase is "soft selection" [17]. In this model, parapatric populations are connected due to "crossover" of pollinators between populations. In crossover events, a pollinator visits a flower in one patch, receiving pollen, and then delivers that pollen to a flower in the other patch (and then returns to its native patch; see Pollination phase for further mechanistic details). These crossover events thus carry pollen between the populations, resulting in gene flow to the extent that that pollen successfully fertilizes flowers in the destination population; this gene flow is the only way in which the two populations interact. The crossover probability, c, ranges from complete allopatry (c = 0.0) to full sympatry (c = 0.5), although even in full sympatry the model separates the flowers into two discrete "patches" with distinct carrying capacity; this case corresponds to a heterogeneous local environment providing two ecologically distinct niches that are sufficiently spatially proximate as to produce no visitation bias in the pollinators. The crossover probability thus produces geographic isolation between populations in a similar manner to the behavioral isolation that would be produced by a pollinator preference for visitation of one type of plant over the other. Indeed, this model is nearly analogous to a fully sympatric model with two types of flowers differentiated by floral traits such as petal color that cause pollinators to exhibit visitation preferences [e.g., 18,19]. However, the distinct carrying capacity of the patches makes sense from a spatial perspective, but would be difficult to justify in a sympatric model of ethological isolation. For this reason, the crossover probability in this model controls geographic isolation, not ethological isolation, and the two framings are not equivalent.

Stochasticity
Stochasticity is present in many aspects of this model. The initial state of individuals is stochastic, such that the particular distribution of trait values slightly affects the speed with which dimorphism develops from the initial unimodal distribution, as well as the "polarity" of the S trait (whether S = 0 represents a pin or a thrum) once stable dimorphism is established.
Demographic stochasticity is present due to finite population size, occasionally resulting in the extinction of one or both populations, and more generally affecting the evolutionary outcome of the model through drift. Stochasticity manifests in many aspects of pollination events (see Pollination phase): the particular flowers visited, the number of pollen grains transported, whether each pollen grain sticks to the pollinator, the precision (or lack thereof) with which pollen grains are delivered at the same height at which they were picked up, whether each pollen grain is delivered to the destination flower's stigma, and whether each delivered pollen grain results in the fertilization of an ovule. Finally, stochasticity in the generation of offspring (due to mutation and to sexual reproduction) affects the phenotypes of offspring relative to their parents.

Observables
Several metrics were observed for each population in each generation of the realizations, including (1) the mean ecological trait value, (2) the mean values of the reproductive-organposition traits, (3) the mean magnitude of herkogamy (mean absolute difference in height between anther and stigma), (4) the mean female reproductive fitness (proportion of ovules fertilized), (5) the mean male reproductive fitness (number of pollen grains that fertilized an ovule, normalized by dividing by the number of ovules per plant), (6) the proportion of pollen taken from anthers, (7) the extent of style clogging due to self pollen, illegitimate pollen (from the same morph, and thus blocked by the dimorphic incompatibility mechanism), and legitimate pollen (from the opposite, compatible morph), and (8) the magnitude of reproductive isolation at fertilization, calculated as the number of ovules fertilized by resident pollen divided by the total number of ovules fertilized, which combines the effects of geographic isolation with the effects of mechanical isolation due to sexual selection against non-local pollen. The mean value of the ecological trait in each patch, z 1 and z 2 , at the end of each realization is particularly important because it allows us to evaluate the extent of local adaptation (or the lack thereof) exhibited by the plants in each patch as a consequence of the evolutionary dynamics experienced in each realization. All individuals were also observable graphically during model runs, including depictions of which flowers they were fertilized by and which flowers they fertilized, in order to allow both testing and exploration of the model [20].

Parameters
Parameter values governing the initial morphological distribution of the population (x i , y i , σ xi 2 , σ yi 2 ) were chosen to generate a unimodal distribution of reproductive-organ positions normally distributed around the center of the corolla tube, with a variance similar to that observed for natural populations [7]. (Initializing the model with no initial variance, σ xi 2 = σ yi 2 = 0, appears to make no difference, however, as the same equilibrium variance is rapidly attained in any case.) Parameter z i was chosen such that the initial population was equally maladapted ecologically (that is, with respect to the optimum) to both patches.
Parameter values governing the mutational variance (μ, α) were chosen to reproduce the same empirically observed variance in reproductive-organ position in subsequent generations. As is typical of individual-based models, a realistically low mutation rate would have resulted in the total loss of genetic variation at all loci, suggesting that processes other than mutation also act to maintain genetic variation in natural populations, but those processes are not well understood [21] and modeling them is in any case beyond the scope of the present research. The high mutation rate here, then, is equivalent to the fixed genetic variance assumed in many analytical models, and was not intended to be a realistic estimate of the mutation rate per se.
Parameter values governing the characteristics of the plants and patches (K, n o , n c , n p , u p , u s , m, v, σ p , σ s , l g , l σ ) were chosen with reference to personal observation of the "typical" values for heterostylous species in the genus Primula, but no attempt was made to precisely measure their values, or to make the model refer to any particular species. The model is not particularly sensitive to their values, although large changes to them can result in dynamics such as extreme pollen limitation that can have large effects.
The parameter values determining the difference in ecology between patches (θ 1 , θ 2 ) were fixed, representing a standardized ecological difference of 1.0 between patches. The difference between θ 1 and θ 2 may be regarded as a scaling factor defining the meaning of both the ecological trait values, z, and the strength of natural selection, ω. Variation of ω thus explores the full dimensionality of the parameter space here (given the initial condition of equal mean maladaptation to both patches).
Five parameters affecting the dynamics of natural selection due to ecology (ω) and sexual selection due to pollination (σ j , c, π 1 , π 2 ) were varied. The values used for these parameters are given in Table 1. Besides the values of σ j listed, "control" realizations of the model were conducted that simulated pollen transfer with no precision in height whatsoever; see Pollination phase for more details.
The parameters π 1 and π 2 represent "pollinator functions" that give the probability that a pollen grain will stick to the pollinator at a given height h in the interval [0,1]. Four pollinator functions were used in the realizations presented, defined as evidence for fine-scale differential stickiness is provided by Washitani et al. [22], in a detailed study of queen bumblebees pollinating Primula sieboldii. They document more than a 50-fold difference in the number of pollen grains stuck to different proboscis regions, with strong spatial segregation of pin and thrum pollen; however, the contribution of differential stickiness remains unclear since pollen was presented only at naturally occurring anther heights. Other studies have generally focused on overall pollinator effectiveness, not fine-scale differential stickiness; nevertheless, these studies do indicate substantial variation in stickiness among body parts and some degree of precision in pollen transfer [23][24][25][26][27][28][29][30]. Studies have also documented evidence for different levels of stickiness at high versus low levels of the corolla tube [29,31], thus modulating pollen transfer differently between high and low reproductive organs, resulting in outcomes such as the loss of style polymorphism [32] and the evolution of dioecy [33], but the generality of these findings is unclear. In short, although evidence for differential stickiness on pollinators exists, very little is known about the details in particular systems, or about how this translates into a quantitative precision of pollen transfer. Since our results demonstrate that these details are important to evolutionary outcomes, further empirical work on pollen transfer dynamics is needed.

Initialization
Each patch was initially seeded with K adult individuals. Each individual's x value was drawn from a normal distribution with mean x i and variance σ xi 2 , and the same was done for y

Germination phase
At the beginning of the germination phase, some of the ovules of the flowers in each patch have been fertilized by pollen in the pollination phase of the preceding year (see Pollination phase). In the germination phase, those fertilized ovules become seedlings, and those seedlings compete for the opportunity to survive to adulthood.
As implemented, a fertilized ovule is essentially a contract between two parent plants to Since one parent has S = 0 and one has S = 1, as enforced by the self-incompatibility mechanism of heterostyly mentioned previously, the value of the S locus of the offspring is randomly chosen as 0 or 1 with equal probabilities (modeling the empirical fact that a cross of Ss × ss produces, on average, 50% Ss and 50% ss).
All of the seedlings generated are added to the patch of their maternal parents. Competition then occurs, through soft selection based upon fitness derived from the ecological trait z, for the opportunity to mature to adulthood. If the total number of adult plants plus seedlings in a patch is less than or equal to the carrying capacity K, all seedlings survive to adulthood; otherwise, exactly enough seedlings will mature to fill the patch to carrying capacity. The probability of survival of each seedling i in patch j is proportional to its fitness W i , as defined by a Gaussian where z i is the ecological trait value of seedling i, θ j is the optimum for patch j, and ω is the strength of natural selection. Seedlings that do not survive to adulthood die, and are removed from the patch.

Mortality phase
At the beginning of the mortality phase, each patch contains only adult plants. During this phase, each plant dies with probability m. Plants that die are removed from their patch.

Pollination phase
At the beginning of the pollination phase, the reproductive state of the survivors of the mortality phase is reset (see Initialization). Pollination events are then conducted consecutively until the end of the pollination phase; the pollination season length v dictates the total number of pollination events conducted. Each pollination event consists of a set of steps executed sequentially: (1) The patch of the donor flower, 0 or 1, is chosen with equal probabilities.
(2) A determination is made as to whether this pollination event is a "crossover", in which a pollinator visits one patch and then the other, using the crossover probability c. If it is not a crossover, the patch of the recipient flower is the same as the patch of the donor flower; if it is a crossover, the patch of the recipient flower is the other patch.
(Note that pollinators that cross over do not then remain in their non-native patch; see (11) For each pollen grain delivered, the probability that it is received by the recipient flower's stigma is determined by a Gaussian function where h pollen is the height at which the pollen grain is delivered, h stigma is the height of the stigma of the recipient flower, and σ p is a scaling parameter determining how rapidly the probability of receipt of a pollen grain falls off with increasing difference between delivery height and stigma height. (16) At the end of each pollination event, a pollinator that has crossed over is assumed to return to its native patch (or to die, which amounts to the same thing since the supply of pollinators in the model is unlimited). This assumption is necessary, since otherwise the pollinators would soon equilibrate at equal frequency in both patches, undermining the very idea of a different pollinator native to each patch.
The preceding description of a pollination event uses Gaussian functions for the probability that a pollen grain at height h pollen will be received by a stigma at h stigma , during both selfpollination (step 4) and pollen delivery (step 11). However, since gravity would tend to pull pollen grains downward, the probability of pollen receipt might more closely resemble a lognormal function with its maximum where h pollen = h stigma . A pollen grain that arrives at or below a threshold l g below the stigma has a probability of zero of being received because pollen doesn't fall upward (the lognormal function is undefined here, but is taken to be zero). Above the maximum value of the function at h pollen = h stigma , the probability of receipt falls off asymptotically toward zero at a rate defined by the scale parameter of the lognormal, l o , reflecting the possibility that pollen can fall downward from any height above the stigma and have some nonzero probability of being received by it.
A version of the model incorporating this lognormal-based pollen delivery was constructed.
Specifically, in this version of the model the formula in step 4 is replaced by P jostle = u s L(l g + h pollen − h stigma , l 2 σ + ln l g ( ) , l σ ) L l g , l 2 σ + ln l g The quotient form of the formulas for P jostle and P receipt serves to normalize the height of the functions to 1 when h pollen = h stigma . The value l 2 σ + ln l g ( ) for the μ parameter of the lognormals serves to locate the peak of the functions at x = l g .
Results from the lognormal version of the model were generally qualitatively similar to results from the Gaussian version of the model. In some ways the lognormal model did appear to better match empirical data; in particular, thrums received more self pollen than pins in this version of the model, increasing their male function and decreasing their female function relative to pins, a phenomenon which has often been observed in nature [37][38][39]. These dynamics did not substantially affect the results of the model presented here, however, and so all results presented are taken from the Gaussian version of the model, since it is conceptually simpler. In general, the dynamics of pollination are remarkably complex [e.g., 40,41,42], and attempting to introduce all of this complexity into a model would be premature; we have strived for a balance that includes only that complexity necessary to pursue the questions at hand.