1. Introduction

Understanding plant sexual reproduction requires quantitative measurements of mating and investigation of the ecological context in which it occurs (Lloyd, 1980; Eckert et al., 2010; Barrett and Harder, 2017). Animal pollinators frequently move among nearby neighbours and pollen dispersal patterns can be highly leptokurtic (Levin and Kerster, 1974; Harder and Barrett, 1996; Ghazoul, 2005). As a result, plant mating is often promiscuous and nonrandom. Both extrinsic (ecological) and intrinsic (plant traits) factors interact to influence nonrandom mating patterns. Where nonrandom mating occurs, sexual partners may be phenotypically more (assortative mating) or less similar (disassortative mating) to one another than would occur by chance (Barrett and Harder, 2017). Both forms of nonrandom mating have important ecological and evolutionary consequences. For example, disassortative mating and negative frequency-dependent selection, which can be driven by discrete floral traits, are major selective forces that maintain several sexual polymorphisms in flowering plants e.g., heterostyly (Lloyd and Webb, 1992; Barrett and Hodgins, 2006). In contrast, positive assortative mating can lead to reduced genetic diversity and increased genetic structure within populations, which can have diverse effects on a population, including limiting the adaptive potential (Turelli and Moyle, 2007; Merrill et al., 2019) or increasing reproductive isolation between phenotypic groups, promoting local adaptation but also possibly leading to genetic divergence and reduced gene flow (Minnaar et al., 2019b).

Previous studies have examined several traits that drive assortative mating in plant populations, including flowering time (Weis et al., 2005; Devaux and Lande, 2008; Ison and Weis, 2017; Rubin et al., 2019; Godineau et al., 2022), flower colour, and polyploidy (Kulkarni, 1999; Jones and Reithel, 2001; Kennedy et al., 2006). Assortative mating may also be caused by the density and spatial distribution of plants in mating neighborhoods, as these factors influence pollinator preferences. For example, mating patterns among floral-tube morphs of Gladiolus longicollis have been shown to be density-dependent, resulting in assortative mating at low density but random mating at high density (Rymer et al., 2010). However, the extent to which floral traits with quantitative variation influence patterns of nonrandom mating has not been investigated in detail.

Nectar spurs are tubular outgrowths from petals or sepals containing nectar at their base (Hodges, 1997; Fernández-Mazuecos et al., 2018). Access to nectar at the base of nectar spurs requires depth-probe feeding by long-tongued pollinators, thus restricting the number of pollinators that can obtain floral rewards, a trait shared with long floral tubes in many species (Nilsson, 1988). These specialized floral structures are often adapted to specific pollinators, with certain plant species relying on long-tongued bees or butterflies to access nectar, thus creating a specialized pollination system (Fenster et al., 2004). Both structures play an important role in determining the sites of pollen deposition on pollinators or pollen export (Campbell et al., 1991; Lloyd and Webb, 1992; Anderson et al., 2010). Nectar spurs and floral tubes generally exhibit quantitative variation; however, cases are known in which floral tubes display bimodal variation within populations (Anderson et al., 2016), and this has been linked to pollen transfer between floral tubes of similar length and potentially assortative mating (Minnaar et al., 2019b). Although nectar spur length is usually distributed unimodally within natural populations, the extent to which quantitative variation in nectar spur length affects mating patterns remains unclear.

Aquilegia (Ranunculaceae) is a genus of approximately 70 animal-pollinated perennial species that usually exhibit prominent nectar spur-lengths, ranging from a few centimetres to 15 cm. This variation has been shown to be heritable (Hodges et al., 2002; Edwards et al., 2021) and associated with pollinator-driven speciation, with the evolution of elongated nectar spurs driven by pollinator shifts from bee to bird to hawkmoth pollination (Whittall and Hodges, 2007); thus, depending on local ecological conditions nectar spur length should be responsive to selection imposed by pollinators.

Genetic markers have provided a wealth of information on the ecological and life-history correlates of plant mating patterns (Barrett and Eckert, 1990; Goodwillie et al., 2005; Munoz et al., 2016; Moeller et al., 2017). Earlier studies in Aquilegia species reported mixed mating, with a wide range of population-level estimates of maternal outcrossing (Herlihy and Eckert, 2004; Brunet and Sweet, 2006; Yang and Hodges, 2010). However, these studies have involved a single axis of mating-system variation involving the mean frequency of cross- and self-fertilization of maternal parents within populations – the classic outcrossing-selfing paradigm (Barrett, 2003). A more comprehensive approach to the measurement of plant mating includes a second axis of variation – mate diversity – encompassing the number, frequency and phenotypic traits of mates that plants couple with as maternal and paternal parents. Measuring both axes of variation enables the examination of ‘mating portfolios’ (Garcia-Gonzalez et al., 2015; Barrett and Harder, 2017), which include all offspring to which plants contribute genetically as maternal or paternal parents. Despite the probable ecological and evolutionary significance of mating portfolios, few studies (e.g., Christopher et al., 2019; Yuan et al., 2023) have attempted to quantify female and/or male mate numbers in plant populations.

Most Aquilegia species exhibit the bee-pollination syndrome, characterized by blue-purple flowers, petals with long blades for landing and relatively short curved nectar spurs (Edwards et al., 2021). A typical example of this syndrome is Aquilegia rockii Munz (1946), which is endemic to the mountains of southwestern China, prevalently in Yunnan province. Our preliminary observations confirmed that this species is primarily pollinated by bumble bees (Fig. 1) and displays a unimodal distribution of spur-length with populations (Fig. S1). Here we use A. rockii to investigate mating patterns in a natural population using open-pollinated seed families from mapped plants and microsatellite markers. A goal of the present study was to extend these findings to include mating at the individual plant level through both female and male functions, and determine whether variation in nectar spur length might influence mating patterns. To our knowledge, no study has investigated whether nectar spur length is associated with mating patterns. Our investigation addressed the following specific questions: (1) How much variation in maternal and paternal outcrossing and mate diversity occurs among plants, and to what extent are mating parameters intercorrelated? (2) Is there evidence that nectar spur length and floral display size influence mating system components, including patterns of non-random mating and mating distance between sexual partners? (3) If quantitative floral traits are associated with nonrandom mating pattern, what is the underlying mechanism? Our analysis of mating portfolios revealed evidence for assortative mating for spur length, and we consider the mechanisms and potential evolutionary consequences of this novel finding.

2. Materials and methods 2.1. Study system and site

Aquilegia rockii flowers between June to August with most plants possessing one shoot with a cymose inflorescence of 1–15 pendulous hermaphrodite flowers with five purple sepals and petals, each with a conspicuous nectar spur (Fig. 1). Spurs generally vary in length from 17 to 27 mm. Flowers are protandrous with a male phase of approximately 12 h and floral longevity ranges from 4 to 7 days. Plants are self-compatible and produce seeds by autonomous self-pollination, ensuring reproductive success when pollinator activity is limited or outcrossing opportunities are scarce.

We investigated a population of Aquilegia rockii located about 0.5 km northwest of Shangri-La Alpine Botanical Garden, NW Yunnan, China (27.89°N, 99.66°E, 3535 m a.s.l.). Plants in the population typically display 2–3 flowers per day and a total of 3–9 flowers during the season. The main pollinators of A. rockii are Bombus sp. bumblebees (Figs. 1b and S1) and Bombus minshanicola Bischoff (Bischoff 1936; Bertsch and Schweer, 2011; Fig. 1c). Both bees behave similarly when probing spurs for nectar, and during this process they touch the sexual organs of flowers.

2.2. Measurement of floral traits and sampling for estimates of mating patterns

We examined two floral traits in this study: spur length, a floral design trait, and the total flowers per plant, an agent of daily display size during the flowering season. In late June 2018, we measured spur length and total flower number per plant for all plants (n = 214) in the population (Fig. S2). We mapped 67 flowering plants in a 50 m × 40 m quadrat within the population to investigate the effects of floral traits on mating patterns. For each plant, we counted all flowers, including withered flowers, open flowers, and flower buds. We then selected 1–3 fresh flowers, depending on the number of open flowers and measured the length of the spurs, from the distal end of the spur to the petal attachment point, using digital calipers (L.S. Starrett Company, 799A-6/150, precision 0.01 mm).

Thirty-one (46.27%) individuals from the quadrat were selected as evenly as possible, with at least 2 m distance between each other (Fig. S3). For each selected individual, we sampled three mature fruits considered as a family, on each plant during early August. Ten seeds per fruit were sown in separate Petri dishes and placed in a growth incubator at 20℃ in the laboratory. After three weeks, we obtained, on average, approximately 90% seed germination per fruit. We then randomly selected nine seedlings (three seedlings per fruit) and collected fresh leaves from each seedling. This sample size was selected because the effective mate number in previous studies generally ranged from two to four (reviewed in Mitchell et al., 2005). In total, we sampled 275 offspring for further parentage analysis, and an additional four offspring from three maternal plants were excluded due to allele errors. Additionally, we collected leaf tissue from the 67 mapped plants within the quadrat and 22 plants located within 10 m surrounding the quadrat to assess their potential as paternal parents. Sampled leaf tissue from each individual was dried in silica gel for genotyping.

2.3. DNA extraction, microsatellite genotyping, and paternity assignment

We extracted total genomic DNA from adult and offspring leaf tissue using the Plant Genomic DNA Kit (TIANGEN, Beijing, China). We used 12 pairs of microsatellite polymerase chain reaction (PCR) primers (Table S1) to score the genotypes of the DNA samples. We conducted PCR amplifications on a thermocycler (PerkinElmer, Foster City, CA, USA) using a 25 μL reaction volume containing 20 ng of genomic DNA template. We then checked the PCR products on 1% agarose gels stained with ethidium bromide. We separated all PCR products and visualized them using an ABI 3730 XL automated sequencer (Applied Biosystems, Foster City, CA, USA). We determined allele sizes using GENEMAPPER 4.0 (Applied Biosystems) before rechecking manually to reduce scoring errors.

To estimate the probability that each seedling from the 31 focal maternal parents resulted from selfing or outcrossing, we used COLONY 2.0.7 (Jones and Wang, 2010; Wang et al., 2012), a program that jointly infers parentage and sib-ships by maximizing the likelihood of the inferred pedigree. The main parameters for each run were as follows: monoecious and diploid species, both sexes polygamous, with inbreeding allowed. The full-likelihood method was employed with medium likelihood precision, using the default of sib-ship prior. Allele frequencies were not updated due to the small family size, and each run was repeated ten times with a default random number of seeds. We permitted an error rate of 0.01 for allelic dropouts and false alleles at each locus.

2.4. Mating portfolios and pollen dispersal

COLONY successfully assigned the most likely paternal parent to 200 (72.7%) of the seeds, with another 75 (27.3%) seeds sired from plants located outside both the quadrat and the adjoining 10 m area. We next used COLONY to analyze components of the mating portfolios of various subsets of the 275 genotyped seeds. The analysis of female outcrossing rate considered all genotyped seeds for which the probability of selfing was > 0.95 (i.e., selfed). Therefore, a total of 118 seeds were self-fertilized. The analysis of the number of outcrossed male mates for maternal plants involved the remaining 157 outcrossed seeds with a selfing probability < 0.05. Potential sires of these seeds included identified genotyped paternal parents (82 seeds) and immigrant paternal parents located outside the quadrat (75 seeds). Specifically, when the probability that the seeds belonged to a full-sib family (same parents) was > 0.95, these outcrossed offspring within a known maternal parent were considered likely to be the same sires. When the probability of seeds belonging to a full-sib family was < 0.95, it indicated lower certainty that these seeds originated from the same pair of parents. In these instances, outcrossed seeds sired by immigrant pollen were assigned to a hypothetical ‘sole dummy’ sire. This method facilitated the estimation of paternal genetic contributions when the exact sire was unidentified. We then assessed the corresponding male parameters, including male outcrossing rate and the number of outcross female mates for paternal plants, using 200 seeds with identified paternal parents in COLONY.

We next quantified the similarity of outcross mate composition between the paternal and maternal functions of each sampled plant using the Jaccard index following Christopher et al. (2019). The index ranges from 0 to 1, with larger values indicating higher similarity between the two sexual functions of the same plant. We then visualized the resulting mating networks using the bipartiteD3 package in R (Terry, 2018).

To determine whether pollen dispersal distances based on parentage analysis were influenced by floral traits and the spatial relations of flowering individuals, we estimated mating distances from the spatial location of maternal and paternal parents for each of the 194 mating events. This included 97% of 200 seeds with identified maternal and paternal patents, as well as six seeds in which the paternal parents were outside the study quadrat.

2.5. Pollinator observations

To investigate whether pollinator foraging behavior is influenced by spur length and daily floral display size, and how this might affect mating outcomes, we established three 3 m × 3 m plots outside the 50 m × 40 m quadrat, with each separated by at least 50 m, which might reduce repeated visits by the same individual and minimize confounding effects on our measurements of foraging behavior. We conducted pollinator observations on sunny days from June 13 to 19, 2019. In each plot, we measured the average spur length (three flowers per individual) of each plant and recorded the number of open flowers on the day, which is referred to as the daily floral display size. Pollinator observations were conducted by two observers from 8:00 to 18:00 on three consecutive days in each plot, with each plot observed for a total of at least 30 h (10 h per day). When a bumblebee arrived in a plot, we recorded the species, the sequence of marked plants visited by each bee, and the number of flowers visited on each plant.

2.6. Spur-length variation and pollination success

We predicted that if spur-length variation influenced mating patterns, we might expect plants to produce different amounts of outcrossed seeds. Conversely, we predicted that variation in spur length would be unlikely to influence the amount of autonomous seed set produced by individuals. To investigate these hypotheses, we performed a manipulative pollination experiment in mid-June 2019 on our study population. We selected 60 individuals with similar daily display sizes and measured the spur length of all open flowers on each individual on five consecutive days. We then evenly and randomly assigned individuals to three treatments with 20 plants per treatment: (1) open (natural) pollination, in which we marked at least three flowers on each individual and left them intact and exposed to pollinators; (2) emasculation, in which at least three flowers on each individual were emasculated before anther dehiscence to prevent self-pollination and then exposed to pollinators; (3) supplemental hand pollination, in which at least three flowers on each individual were marked and exposed to pollinators but also received hand cross-pollination. After five weeks, we collected mature fruits derived from the three manipulated flowers per plant before follicles started to split and counted all fully developed seeds, aborted seeds, and possible unfertilized ovules in each fruit. We calculated seed set per flower as the number of fully developed seeds divided by the number of total ovules (including fully developed seeds, aborted seeds, and unfertilized ovules).

2.7. Statistical analyses

We conducted all statistical analyses in R v.4.2.2 (R Core Team, 2022). We assessed the relationships between floral traits (spur length and total flower number per plant) and mating parameters (female or male outcrossing rate, female or male outcross mate number, and mating distance) using generalized linear mixed models (GLMMs) implemented with the glmmTMB function from the glmmTMB package. The models were constructed as follows: Model-1, in which spur lengthˆ2 was included as a fixed effect and total flower number per plant as a random intercept; Model-2, in which only spur lengthˆ2 was used as a fixed effect; and Model-3, which included spur lengthˆ2, total flower number, and their interaction as fixed effects. Additionally, to account for the potential influence of total flower number per plant on mating parameters, we constructed two further models: Model-4, where total flower number per plant was included as a fixed effect and spur length as a random intercept, and Model-5, which included only total flower number per plant as a fixed effect. In all models, mating parameters served as the dependent variable. To test whether the mating parameters were affected by the frequency of a given spur length in the population, we first described the distribution of quantitative traits and calculated frequency distributions of spur length. We divided the spur length data into 40 equal-interval classes and visualized the distribution using histograms generated with the stat.test function. Next, we performed generalized linear models (GLMs) to evaluate whether the frequency of a given spur length affected mating parameters. We used the glmmTMB function from the glmmTMB package to fit the models, with the spur-length frequency as an independent variable and mating parameters as the dependent variable. The specific sampling distributions and link functions used for each model are detailed in the supplementary materials (Tables S2 and S3).

To evaluate whether floral traits affected pollen dispersal, we examined absolute differences in standardized traits between paternal and maternal parents (simulated a pollinator moving from paternal pollen donor to maternal pollen recipient) with a paired Wilcoxon test using the wilcox.test function from the dplyr package. To investigate the potential for assortative mating between mate pairs on floral trait classes, we performed Spearman rank correlation analysis of floral traits for each mate pair for the 82 outcrossed mating events and for all 200 mating events (including selfing and outcrossing) with the cor.test function in R.

To investigate the relationship between floral traits and variation in bumblebee visitation at the plant and flower levels, we calculated visitation rate (the number of plants or flowers visited in an hour). As a preliminary screening of potential associations, we first performed Spearman correlation analyses to examine the relationships between spur length and visitation rate, and between daily floral display size and visitation rate by using the cor.test function. This non-parametric approach allowed us to quickly assess whether these floral traits showed significant monotonic correlations with visitation rate. We then employed GLMMs to more rigorously evaluate how floral traits influence visitation rate while controlling for other traits and pollinator species as random effects. Specifically, we constructed three separate models in glmmTMB function with Tweedie distribution (link = “log”) from the glmmTMB package: 1) Model-1 treated spur length as the fixed effect, with daily floral display size and bumblebee species as random intercepts; 2) Model-2 treated daily floral display size as the fixed effect, with spur length and bumblebee species as random intercepts; 3) Model-3 treated spur length and daily floral display size as the fixed effect, with bumblebee species as random intercepts. Details of these models are provided in the supplementary materials (Table S4). Additionally, after confirming normality, we used the lm function to analyze the number of flowers visited per plant per foraging bout (across three plots and for each bumblebee species) and the total flower number per plant. These analyses were conducted to explore how daily floral display size relates to short-term pollinator foraging behavior and to assess whether total flower number over the season serves as a suitable indicator for evaluating longer-term reproductive outcomes associated with the relationship between daily floral display size and short-term pollinator foraging behavior.

To evaluate the relationship between spur length and reproductive success, we performed generalized linear models (GLMs) using the glmmTMB function with Beta_family distribution (link = “log”) from the glmmTMB package. First, we assessed the effects of different pollination treatments on seed set by constructing models in which spur length and/or spur lengthˆ2 were included as independent variables, with seed set under different pollination treatments as the dependent variable. Second, we examined the relationship between spur length frequency (independent variable) and seed set across different pollination treatments (dependent variable). The detailed models are provided in the supplementary materials (Table S5).

3. Results 3.1. Variation in floral traits in the natural population and experimental plants

Variations in spur length and total flowers per plant for the natural population, the 50 m × 40 m plot, and the 31 individuals used as parents exhibited similar mean values and ranges (Figs. S1 and S2; Table S6). Here we provide values only for the natural population. The mean (±s.e., hereafter) spur length was 23.66 ± 0.14 mm (n = 214, range 17.58–27.83 mm), and variation in spur length within plants was lower than among the 19 observed plants (Table S7). Spur length exhibited a leptokurtic distribution and was skewed towards longer spurs (Fig. S1a). The mean total flower number per individual was 6.1 ± 0.4 (n = 214, range 1–45) and exhibited a possible log-normal or gamma distribution skewed towards a higher number of flowers per plant (Fig. S1d). Correlation analysis revealed no significant relationship between spur length and total flower number at the individual level (r = −0.041, P = 0.823).

3.2. Female and male outcrossing rate and floral traits

Among parental plants, the female and male outcrossing rates averaged 0.573 ± 0.049 (range 0.111–1) and 0.477 ± 0.066 (range 0–1), respectively. A significant negative quadratic relationship was observed between spur length and female outcrossing rate (P = 0.002; Fig. 2a and Table 1), but was not observed between spur length and male outcrossing rate (P = 0.919; Fig. 2b and Table 1). In addition, a significant negative relationship was observed between total flower number and female outcrossing rate (P = 0.006; Fig. S4a and Table S2), whereas no significant relationship was observed between total flower number and male outcrossing rate (P = 0.224; Fig. S4b and Table S2). Furthermore, the GLM results revealed that spur length was significantly associated with female outcrossing rate but not male outcrossing rate when total flower number was also assigned as the independent variable (Table S2).

3.3. Female and male mate number and floral traits

The mate composition between male and female functions varied significantly among the majority of the plants (Fig. 3). The mean Jaccard similarity between the outcross male and female function was 0.043 (range = 0–1, n = 31), indicating a very low level of mate overlap between alternate sexual functions of the same individual in the 31 sampled plants. Mean female and male outcross mate numbers were 3.5 ± 0.3 (range 1–8) and 1.4 ± 0.2 (range 0–4), respectively. The relationships between floral traits and mate number paralleled those for outcrossing rate. Indeed, there were significant positive linear relationships between female outcrossing rate and female outcross mate number (P < 0.001; Fig. S5a) and also between male outcrossing rate and male outcross mate number (P < 0.001; Fig. S5b).

A significant negative quadratic relationship was observed between spur length and female outcross mate number (P = 0.009; Fig. 2c and Table 1), but the quadratic relationship between male outcross mate number and spur length was not significant (P = 0.928; Fig. 2d and Table 1). Similarly, significant relationships were not evident between total flower number and female (P = 0.083; Fig. S4c and Table S2) or male outcross mate number (P = 0.268; Fig. S4d and Table S2), respectively. The GLM results indicated that the spur length was significantly associated with female outcross mate number but not male outcross mate number when the total flower number was also assigned as the independent variable (Table S2).

3.4. Analysis of distance between mating pairs

We detected a total of 200 mating events by paternity assignment in the 67 individuals mapped in the experimental quadrat. Most of these events (194, 97%) occurred within the quadrat. The frequency distribution of mating distance based on 194 mating events indicated that most (164, 82%) involved parents that were 10 m or less apart (Figs. 2e and S4e). After excluding selfing, the mean mating distance was 10.778 m (±0.125 m; range: 0.462–38.773 m). Regression analysis indicated a significant negative quadratic relationship between spur length and mating distance (P < 0.001; Fig. 2e and Table 1) and a significant negative linear effect of total flower number on mating distance (P = 0.012; Fig. S4e and Table S2). When total flower number was used as a random intercept, the GLMM results revealed that spur length significantly affected mating distance (Table S2).

3.5. Mating patterns and variation in spur-length frequency

In our plots, there was a significantly positive relationship between spur-length frequency and female outcrossing rate (P < 0.001; Fig. 4a and Table S3), but not with male outcrossing rate (P = 0.367; Fig. 4b and Table S3). Moreover, there was a significant positive linear relationship between spur-length frequency and female outcross mate number (P = 0.008; Fig. 4c and Table S3) and with distance between mating pairs (P = 0.001; Fig. 4e and Table S3), but not with male outcross mate number (P = 0.520; Fig. 4d and Table S3).

Paternity analysis of 200 mating events revealed a significant difference in spur length (V = 12,958, P < 0.001; Fig. 5a) and total flower number (V = 13,420, P < 0.001; Fig. 5b) between mating pairs. Spearman rank correlation analysis indicated a significant positive correlation in spur length (r = 0.759, P < 0.001; Fig. 5c) and total flower number (r = 0.552, P < 0.001; Fig. 5d) between mating pairs regardless of whether the mating was selfing or outcrossing. When selfing events were excluded, a significant positive correlation in spur length between mating pairs still persisted (r = 0.270, P = 0.014; Fig. 5c), but no significant correlation was found for total flower number (r = 0.013, P = 0.907; Fig. 5d).

3.6. Pollinator observation

Two species of bumblebees (Fig. 1b and c) made 102 visits to Aquilegia rockii in our 30–40 h of observations. Bombus sp. and B. minshanicola accounted for 78.43% and 21.57% of total visits, respectively. The average spur length of plants had a significant positive correlation with visitation frequency when data from both bee species were combined at both flower and plant levels (Fig. S6a, c). However, the GLMM results indicated that visitation rates were not influenced by spur length when daily floral display size and pollinator species were treated as covariates (Tables 2 and S4). Instead, visitation rates were primarily influenced by daily floral display size (Table S4), suggesting that spur length is not the main factor affecting pollinator choice. Daily floral display size per plant was also positively correlated with visitation rate at both flower and plant levels (Fig. S6b, d). We found a positive relationship between the number of flowers visited per plant by bees and daily floral display size (P < 0.001; Fig. S7).

3.7. Pollination success and spur length

Spur length did not significantly affect the seed set of naturally-pollinated or hand-pollinated flowers (Fig. 6a and Table S5). However, a significant negative quadratic relationship was observed between the average spur length per plant and seed set when flowers were emasculated (P = 0.012; Fig. 6a and Table S5). We also found that spur-length frequency per plant was positively related to seed set in the emasculation treatment (P < 0.001; Fig. 6b and Table S5) but not in open-pollinated flowers (P = 0.575; Fig. 6c and Table S5).

4. Discussion

Consistent with earlier population-level studies of Aquilegia species (Herlihy and Eckert, 2005; Brunet and Sweet, 2006; Yang and Hodges, 2010), our analysis of maternal outcrossing in A. rockii revealed mixed mating. This description, however, masks considerable mating complexity at the plant level resulting from differences between sexual functions in outcrossing rate and mate number, some of which we show are associated with nectar spur length. Moreover, our parentage analysis provides evidence of assortative mating for spur length. Below we consider the ecological and evolutionary consequences of the mating complexity. Then, we report and evaluate several hypotheses concerning the potential mechanisms influencing assortative mating.

4.1. Associations between floral traits and outcrossing rates

In self-compatible angiosperms, maternal outcrossing rates exhibit nearly continuous variation, with several floral traits playing crucial roles in shaping this variation (Whitehead et al., 2018). For instance, the spatial separation between stigmas and anthers, a key trait known as herkogamy, has long been recognized for its influence on outcrossing rates (Barrett and Shore, 1987; Brunet and Eckert, 1998; Takebayashi et al., 2006; Opedal, 2018). A comparative analysis across species has shown a positive correlation between floral attraction, measured as the product of flower size and flower number, and outcrossing rates (Goodwillie et al., 2010). However, in Aquilegia species, no clear association has been observed between flower size and outcrossing, either within (Brunet and Eckert, 1998) or among populations (Herlihy and Eckert, 2005). Additionally, while large daily floral displays are often thought to promote selfing due to geitonogamy (Harder and Barrett, 1995; Karron and Mitchell, 2012), empirical evidence on this relationship remains mixed. For example, in Aquilegia coerulea, Brunet and Sweet (2006) reported a positive association between outcrossing rate and floral display size, probably as a result of increased geitonogamy in large displays. In contrast, Herlihy and Eckert (2004) found no such effect in Aquilegia canadensis, where geitonogamy was minimal. Our study similarly uncovered a significant negative correlation between total flower number per plant and the outcrossing rate of individuals (Fig. S4a). While we did not directly measure daily display size for the plants from which these estimates were obtained, flower number per plant was strongly correlated with daily floral display (P < 0.001; Fig. S8), a relationship also observed in other Aquilegia species (Brunet and Eckert, 1998). It suggests that flower numbers per plant can act as an agent of the floral display size for a particular individual, as our interest was not to estimate the effect of the floral display size on mating system parameters over a few days.

Despite the relatively limited variation among plants in the spur length of Aquilegia rockii (Fig. S2), especially when compared to interspecific patterns, our study highlights that spur length is an important trait that influences mating patterns. Correlation analysis revealed significant associations between spur length and outcrossing rate at the plant level, with individuals of intermediate spur length exhibiting higher outcrossing rates than those with either shorter or longer spurs (Fig. 2a). Furthermore, we found that variation of spur length influences the distance among mating pairs (Fig. 2e). When we included total flower number as a random intercept in the model, although the generalized linear mixed effects model failed to converge, the random effect variance for total flower number was negligible. This suggests that the observed effects were largely driven by spur length rather than flower number (Table S2). Significantly, the correlation between spur length in mating pairs persisted even after excluding selfing from the analysis (Fig. 5c), further supporting the importance of spur length in mating dynamics. While the pronounced variation in spur length among Aquilegia species and its functional role in promoting macroevolutionary diversification has been well documented (Whittall and Hodges, 2007), our study is the first to demonstrate that even restricted population-level variation in spur length can potentially influence mating patterns.

Our estimates of female (t = 0.57) and male (t = 0.48) outcrossing rates suggest that nearly half of all seeds in the Aquilegia rockii population result from self-fertilization, regardless of whether ovule or pollen function is considered. This substantial selfing rate raises important questions about the underlying mechanisms that drive selfing and the survival of selfed offspring. Selfing in A. rockii is likely facilitated by multiple mechanisms, including intrafloral and geitonogamous pollination. To fully assess the relative contributions of these mechanisms to the overall selfing rate, manipulative experiments, such as emasculation, are necessary. These experiments, similar to those conducted in A. canadensis (Herlihy and Eckert, 2004), allow us to decompose the selfing rate into its component parts, distinguishing between intrafloral selfing and geitonogamy. Additionally, evaluating the potential role of early-acting inbreeding depression is crucial for understanding the long-term effects of selfing on the mating system and the survival of selfed progeny. In several Aquilegia species, early-acting inbreeding depression has been documented, with selfed progeny failing to reach flowering due to strong selective pressures against inbred individuals (Herlihy and Eckert, 2004; Karron and Mitchell, 2012). This phenomenon may also occur in A. rockii, which would significantly affect the role of selfed offspring in the mating system. Thus, assessing the fitness of selfed progeny and their potential participation in future mating events is essential for understanding the evolutionary consequences of selfing. If inbreeding depression is present, selfing could reduce genetic diversity and limit adaptive potential, particularly in small or fluctuating populations, as has been observed in other species (Cheptou and Donohue, 2010; Cutter, 2019).

4.2. Variation in female and male mate number

Relatively few studies of plant mating have measured the number of female and male mates for individual plants within natural populations, although these have been proposed to likely vary considerably (Pannell and Labouche, 2013). Some information on the average female mate number, m_e, is available from analyses that measure outcrossed correlated paternity (Ritland, 1989). This parameter involves the average number of mates per female, assuming all males contribute independently and equally to siring ovules of offspring. A survey of 66 angiosperm species revealed m_e increases with female outcrossing rate, as one might predict (Fig. 3a in Barrett and Harder, 2017). Although there is abundant evidence, beginning with Darwin's seminal work on inbreeding depression, on the fitness benefits of outcrossing as opposed to selfing (Charlesworth and Charlesworth, 1987), there is a paucity of data on the fitness benefits of mate diversity. Multiple mating should increase mean family quality and thus parental fitness (Nora et al., 2016) and may serve as a bet-hedging strategy in variable environments (Garcia-Gonzalez et al., 2015), but empirical work on the fitness consequences of family variation in mate diversity, and the extent to which it differs between the sex functions as we have shown (Fig. 3), has yet to be examined in detail.

We found that although mate number for plants was relatively low, it differed between the sexual functions (female = 3.5 ± 0.3, range 1–8; and male = 1.4 ± 0.2, range 0–4), with on average more than twice as many mates resulting from female than male outcrossing. The generally low mate number is not unexpected because nearly half of all ovules were self-fertilized (s = 0.43). Indeed, both female and male outcrossing rates were strongly correlated with female and male mate numbers (Fig. S5). Christopher et al. (2019) reported relatively similar numbers of mates through female (3.8) and male (4.4) functions in an experimental population of Mimulus ringens. In their study, all plants acquired fitness through both sex functions, whereas we found some individuals were unsuccessful at male outcross siring success. Consistent with their study, we found that there was minimal overlap between the number of outcross mates that were shared between the female and male functions of individuals (Fig. 3), demonstrating how hermaphroditic sex expression can serve to increase mate diversity. The mate number and the overlap between male and female function were probably underestimated in our study due to sampling a relatively small proportion of the total seeds produced per plant. Similar to Mitchell et al. (2005) and Christopher et al. (2019), it is likely that mates have been largely sampled in our study, as pollen carryover may be very limited in species pollinated by bumblebees due to behaviors such as grooming and trapline foraging, which result in greater loss of conspecific pollen (Holmquist et al., 2012; Gamba and Muchhala, 2023).

Few studies have investigated associations between floral traits and mate diversity. Our analyses of the relation between mate number and spur length in A. rockii revealed a complex pattern mirroring our findings for female and male outcrossing rate. There was a significant positive relationship between female mate number and spur length for plants with short to near medium length spurs but a negative relation for plants with very long to intermediate values (Fig. 2c). In contrast, there was no relationship between spur length and male mate number (Fig. 2d). A previous study of the potential influence of herkogamy on mate diversity in Narcissus longispathus reported an increasing number of mates with greater herkogamy despite similar outcrossing rates among the three herkogamy classes that were investigated (Medrano et al., 2012). Although it is intuitive that mate diversity should increase with outcrossing rate, further studies are required to determine the pollination mechanisms and floral traits that jointly influence these two key mating parameters.

4.3. Assortative mating – mechanism and consequences

Assortative mating based on spur length was evident in our study population of Aquilegia rockii (Fig. 5c). Several mechanisms can lead to assortment. For example, selfing is a form of assortative mating, and on average close to half of all seed in our study population were self-fertilized, although plants with mid-length spurs were more outcrossed than those with shorter or longer spurs (Fig. 2a); however, we still detected significant assortative mating when selfing events were excluded from our analysis (Fig. 5c), indicating that most outcrossed mating events occurred between individuals with similar spur lengths.

Assortative mating in Aquilegia rockii may occur if the types of pollinators differed among flowers with different spur lengths, or pollinators discriminated among spur-length phenotypes exhibiting pollinator preference, as occurs for other floral traits (e.g., Waser, 1986; Epperson and Clegg, 1987; Kennedy et al., 2006; Rymer et al., 2010). Our observations revealed only two effective pollinators of similar size, Bombus sp. and B. minshanicola (Fig. 1b and c) but, on average, they showed no apparent preference for particular spur-length phenotypes (Table 2 and Fig. S6a, c). However, bumble bee species show considerable body-size variation (Fitzgerald et al., 2022) and we cannot rule out the possibility that individuals of different sizes might display preferences for floral traits.

Another mechanism that can lead to assortative mating includes differences in pollen placement on the bodies and tongues of bumble bees that varies consistently with spur length. Indeed, studies have suggested that this mechanism cause assortative pollen receipt by phenotypes of contrasting floral-tube length (Minnaar et al., 2019b), as well as phenotypic divergence leading to speciation (Maad and Nilsson, 2004). However, observations from our unpublished study indicate that pollen is deposited on numerous sites on bumble bees visiting flowers of Aquilegia rockii, even in a population with bimodal spur length (our unpublished data). Therefore, it seems unlikely because of the diffuse pollen placement on bumble bees that this mechanism promotes assortative mating for spur length.

Two additional mechanisms that may promote assortative mating can be ruled out as well. First, we observed no significant spatial autocorrelation for spur length among mating pairs (Fig. S10), suggesting that nonrandom mating is unlikely to be caused by the clumped distributions of spur-length phenotypes. Second, although differences in flowering time have been commonly implicated in assortative mating (Weis et al., 2005), this mechanism is not relevant because all plants chosen for mating system analysis flower during the same period. Mating outcomes result from multiple factors acting sequentially during pollination and post-pollination processes (Anderson et al., 2016; Barrett and Harder, 2017; Minnaar et al., 2019a). Further studies on pollen dispersal and pollen-tube competition among spur-length phenotypes of A. rockii are necessary to determine how assortative mating is achieved.

Our study showed significant positive relationships at the plant level between the frequency of spur length and maternal outcrossing, siring success, mate distance (Fig. 4), and seed set of emasculated flowers (Fig. 6a and b). Thus, mating patterns in our study population were influenced not only by spur length itself, but also by the frequency of spur-length phenotypes. This suggests that differences in mate availability among individuals may play a role in mating. Our finding that opened flowers generally set less seed compared to supplemental pollinated flowers (Fig. 6a) suggests reproductive success is limited by mating opportunity, with autonomous self-pollination providing a certain degree of reproductive assurance. Thus, the mechanism underlying frequency-dependent relationships between spur length and mating traits may be mating availability, as spur length was not found to affect floral constancy or pollinator visitation (Tables 2 and S4). Although spur length is unlikely to influence the amount of autonomous seed set, plants with the most frequent mid-length spurs appear to benefit by receiving more outcrossed pollen than rarer spur-length phenotypes (Fig. 6), suggesting mid-length spurs have more outcross mating opportunities. This is because the largest number of individuals with mid-length spurs within populations results in relatively more mating opportunities when positive associative mating occurs.

Evidence of associations between the frequency of spur length and female outcrossing rate, mate number, mating distance and outcross opportunity are consistent with the hypothesis that mate availability may be a cause of positive frequency-dependent mating in Aquilegia rockii (Figs. 4cee and 6a, b). Furthermore, this process causes assortative mating. A prediction from this form of non-random mating is that stabilizing selection for the average trait value is likely to occur. Although we did not explicitly study stabilizing selection, our findings may have significant implications for understanding the mechanisms behind it. Future studies of phenotypic selection on quantitative floral traits with unimodal distribution in other plant taxa could confirm this prediction and contribute to a general understanding of plant diversity.

5. Conclusion

Our study applies the mating portfolio approach to plant sexual reproduction, offering the first quantitative estimates of female and male mating parameters at the plant level in a natural population. We show that variation in and frequency of a key floral trait – petal nectar spur length – can significantly influence maternal outcrossing rate, mate number, and mate distance. Furthermore, our investigation into the phenotypic traits of mating pairs provides evidence of positive frequency-dependent assortative mating, which has broader implications for our understanding of plant mating systems. This finding contributes to the growing theoretical interest in the relationship between assortative mating and frequency-dependent mating (Otto et al., 2008; Peischl and Schneiderr, 2010). Most of this work has considered simply inherited Mendelian traits and how these processes might maintain variation or lead to divergence. In contrast, our study involves quantitative trait variation, and it is hoped that our study will ‘spur’ interest in future theoretical and empirical work on this topic as well as on the mating success of individual plants through both female and male functions.

Acknowledgements

This paper is dedicated to the memory of Prof. Qing-Jun Li, who passed away on December 1, 2022. We are grateful to Spencer CH Barrett for critical reading and helpful editing of an earlier version of the manuscript, as well as to the handling editor and two anonymous reviewers for their constructive suggestions and comments on the manuscript. We also thank Dr. Hao Wang, Dr. Hui Huang and Ms. Xia Zou for field assistance; Mr. Yulin Weng for help with the extraction of plant DNA. This work was financially supported by the National Natural Science Foundation of China (32271693) and the Cultivating Plan Program for the Leader in Science and Technology of Yunnan Province (202405AC350111) to ZQZ.

CRediT authorship contribution statement

Mingliu Yang: Writing – original draft, Visualization, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation. Zhi-Qiang Zhang: Writing – review & editing, Validation, Supervision, Resources, Project administration, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.01.005.