1. Introduction

Recent studies highlight the crucial role of animal pollination in angiosperms, with about 90% of angiosperms depending on animals for pollination (Ollerton et al., 2011; Tong et al., 2023). For a given plant species, pollinators often transfer both conspecific and heterospecific pollen, leading to heterospecific pollen deposition on stigmas. Heterospecific pollen deposition, observed in 245 species across 52 families on five continents (Arceo-Gómez et al., 2019), can clog stigmas (Caruso and Alfaro, 2000), usurp ovules (Burgess et al., 2008) and potentially reduce female reproductive success (Campbell and Motten, 1985). The negative effects of heterospecific pollen deposition can be mitigated, as conspecific pollen has a germination and fertilization advantage over heterospecific pollen, which creates a reproductive barrier (Alarcón and Campbell, 2000; Campbell et al., 2009), although one that is not always effective (Morales and Traveset, 2008; Jakobsson et al., 2009). For example, in Phacelia parryi, heterospecific pollen deposition was found to negatively affect pollen tube growth and reduce seed set, especially when heterospecific pollen was applied to the stigmas before conspecific pollen (Bruckman and Campbell, 2016). The impact of heterospecific pollen on reproduction varies, and has been shown to be influenced by multiple factors, including the amount of conspecific pollen and heterospecific pollen, and the specific combination of pollen donor and recipient species (Briggs et al., 2016; Lanuza et al., 2021; Fang et al., 2023). In addition, heterospecific pollen deposition may have positive effects on reproduction (see Tur et al., 2016; Bi et al., 2024). These adverse and beneficial effects are important for understanding the dynamics of interspecific interactions and their evolutionary strategies (Morales and Traveset, 2008; Sargent and Ackerly, 2008; Moreira-Hernández and Muchhala, 2019).

Heterospecific pollen transfer is common among co-flowering species, but both the amount of heterospecific pollen receipt and the proportion of heterospecific pollen (HP%, calculated as heterospecific pollen divided by total pollen amount, represents the degree of heterospecific pollen interference) can vary greatly (Montgomery and Rathcke, 2012; Fang et al., 2019; Wei et al., 2021). For instance, HP% is less than 0.1% in several species, whereas in others it exceeds 50% (Montgomery and Rathcke, 2012; Fang et al., 2019). This variation in heterospecific pollen receipt among species suggests two primary evolutionary strategies to mitigate the impacts of heterospecific pollen deposition on female reproductive success: (1) pre-pollination avoidance or reduction of heterospecific pollen deposition, and (2) post-pollination tolerance of heterospecific pollen deposition (Ashman and Arceo-Gómez, 2013; Arceo-Gómez et al., 2016). Plant species may diverge in floral phenology, develop pollinator specialization, promote flower-visiting constancy, and use multiple pollen placement sites on pollinators to minimize heterospecific pollen deposition (Armbruster et al., 1994; Muchhala and Thomson, 2012; Huang and Shi, 2013). In contrast, other species may have evolved tolerance to heterospecific pollen to minimize its deleterious effects (Arceo-Gómez et al., 2016; Hao et al., 2023). A necessary prerequisite for the tolerance-avoidance mechanism is spatial–temporal stability of pollen load on stigmas. If HP% varies significantly among species but remains relatively stable over years and communities, this indicates that the heterospecific pollen deposition might exert different persistent effects on different species, potentially leading to the tolerance or avoidance mechanisms. Conversely, if variation in HP% is greater between years and communities than between species, this suggests that pollen deposition is more influenced by community species composition and dynamics. In this case, HP% would not be a stable attribute of the species, making the differentiation between tolerance and avoidance mechanisms less likely (Arceo-Gómez, 2021).

Patterns of heterospecific pollen deposition within species can create complex mosaics of heterospecific pollen transfer and receipt, potentially leading to varied adaptive landscapes if these population differences in heterospecific pollen receipt persist over time (Arceo-Gómez et al., 2016; Arceo-Gómez, 2021). However, there is a scarcity of studies evaluating the extent of these spatial–temporal variations in heterospecific receipt across multiple species, limiting our understanding of the consequences of heterospecific pollen deposition. For instance, a three-year study of 34 co-flowering species in a single subalpine meadow revealed that while most species experienced low heterospecific pollen proportions (indicative of avoidance), several species received higher quantities of heterospecific pollen, surpassing conspecific pollen (indicative of tolerance) (Fang et al., 2019). To our knowledge, this spatial variation in heterospecific pollen receipt across multiple species was largely unknown (but see Ashman and Wei, 2024), limiting our understanding of the consequences of heterospecific pollen deposition. Addressing this gap, we collected stigma samples from 19 plant species across six communities over four consecutive years. We identified and counted the pollen species deposited on stigmas and analyzed variations in conspecific pollen and heterospecific pollen deposition across species, communities, and years. Using these data, our research primarily focused on the variation in pollen receipt on stigmas among co-flowering species, as well as its spatial–temporal variation. We predicted that significant among-species variation in heterospecific pollen deposition would be observed across different years and communities. Moreover, we investigated spatial–temporal variations of HP%, to test whether the variations in HP% between species is greater than that between years and between communities, as would be required as a precondition for the avoidance-tolerance response.

2. Materials and methods 2.1. Study sites and species sampling

We conducted our study from 2019 to 2022 around Shangri-La city in Yunnan Province, southwestern China, one of the world's biodiversity hotspots (Myers et al., 2000). In this study, six subalpine meadow communities were selected, and a quadrat of about 50 m × 50 m was marked in each meadow. The minimum distance between these communities was about 4 km and the maximum distance was 25 km (Fig S1; See Table S1 for details). Every year, from July 25 to 30, we counted the composition of flowering plant species in each community. To elucidate the differences in pollen deposition among communities, we selected 19 flowering plant species that were relatively abundant in terms of flower quantity. These species belong to 12 different families, and each species bloomed in at least two different communities (Table S2). Several species were not sampled in certain communities for consecutive years due to annual fluctuations in the composition and abundance of flowering species. Each species was observed in at least two communities within each year and had a minimum of three years of repeated observations. An exception was Galeopsis bifida Boenn., which was observed in two communities in 2019 and 2020, in one community in 2022, but not in 2021 because of insufficient flowering individuals (Table S2).

On average, we sampled 12.3 plant species per community annually, with a range of 9–19 species (see Table S1). The sampling species composition among different communities varied (Table S2). The Bray–Curtis dissimilarities of species compositions between any two communities averaged 0.33 ± 0.14 (mean ± SD), ranging from 0.19 to 0.47, with a lower value indicating more similar flowering species shared between two communities. The annual collection process for each species in each community was termed as an annual sampling. In total, 281 annual samplings (2019, 70; 2020, 69; 2021, 70; 2022, 72, respectively) were conducted across six communities for 19 plant species averaging 14.8 annual samplings per species, with a range from 5 to 24.

2.2. Stigma collection and stigmatic pollen counting

To investigate the pollen deposition on stigmas, we collected stigmas from flowering plants of all communities across four consecutive years. The collection adhered to the following criteria: (1) following two consecutive sunny days; (2) at the end of their flowering stage and after achieving saturated pollination; (3) thirty stigmas per species for each community; (4) typically one, and rarely two stigmas from each individual, regardless of whether the plant had a single flower or an inflorescence with multiple flowers; (5) exclusion of individuals from the same patch; (6) completion of the collection in each community within two consecutive days. We further evaluated the sampling sufficiency using rarefaction curves. The rarefaction curves of most species tended to plateau at the collected sample size of 30 stigmas, indicating that the sample size was sufficient to adequately represent heterospecific pollen abundance (Fig S2). Stigmas were preserved in clean microcentrifuge tubes with 70% ethanol. In the lab, each stigma was processed on a glass slide to identify and count conspecific and heterospecific pollen grains. Pollen that detached during storage was centrifuged (7000 rpm, 10 min). The bottom 40 μl was then transferred to a slide for examination. A 400× magnification microscope was used to identify the conspecific (CP) and heterospecific (HP) pollen on each stigma based on size, shape, color, and surface features. The proportion of heterospecific pollen (HP%) was calculated as HP/(CP + HP) × 100%.

2.3. Data analysis

Before analysis, we excluded a small number of stigmas (n = 7, 0.1% of total) with no pollen grains (either conspecific or heterospecific). We calculated the average pollen receipt (conspecific pollen, heterospecific pollen and HP%) for all species within each community across different years, as well as the cumulative totals for all communities each year and for each community over years. Person correlations were used to compare the heterospecific pollen proportions for each species across different communities in different years. To evaluate the among-species dynamic of pollen deposition, the coefficient of variation (CV) of total stigmatic pollen, conspecific pollen and heterospecific pollen loads and heterospecific pollen proportion were calculated for each species based on stigma level. The CV (SD/mean) was calculated based on HP% and conspecific pollen among the 30 stigmas for each species within each community and year, treating each stigma as one estimate. Due to the high number of zero values in the original data for heterospecific pollen and HP% (3962, 47.0% of all), we used zero-inflated Generalized Linear Models (GLMs) with binomial distribution (HP% as dependent variable) and GLMs with Gaussian distribution (CV of HP% as dependent variable) to test whether the variation of HP% and its CV are driven primarily by conspecific pollen or heterospecific pollen. To avoid the potential impact of zero-inflation in heterospecific pollen deposition, we reanalyzed the non-zero dataset using a GLM model with binomial distribution (HP% as dependent variable). Considering the positively skewed distribution of the CV data, we applied a log transformation to CV. To evaluate the effects of year, species and community (all treated as fixed factors, including the interactions) on patterns of stigmatic pollen receipt, zero-inflated GLMs with Poisson distribution (conspecific pollen, heterospecific pollen as dependent variable) and binomial distribution (HP% as dependent variable) were conducted. For each response variable, there are four candidate models: (1) an additive model (only three fixed factors); (2) Model 1 plus the interaction between Species and Year; (3) Model 1 plus the interaction between Species and Community; and (4) a full model including interactions between Species and Year, and between Species and Community (Table S3). We used Akaike's Information Criterion (AIC) and Bayesian Information Criterion (BIC) with the lowest values to determine the best model. Consequently, we selected Model 4, as the optimal model for conspecific pollen and heterospecific pollen, and selected Model 1 for HP% (Table S3). We used the 'glm' and 'glmmTMB' packages in R v.4.1.2 to perform GLMs analyses. The model comparison was conducted using the 'glm' and 'anova' function from the package in R v.4.1.2. Further, the R package 'glmm.hp' was used to compute the relative importance of fixed-effect variables and interaction terms based on hierarchical variation portioning (Lai et al., 2022, 2023).

Additionally, we calculated the average HP% for each annual sampling over four years and six communities and used hierarchical cluster analysis (average linkage, Euclidean distance) to categorize them into clusters, focusing on the relationship of HP% levels from most to least similar. To perform the hierarchical clustering of HP%, we constructed a matrix where the rows represented the 19 co-flowering species, and the columns represented the 24 combinations of six communities over four years. Each element in the matrix contained the average HP% in a specific annual sampling. This matrix could assess similarities and differences in HP% among different species, and also compare variations in HP% within the same species across different spatial (communities) and temporal (years) scales. We present bar charts showing the average HP% for specific annual samplings, focusing on those species whose classifications varied in our clustering analysis across different community–year combinations.

To examine whether the HP% has a systematic pattern, we constructed a phylogenetic tree relating all study species based on sequences of the Internal Transcribed Spacer (ITS) gene, sourced from GenBank (https://www.ncbi.nlm.nih.gov/genbank/, accessed November 2023). Due to the absence of the Aster oreophilus sequence, its close relative species, Aster souliei, was used as a substitute (Table S2). The phylogenetic tree was constructed using the maximum likelihood method with IQ-tree v.2.2.6, based on sequence alignment results obtained with MAFFT v7.450. Blomberg's K value is used to assess the phylogenetic signal of the HP%. The Blomberg's K value is used to assess whether the distribution of HP% across species is associated with phylogenetic signals or driven by ecological pressures. A K value less than 1 indicates a weak association of HP% with phylogeny, suggesting that variations in HP% are more influenced by environmental factors or pollinator behavior. This analysis was conducted using the 'phylosig' function from the 'phytools' package in R (Blomberg et al., 2003). The rarefaction curves based on the heterospecific pollen quantity on the stigmas of each species were used to evaluate the relationship between sample size and pollen richness deposited on the stigmas for each community and year. The rarefaction curves included 95% confidence intervals and were generated using 500 randomizations with replacement. Each stigma was considered as a sampling unit, and the total number of stigmas represented the overall sampling effort. This analysis was conducted using the 'ggiNEXT' function from the 'iNEXT' package in R (Chao et al., 2014).

Other analyses were performed using R v.4.1.2 and involved packages such as 'pvclust', 'ggplot2', 'car', 'Matrix' and 'tidyverse'. These comprehensive approaches allowed us to assess the dynamics of pollen deposition and its variations across species and communities.

3. Results 3.1. Variation in pollen receipt

In total, we collected 584, 631 pollen grains from 8432 observed stigmas, among which 562, 632 (96.2%) were conspecific and 21, 999 (3.8%) heterospecific (Table 1). Each stigma, on average had 69.3 pollen grains, of which 66.7 (92.3%) were conspecific and 2.6 (7.7%) were heterospecific pollen. The stigmatic pollen loads did not follow similar patterns across communities. Specifically, the average conspecific pollen load, heterospecific pollen load, and HP% differed significantly across communities and years (all p < 0.001; Fig S3; Table 2), except that the difference in HP% between communities was not significant (z = 1.12, p = 0.263; Table 2c). The amount and proportions of pollen loads varied considerably among the 19 studied species (all p < 0.001; Fig. 1; Table 2). Cynoglossum amabile Stapf & Drummond received the largest amount of conspecific pollen, with mean values of 349.0, 194.1, 349.0 and 260.8 from 2019 to 2022, respectively. By contrast, the lowest amount of conspecific pollen was received by Spenceria ramalana Trimen in 2019 (4.8; mean number) and 2020 (4.5), Dipsacus asper Wall. in 2021 (5.4), and Ajuga forrestii Diels in 2022 (6.9). Mean numbers of heterospecific pollen grains on stigmas varied significantly among species, ranging from 0.1 (2021, CZV, Tibetia himalaica (Baker) H.P. Tsui) to 39.4 (2022, SABG, Silene yunnanensis Franch.).

The optimal models explained 97.4%, 69.3% of variations in conspecific pollen load and heterospecific pollen load, respectively (Table 2). For conspecific pollen load, Model 4 showed that the most important factors in explaining variation were Species (z = −31.27, p < 0.001, accounting for 29.79% of the variance) and the interaction between Species and Year (z = 31.04, p < 0.001, accounting for 29.79% of the variance; Table 2a). For heterospecific pollen load, Species was the most important factor (z = 20.43; p < 0.001), explaining 29.79% of the variation (Table 2b). Community (z = 5.83, p < 0.001, 5.92% of the variance) and Year (z = −3.86, p < 0.001, 4.51% of the variance) had small but significant impacts on heterospecific pollen load variation (Table 2b). Additionally, the interaction between Species and Year had a modest but significant effect on heterospecific pollen load variation (z = 8.67, p < 0.001, 18.47% of the variance) (Table 2b). Our findings demonstrate that while conspecific pollen and heterospecific pollen deposition on stigmas vary among species and are affected by spatial (community) and temporal (year) factors, the variation is predominantly driven by species identity. This highlights the importance of species-specific characteristics in determining pollen deposition patterns.

3.2. Variation in heterospecific pollen proportion

The average proportion of pollen grains that were heterospecific varied greatly among species (Fig. 1), ranging from 0.1% (2021, CZV, Tibetia himalaica) to 41.8% (2022, SABG, Silene yunnanensis). The optimal models (Model 1) explained 97.4% of variations in HP%. The most important effects on HP% were Species, accounting for 88.98% of the variation (z = 120.81, p < 0.001; Table 2c). In addition, Year had a small but significant effect on HP% variation (z = −4.48, p < 0.001, accounting for 0.04% of the variance), whereas Community did not significantly affect HP% variation (z = 1.12, p = 0.263, accounting for 0.04% of the variance; Table 2c). These results indicate that there are significant differences in HP% among species and the most important predictor of this variation is Species identity, which accounts for the majority of the variance in HP%. Furthermore, the results of zero-inflated GLMs showed that HP% is primarily driven by heterospecific pollen load (z = 71.17, p < 0.001, accounting for 34.21% of the variance) with a lesser but still significant contribution from conspecific pollen load (z = −36.22, p < 0.001, 8.87% of the variance; Table S4a). Analysis of the non-zero dataset showed similar results (Table S4b).

3.3. Preconditions for natural selection

HP% for each species was positively related in all pairs of years and across all pairwise community comparisons (Figs. 2 and S4), which indicates that among-species variation in HP% persists and is similar across years and communities. Our separate hierarchical cluster analysis based on annual samplings showed that most species were consistently classified in distinct HP% groups. Two plant species (S. yunnanensis and D. asper) were consistently classified into a high HP% group (Fig. 3a). Eight plant species were categorized into a medium HP% group in different communities across years. Seven species were consistently classified into a low HP% group (Fig. 3a). Two species (Pedicularis cephalantha Franch. and Salvia przewalskii) were classified into different groups across communities and years. P. cephalantha was classified in the medium HP% group in four annual samplings in 2019 and 2022, but in the low HP% group in nine other annual samplings (Fig. 3b). S. przewalskii was classified in the medium HP% group in six annual samplings during 2020, 2021 and 2022, but in the high HP% group in fifteen other annual samplings (Fig. 3c).

Furthermore, the average of HP% across all annual samplings exhibited a low phylogenetic signal across the entire tree (K = 0.111; Fig S5), and ranged from 0.095 to 0.137 for each year. This suggests that the phylogenetic signal of HP% was weaker than expected, indicating that HP% is a highly plastic trait in evolutionary terms. It is more likely to be influenced by environmental factors or other ecological drivers rather than being constrained by phylogenetic relationships. Additionally, the CV of HP% was lower in the high HP% groups compared to the low HP% group (Fig. 4; Wilcoxon test, p < 0.001). This implies that species with higher HP% display more consistent HP% across individuals, indicating uniformity in heterospecific pollen exposure, as opposed to those in lower HP% groups which exhibited greater variability. For species with high HP%, heterospecific pollen was deposited on 71.8%, 68.1%, 71.3% and 69.4% of stigmas of species with high HP% in the four years, respectively. In comparison, for species with low HP%, it was deposited on of 21.8%, 16.6%, 18.8% and 22.7% stigmas in the four years, respectively. Additionally, GLMs revealed a significant effect of both conspecific pollen load (z = −16.40, p < 0.001) and heterospecific pollen load (z = 7.49, p < 0.001) on the CV of HP%, with heterospecific pollen (accounting for 44.78% of the variance) explaining more of the variation than conspecific pollen (10.62% of the variance; Table S4c).

4. Discussion

Our study focused on stigmatic pollen loads, the quantities of conspecific and heterospecific pollen grains, illuminating potential interactions among plant species within a generalized pollination system. Our key finding is that variation in proportions of heterospecific pollen (HP%) is spatially and temporally consistent among co-flowering species. Species with high HP% showed greater heterospecific pollen deposition annually within any given community, whereas those with low HP% showed the opposite trend. Despite significant interannual and intercommunity variations in the amounts of conspecific and heterospecific pollen at the species level, HP% remained relatively stable, fluctuating within a narrow range.

4.1. Spatial–temporal variations in pollen receipt

Theoretical models suggest that the characteristics of plant-pollinator networks develop randomly, implying that heterospecific pollen received by plants might be distributed in a random manner (Kallimanis et al., 2009). However, an observational study examining pollen receipt in 34 species within a meadow community revealed a temporally consistent difference among species in levels of HP% across different species over three years (Fang et al., 2019). Our findning that most species were consistently grouped into different HP% categories supports this work (Fig. 3a). GLMs models also showed that one of the most important factors to explain HP% variation is species identity. Moreover, we found that the interaction between Species and Year has a similar contribution, indicating that the inter-annual variation of HP% is not consistent across different species (Table 2c). Most species exhibited small variations in HP% between years, with several species showing significant inter-annual variation, resulting in their classification into different HP% groups. Given that the species abundance and pollinator composition within communities may change from year to year, the significant effect between Species and Year is to be expected (Arceo-Gómez and Ashman 2014; Emer et al., 2015; Tur et al., 2016; Arceo-Gómez 2021). The composition of the pollinator population has been suggested to be a critical factor that influences pollen transfer, as there is considerable variation in the size and diversity of heterospecific pollen loads (Arceo-Gómez and Ashman 2014; Kay et al., 2019; Smith et al., 2021). For instance, one study in Hawaii demonstrated that spatial variation in pollen receipt within a species primarily stems from variation in abundance of a dominant invasive pollinator (Johnson and Ashman 2019). Despite this, the between-year correlations within the same community demonstrated a positive relationship in HP%, suggesting that even with annual variations in HP%, these species tended to remain within a consistently high or low range of HP%.

Studies that have analyzed pollen deposition at one site over specific periods of time (July and August) for years have suggested that heterospecific pollen deposition is relatively stable across species (Fang et al., 2019). However, a single-site survey may not fully capture the stability of pollen deposition (Arceo-Gómez 2021). Our study expands on this by examining spatial variations in both intraspecies and interspecies pollination. We observed notable differences in species composition and richness across various communities (as detailed in Table S2). For example, one community (SABG) had the highest species richness, with less fluctuation in species composition across years. Another community (SNR) exhibited a relatively simple species composition. Additionally, sustained grazing has drastically reduced the numbers of certain species in particular years (Table S1). These community differences are also reflected in the observed pollen grain counts on stigmas. In 2019, for instance, the average conspecific pollen count for Astragalus pullus N.D. Simpson in SABG had an average conspecific pollen count of 52.3 was 52.3 in the community with the highest species richness (SABG), but much lower (27.8 ± 1.4) in a community with simple species composition (SNR). Similarly, in 2020, heterospecific pollen receipt in S. przewalskii was on average 4.7 HP grains per flower one community (JFV), but only 2.3 in another (ZCV). Although great variation was observed across communities between a pair of years, HP% exhibited a similar spatial–temporal pattern. That is, species experienced the same high or low HP% on the stigmas, between pairs of years within each community (Table S5) and between communities within each year (Table S6). Our results broaden the understanding of the generality and variations of multiple species differences in heterospecific pollen deposition, compared to previous studies that evaluated heterospecific pollen deposition for a single species or species pairs, or focused on pollen transfer within limited temporal and spatial ranges (Montgomery and Rathcke 2012; Fang and Huang 2013; Zhang et al., 2021).

4.2. Heterospecific pollen-avoidance or tolerance hypothesis

Although previous studies provide insights into the coexistence and adaptation of flowering species to interspecific pollen transfer, there has been limited data and experimental validation over the years to support the tolerance-avoidance mechanism (Ashman and Arceo-Gomez., 2013, Arceo-Gómez, 2021). A predictable outcome of this mechanism is that avoidance species have a lower HP%, tolerance species have a relatively higher HP%, and these among-species variations should be consistently detectable across different spatial and temporal scales. Only one long-term study (3 years) of a natural plant community has demonstrated that the patterns of HP% among species are similar across different years, but this study has yet to be replicated across multiple communities (Fang et al., 2019). Species that employ avoidance mechanisms may feature specialized floral traits to minimize mismatches between pollen output and stigma contact. For example, the typical floral isolation in Pedicularis species demonstrates a high degree of constancy between pollinators and plants, as well as a strategy employed by plants to reduce heterospecific pollen transfer (Huang and Shi 2013). In our study, the average HP% for three Pedicularis species was consistently low (mean, 1.8%), ranging from 0.4 to 4.4, and these low levels were maintained over the four-year survey period, and even across different communities with significant variations in species composition and abundance. This aligns with our expectations for the avoidance mechanism, where plants minimize the possibility of heterospecific pollen deposition at the pre-pollination stage.

Species that regularly receive higher levels of heterospecific pollen might develop tolerance mechanisms, with larger stigmas and longer styles potentially more inclined to develop tolerance mechanisms (Ashman and Arceo-Gómez 2013; Tong and Huang 2016; Arceo-Gómez et al., 2019). In our study, the results of hierarchical cluster analysis indicated that three plant species were classified into the high HP% group, and this was mostly consistent across different years and communities (Fig. 3a). It is important to note, one species, S. przewalskii, had greater variation of HP% in different communities. Specifically, from 2020 to 2022, the average HP% in several communities was below 12.9 %, placing the species into the intermediate group (Fig. 3c). A recent study of S. przewalskii reported that HP% substantially varied among sites, ranging from 3.4% to 51.3% (Fang et al., 2023). Although HP% varied among communities, it was still higher than the average HP% of other species within communities, indicating that significant variations among species still remained. In summary, the spatial–temporal stability of variations in HP% aligns with the expected outcomes of the avoidance-tolerance mechanism and provides a necessary, though not sufficient, precondition for the evolution or maintenance of a spatial–temporal mechanism in plants.

We also found that species with higher HP% exhibited higher consistency in heterospecific pollen receipt among individuals, as indicated by the lower CV (Fig. 4). The GLM results further support this, showing a negative correlation between heterospecific pollen and CV, implying that as heterospecific pollen increases, the variability in HP% decreases. This indicates that individuals within these species are more consistently exposed to heterospecific pollen (Table S4). In contrast, species with higher conspecific pollen deposition, such as Pedicularis species, typically employed that minimized strategy, often resulting in low HP%. Most individuals in these species received only small amounts of heterospecific pollen, although a few occasionally received larger amounts, leading to a higher CV for heterospecific pollen. This pattern may explain the positive effect between conspecific pollen and CV of heterospecific pollen. Conspecific pollen had a relatively lower explanatory power (12.11% of the variance) compared to heterospecific pollen (23.42% of the variance), which indicates that the variation in CV is more strongly driven by heterospecific pollen than by conspecific pollen. Previous studies have indicated that greater variability in heterospecific pollen receipt among individuals lowers the opportunity for natural selection to act on traits that reduce heterospecific pollen effects (Arceo-Gomez 2021). In our study, species with high heterospecific pollen loads exhibited lower variability in HP%. This stability in heterospecific receipt reduces the unpredictability of pollination events and individual differences, providing greater opportunities for natural selection. Converesely, species with high heterospecific pollen receipt are more likely to evolve tolerance mechanisms to mitigate the effects of heterospecific pollen. These results support the premise that the spatial–temporal stability of pollen receipt patterns is a prerequisite for adaptive evolution under natural selection, and also suggests that the variability in pollen receipt can influence the opportunity for natural selection, such as in the adaptive evolution to heterospecific pollen transfer.

Our study evaluated the spatial–temporal stability of heterospecific pollen deposits, but our dataset still covers only a limited spatial scale. The absense of seasonal variation has likely lead to potential homogeneity in plant species pools, which might limit the generalizability of our results. In addition, our dataset does not include pollination environment dynamics within communities, such as the composition and abundance of plant and pollinator species. This limitation restricts our evaluation of whether these dynamics significantly impact the results of stigmatic pollen deposition in plants. Future studies that include dynamic data of the pollination environment would be helpful in demonstrating the spatial–temporal stability of interspecific pollination. Experimental study, such as hand-pollination, could provide direct evidence that would allow us to assess the impact of heterospecific pollen deposition and the adaptive outcomes of different plants. For instance, a study on Clarkia xantiana showed that populations previously exposed to heterospecific pollen experienced reduced interference compared to unexposed ones, indicating that historical interactions between heterospecific pollen donor and recipient species could encourage the evolution of heterospecific pollen tolerance (Arceo-Gómez et al., 2016). This tolerance might be specific to long-term co-existing species pairs. For instance, Fang et al. (2023) found a negative correlation between heterospecific pollen reception and seed set in S. przewalskii Maxim., with seed set being more affected by rare heterospecific pollen species than common ones, which suggests varying effects of heterospecific pollen from multiple donors on species with tolerance mechanisms. However, most hand-pollination studies lack surveys of natural pollen deposition, which are essential for selecting appropriate donor and recipient species. A recent study that chose three Silene species as pollen recipients and donors based on previous surveys found that seed set was unaffected by heterospecific pollen deposition, irrespective of excess or insufficient conspecific pollen (Kai et al., 2023). This supports the tolerance hypothesis, filling a gap in experimental research on tolerance mechanisms.

5. Conclusions

Our study investigated pollen receipt and its spatial–temporal variation across 19 species within six sub-alpine meadow communities for four consecutive years. We found that variations in HP% among co-flowering species showed remarkable stability. This supports previous research that indicated heterospecific pollen receipt may be a characteristic trait of certain species (Fang et al., 2019). Notably, our findings bolster the tolerance-avoidance hypothesis, suggesting plants may evolve distinct adaptations to frequent interspecies pollination. Future research should not only investigate spatial–temporal variations in heterospecific pollen deposition but also the factors that influence these variations, including co-flowering community and pollinator community composition. Our research contributes to understanding the dynamics of interspecific pollen transfer among plant communities and strategies to mitigate its negative effects, enhancing our knowledge of the mechanisms maintaining species diversity in biodiversity hotspots.

Acknowledgments

We thank Shiyun Guo, Jiajun Wu, Suyan Ba and Chunyan Zhao for participating in the fieldwork and collection. We thank the editors and the reviewers for their valuable advice and comments that greatly improved out manuscript. This work was supported by the National Natural Science Foundation of China (No. 32071535 and 32371602), the Key Scientific Research Projects of College and University in Henan Province (No. 24ZX001) and the Henan Province Foundation for University Key Teacher (No. 2020GGJS074).

CRediT authorship contribution statement

Tao Zhang: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Conceptualization. Qiang Fang: Writing – review & editing, Writing – original draft, Supervision, Investigation, Funding acquisition, Concept-ualization.

Data availability statement

The data that support the findings of this study are openly available in the Science Data Bank at https://www.scidb.cn/en/anonymous/STdSUmJp.

Declaration of competing interest

The authors declare that they have no conflicts of interest.

Appendix A. Supplementary data

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