P. caeruleum L. (Jacob’s ladder; Polemoniaceae) is a perennial plant distributed in the temperate zone of the Northern Hemisphere. Poland represents the southernmost point of its distribution range, with most populations distributed in the northeast of the country (Rutkowski, 2000). This region is characterised by continental climate with lower mean annual temperatures compared to other regions of the country.

P. caeruleum produces a corymbose inflorescence with a few to over a dozen simultaneously opening, fragrant flowers with radial symmetry. The blue corolla is composed of five petals, surrounding five stamens with orange anthers and the pistil terminating in a three-lobed stigma. Most of the flowers are hermaphrodite and protandrous; however, overlapping sexual phases are common (Zych et al., 2013). The male phase is shorter than the female phase (1.7 ± 0.9 versus 2.0 ± 0.8 days; Zych et al., 2013). Previous studies on three P. caeruleum populations, which were also included in the present study, revealed that the Polish populations of this species are characterised by different reproductive systems, two of which are self-incompatible and one is self-compatible (Zych et al., 2013; Ostrowiecka et al., 2017). According to Pigott (Pigott, 1958), individual plants live for at least 10 years and reproduce by seeds; vegetative reproduction does not occur.

Insects visiting P. caeruleum flowers are attracted to pollen, nectar, and sweet scent. The study species is characterised by a generalist pollination system – flowers are usually visited by a broad range of insects, mainly social bees, flies, beetles, and butterflies (Zych et al., 2013; Ostrowiecka et al., 2017). P. caeruleum is a pollinator-dependent species, as pollinators are the key factor determining reproductive success (Ostrowiecka et al., 2017); consequently, lack of pollinators may lead to pollen limitation.

Land drainage, climate change, and low abundance and diversity of pollinators negatively affect P. caeruleum populations. Consequently, the number of plants tends to decrease, which can lead to the complete disappearance of populations (Rutkowski, 2000; Zych et al., 2013). Additionally, high temperatures and periodic droughts, which are particularly detrimental to seedling survival (Pigott, 1958), have become common in Poland in recent years (Kundzewicz and Matczak, 2012). As a result of a significant decline in populations, P. caeruleum is included in the Polish Red List of Plants (VU category; Kaźmierczakowa et al., 2016).

The experiments involved 15 P. caeruleum populations of various sizes (both in terms of the number of individuals and area occupied; Figure 1) and comprising different types of plant communities distributed across the country. In the present study, we included all Polish populations of P. caeruleum, known from the literature, which were within the walking distance (<2 km) from roads accessible by a 4WD car. The population size was defined as the number of flowering shoots estimated every year at the peak of the flowering season (10 to 15,000 individuals per population; Supplementary Table S1).

Populations occupied different habitats: majority occurred in floristically rich wet meadows (BBR, BIA, BOB, CZL, DRO, KCZ, KLE, MAL, ORZ, WPN, and ZED), two occurred in forests or in their immediate vicinity (KOP and SIE), one occurred in a fen (ROS), and one occurred on a peninsula surrounded by sedgelands and a lake (SPN; Supplementary Table S1). Of these, three populations were located in national parks, where human activities are restricted (KOP, SPN, and WPN).

To assess the extent to which pollen limits seed production in P. caeruleum and spatiotemporal variations in this phenomenon, if any, we performed supplementary pollination experiments in all populations. During 2014–2018, at the peak of flowering (June to July in Poland), we randomly selected and marked approximately 40–50 flowers with receptive (open) stigmas from different individuals in each population. There were fewer experimental flowers in smaller populations (<100 flowering shoots). Half of the marked flowers were left for open pollination, and the rest were hand-pollinated with pollen obtained from another individual growing at a minimum distance of 3 m. Hand pollination was performed by brushing the dehisced anthers against the stigma using tweezers. The procedure was continued until the stigma was covered with pollen.

At 4–5 weeks after pollination, mature fruits were collected, and seeds were counted under a stereomicroscope. Empty seeds (aborted seeds) were excluded from the analysis. Depending on the year and population, we noted the different degree of fruit damage caused by herbivores. In some populations, we could collect almost all experimental fruits, while in others, only 30% of the fruits were undamaged. Damaged fruits and seeds were excluded from the analysis.

The pollen limitation (PL) index – a measure of the magnitude of pollen limitation – was calculated using the following formula:

where Ps is the number of seeds from the hand-pollination treatment, Po is the number of seeds from the control treatment, and Pmax is the larger of the two values (Baskin and Baskin, 2018).

The PL index was calculated for each population in a given year; thus, the values are always expressed at the population level. The PL index should range from zero, when both treatments produce the same number of seeds, to one, when natural pollination (control) does not produce seeds (complete pollen limitation). In the present study, however, in 23 of the 56 cases, the PL index was below zero, and we decided to include the negative values in further analysis, as we assumed that raising only the negative values to zero will reduce the reliability of the data.

To assess the assemblage composition of insects that visited P. caeruleum populations, we performed observations of insect activity on flowers. During 2015–2018, we recorded insect activity simultaneously with the experiments set up to assess pollen limitation. To determine the frequency of insect visitation to flowers, we applied the method described by Zych et al. (2013). We randomly selected a patch of flowering plants, with two to five shoots in full bloom, and recorded insects visiting flowers for 15 min using a digital camera (HC-VX870; Panasonic Corp.) set on a tripod at ~1.0–1.5 m from the plants; if possible, we attempted to record a different patch each time. In each year, we obtained 12 recordings of insect activity for at least 2 days (180 min⁻¹·population⁻¹·year⁻¹) in each population. Observations were recorded from 10.00 to 16.00 h (peak pollinator activity) under appropriate weather conditions; we avoided recording on rainy and windy days. In the laboratory, each insect recorded visiting a flower was counted and classified into the following groups: (i) honeybees (Apis mellifera), (ii) bumblebees (Bombus spp.), (iii) solitary bees, (iv) hoverflies (Syrphidae), (v) other flies (other Diptera), (vi) butterflies (Lepidoptera), (vii) beetles (Coleoptera), and (viii) others (e.g., wasps and members of orders Neuroptera, Hemiptera, and Orthoptera). We only considered insects that touched the flower sex organs as potential pollinators. The frequency of insect visitation was calculated per census (15 min) relative to the number of inflorescences in a recorded patch.

To determine the effects of meteorological conditions on seed production and the PL index, we collected meteorological data, including mean monthly precipitation and ambient temperature, from the nearest weather station (Institute of Meteorology and Water Management).¹ We included meteorological data from May to July of the same year in which the seed experiment was conducted. Additionally, we included meteorological data from August to October of the previous year to determine the effect of weather on the condition of plants and accumulation of resources. For instance, favourable conditions for the growth of perennials in autumn may delay the initiation of cold hardiness, which triggers the storage of assimilates in roots (McKenzie and McLean, 1980).

We used R 4.0.3 for data analyses. First, data were assessed for normality (Shapiro–Wilk test) and transformed when necessary and possible. To determine the differences in seed production and frequency of insect visitation to flowers across populations and years, we applied two-way ordinal ANOVA using the “clm” function in the “ordinal” package, followed by the “Anova” function in the “car” package, because the assumption of data normality for parametric tests was not met, even when the data were appropriately transformed.

To determine the significance of differences in seed set between hand-pollinated and naturally pollinated flowers among populations in a given year (i.e., populations experiencing pollen limitation), we applied a nonparametric Mann–Whitney U-test.

Prior to running the models, we checked for multicollinearity among the variables. As the assumption of data normality was not met for most variables, we performed Spearman’s pairwise correlation analysis, and excluded variables that were strongly correlated (r ≥ |0.7|) (Hosmer and Lemeshow, 2000) followed by variance inflation factor (VIF) stepwise selection, and excluded variables with VIF ≥ 5 (Zuur et al., 2010). Finally, variables were centred and scaled to adopt a similar range of values to improve model performance.

Generalised linear mixed model (GLMM) was used to evaluate the effects of the frequency of visitation by specific insect groups and meteorological variables on seed production in control treatment. In the model log-transformed population size, the mean frequency of visits by specific insect groups, and temperature and precipitation in selected months were included as the fixed effects and population was included as the random effect. The effect of year was not included in the model, due to near-zero variance and negligible role of the random effect of the year in relation to the random effect of the population.

During the selection of the variables to GLMM, characterizing meteorological conditions, we avoided rejecting variables describing precipitation and temperature during months in which the seed experiments were set, as we found them more influential. Finally, all variables were centred and scaled to adopt a similar range of values to improve model performance. As the data on seed set included many zero values and were characterised by negative binomial distribution, we applied zero-inflated negative binomial GLMM using the “glmmTMB” function in the “glmmTMB” package (Brooks et al., 2017), as it provided a better fit than the other models [as evidenced by a lower Akaike information criterion (AIC) value] and prevented overdispersion. We selected the best models based on AIC. The models developed during the analyses are presented in Supplementary Data (Supplementary Table S2).

Furthermore, we applied a linear mixed model using the “lmer” function in the “lme4” package to analyse the possible predictors affecting the PL index, as it was characterised by normal distribution. In two models, the calculated values of the PL index for each population in each year were combined with the (i) average frequency of insect groups visiting flowers and (ii) average monthly meteorological variables for each month, including population and year as the random effects. As zero models in both cases were characterised by the lowest AIC value, we averaged the models. For this purpose, we compared the AICc values of the candidate models, and averaged those including factors with ΔAICc <2 (compared with the zero model, characterised by the lowest AICc). Finally, only one factor was included in each model: frequency of visits by honeybees and temperature in October (ΔAICc values of the candidate models are presented in Supplementary Tables S3 and S4).

To determine whether the PL index is shaped by the population size and the frequency of insect visits, we performed linear regression analysis using the “lm” function in the “stats” package. We also performed linear regression analysis to determine whether the frequency of insect visitation was affected by population size, after transforming the dependent variable (frequency of insect visitation) to a normal distribution by adding 0.5% to each score, followed by log-transformation.

The number of seeds produced differed among populations and years in both control (χ² = 803.70, df = 14, p < 0.001 and χ² = 798.79, df = 4, p < 0.001, respectively) and hand-pollinated flowers (χ² = 809.26, df = 14, p < 0.001 and χ² = 758.62, df = 4, p < 0.001, respectively).

The mean number of seeds per fruit was, respectively, 11.9 ± 8.0 and 12.7 ± 7.9 in the control and hand-pollinated treatments, indicating 6.3% increase with the latter, although this difference was not significant. However, seed production significantly differed among populations and years, indicating that pollen limitation was temporally present in at least three populations (Figure 2). Pollen limitation primarily affected small populations; however, this phenomenon was also observed in one large population (MAL). Additionally, in the three affected populations, control flowers produced significantly more seeds than hand-pollinated flowers (Figure 2). The degree of pollen limitation varied across populations and years, ranging from −0.77 to 0.96, with mean of 0.04. Moreover, 23 of the 56 PL index values in P. caeruleum populations were below zero.

The overall insect visitation frequency was 6.9 ± 8.9 visits per census (15 min) per inflorescence. The frequency of insect visits significantly differed among populations (χ² = 151.676, df = 14, p < 0.001) and years (χ² = 52.133, df = 3, p < 0.001; Figure 3).

During 143 h of recording, we observed 13,553 insect visits on P. caeruleum flowers, and honeybees and bumblebees accounted for over 70 and 17% of all visits, respectively. Dipterans, lepidopterans, and coleopterans accounted for only 5, 1.2, and 2% of all visits, respectively (Figure 4).

Honeybees were the most frequent flower visitors among the recorded insects, and some populations were visited almost exclusively by honeybees (Figure 4). Meanwhile, in some populations, bumblebees were the predominant visitors in some years. The populations located in national parks (KOP, SPN, and WPN) were visited by more diverse insect groups.

The highest mean monthly temperatures were recorded in populations in the south of the country (BBR: 15.2°C; MAL: 15.0°C; and CZL: 14.9°C) and in the population in the northeast (DRO: 15.0°C). In contrast, the lowest temperatures were recorded in populations in the north of the country (WPN: 13.9°C; KCZ: 14.0°C; and ZED: 14.2°C).

Regarding precipitation, less pronounced trends were noted. The highest precipitation was recorded in two population in the north (SPN: 63.5 mm and BOB: 61.8 mm) and one population in the northeast (ZED: 61.0 mm), and the lowest precipitation was recorded in populations in the northeast (SIE: 44.2 mm; KLE: 46.3 mm; and KOP: 49.8 mm). Across years, plants that produced seeds in 2018 were exposed to the highest temperatures, whereas plants that produced seeds in 2015 were exposed to the lowest precipitation.

The number of flowering shoots varied among populations, and while this number remained relatively stable in most of the populations, there were marked fluctuations in some populations (e.g., ZED, MAL, KCZ, DRO, and BOB) during the study period (Supplementary Table S1).

The GLMM testing the response of number of seeds in control treatment showed a positive effect of population size (GLMM, p = 0.014), the visitation frequency of hoverflies (GLMM, p = 0.02), and average precipitation in June (GLMM = 0.046). On the other hand, the visitation frequency of honeybees, precipitation in August of the previous year, and average temperature in September of the previous year (GLMM, p = 0.041, p = 0.042, p = 0.012, respectively) negatively affected seed production in control treatment (Table 1). The difference between the value of R2 marginal (describing the proportion of variance explained by the fixed variables), and R2 conditional (describing the proportion of variance explained by both the fixed and random variables), indicate that the share of variance explained by the random variable (populations) in the model was moderate.

The linear mixed models demonstrated that the PL index was negatively affected by the visitation frequency of honeybees (Table 2A) and positively affected by ambient temperature in October (Table 2B). In both models, among the random variables, population explained more variation than year; however, in general, the random variables had little effect. The value of the R² marginal was much lower than that of R² conditional, indicating that the random variables (years and populations) in the model explained most of the variation.

Overall, the PL index was negatively affected by population size (F-statistic = 4.54, DF = 53, p = 0.038; Figure 5) and the frequency of insect visitation (F-statistic = 5.93, DF = 45, p = 0.019; Figure 6). However, the frequency of insect visitation was not associated with population size (F-statistic = 0.002, DF = 47, p = 0.962).

The present study demonstrated that P. caeruleum is pollen limited, but only to a small extent, at least in the Polish range of its distribution. Moreover, we observed spatiotemporal variations in pollen limitation, but with no obvious trends. These results are consistent with previous reports on pollen limitation in the Polish populations of this species, indicating that some populations suffer from pollen limitation; however, this phenomenon depends on population and year (Zych et al., 2013; Ostrowiecka et al., 2017). In the present study, the PL index markedly varied among populations and years, ranging from −0.77 to 0.96, and over 40% of the values were below zero. Surprisingly, in some populations, hand-pollinated flowers produced fewer seeds than naturally pollinated flowers. According to Young and Young (1992), this result can be explained by several facts, such as stigma damage during hand pollination, negative effects of high pollen density on the stigma on pollen tube growth in the style (stigma clogging), lower diversity of pollen donors in hand-pollinated flowers, and pollen removal by insects from the stigma of hand-pollinated flowers. Although we cannot rule out any of the above possibilities, our observations of insect activity suggest that insects visiting flowers treat the stigmas covered with pollen as a food source.

In three of the analysed populations (in one of them during 2 years), control flowers produced significantly fewer seeds than hand-pollinated flowers. Among these populations, two were small (BBR and SPN, with 10 and 90 flowering shoots, respectively, in the corresponding year) and one was large (MAL, with 6,000 flowering shoots in the corresponding year). The presence of pollen limitation in the large population may be the result of a sudden increase in the number of flowering individuals (probably as an effect of mowing) and consequent increase in the visitation frequency of insects, particularly honeybees, collecting large amounts of pollen, which may not left enough to pollination of other flowers (for details see next subsection “Flower visitors and their effects on seed set and pollen limitation”).

In contrast, in three other populations, control flowers produced significantly more seeds than hand-pollinated flowers. Among these populations, two were small (SIE and BIA, with 20 and 100 flowering shoots, respectively, in the corresponding year) and one was medium (ROS, with 120 flowering shoots in the corresponding year). In these three populations, insects likely treated the profusely hand-pollinated stigma as a food source. Additionally, apiaries are present close to the two small populations; and in BIA alone, honeybees were almost exclusive flower visitors and the frequency of their visitation was also high. The medium population was surrounded by a fen, with few other co-flowering plants; thus, stigmas covered with pollen may have been an attractive food source.

The magnitude of pollen limitation depends, among other factors, on the plant mating system; as such, self-compatible plants and generalists are less prone to pollen limitation than self-incompatible plants and specialists (Ashman et al., 2004; Knight et al., 2005; Wolowski et al., 2014). The generalist pollination system of P. caeruleum may reduce the probability of pollen limitation in populations, as has been shown in other plant species (Knight et al., 2005; Gómez et al., 2010; Wolowski et al., 2014). The P. caeruleum populations analysed in the present study are overall characterised by the presence of mixed mating systems (unpublished data), which rather does not predispose them to a high degree of pollen limitation.

The pollination system of P. caeruleum is described as generalist. In this plant, the cup-shaped flowers, which offer nectar and pollen, are easily accessible to a wide spectrum of insects. In the present study, we recorded visits by diverse groups of pollinators. In 3 of the 15 study populations, the assemblages of insects visiting flowers have been previously analysed. According to Zych et al. (2013), bumblebees were the most frequent and efficient visitors in KLE. Meanwhile, Ostrowiecka et al. (2017) observed that while KOP was visited by a diverse assemblage of insects, ZED was mainly visited by honeybees. In the present study, we observed a similar assemblage composition of insect visitors in these three populations in all seasons examined. However, this was not true for the rest of the studied populations. As such, the composition of insect assemblages varied across years, indicating spatiotemporal variations in the potential pollinator assemblages of P. caeruleum.

Considering data across all seasons and populations, honeybees were the most frequent visitors of P. caeruleum flowers, and together with bumblebees, they accounted for nearly 88% of total insect visits. Other groups of insects, although present in most populations, accounted for over 12% of all visits to flowers. Among these, flies (including hoverflies) were the largest contributors, whereas solitary bees, beetles, butterflies, and other insects contributed to a lesser extent. These trends support the previous assumption of Zych et al. (2013) that in terms of the pollination system, regardless of the high apparent generalisation, the realised generalisation is rather low among the P. caeruleum populations studied. Several groups of insects willingly used the food sources offered by P. caeruleum flowers, but only when honeybees and bumblebees were sparse.

The high proportion of honeybees (70.5%) among insects visiting P. caeruleum flowers is a consequence of the widespread introduction of colonies of these insects by beekeepers. Populations distant to settlements, such as those located in national parks (KOP, SPN, and WNP), were characterised by a lower frequency of honeybees visits and a more diverse composition of flower visitors assemblage. This result suggests the negative effect of honeybees on the diversity of other pollinators visiting flowers, corroborating previous reports (Aslan et al., 2016; Mallinger et al., 2017). Honeybees are highly competitive and may disrupt interactions of other native pollinators with plants (Goulson, 2003; Thomson, 2016; Valido et al., 2019; Milner et al., 2020; Angelella et al., 2021). The authors’ personal observations are that in populations characterised by a high visitation frequency of honeybees, pollen was removed from anthers very efficiently and rapidly, often resulting in difficulties in finding pollen for the hand-pollination of experimental flowers. A negligible amount of available pollen reduces the effectiveness of pollen transmission by insects and may thus reduce the seed set.

The GLMM showed that the visitation frequency of honeybees negatively affected seed production in open-pollinated flowers. However, despite the negative effect of honeybees on seed number in control treatment, among other insect groups visiting flowers, only honeybee visitation frequency showed a significant negative effect on the PL index, suggesting an important role of these insects as pollinators of P. caeruleum. This result is confirmed by previous findings that honeybees, in addition to bumblebees, are the key pollinator groups based on the pollen load carried by individuals and the frequency of their visitation (Zych et al., 2013). We assumed that the influence of honeybees on pollination of P. caeruleum may change with the number of insects visiting the population. In higher densities, those insects remove pollen from anthers very efficiently, and probably at some point, due to lack of available pollen to pollination, it may negatively influence the reproductive success of P. caeruleum. It is also possible that the observed negative impact of honeybees on seed set in control treatment is related to the collection of pollen grains from the stigmas by those insects, as evidenced by the fact that the PL index often took negative values. However, the negative PL index values may arise from the delivery of pollen of inadequate quality to the stigma (e.g., incompatible) of hand-pollinated flowers, or stigma damage during the hand-pollination by the experimenters. Consistently, honeybees appear to play a dual role as insects that on the one hand may in fact reduce pollen limitation, as reported previously (González-Varo et al., 2009; Tscheulin and Petanidou, 2011), while on the other hand may be a competitive species that collect and remove pollen from the population and in consequence preclude efficient pollination of plants.

Furthermore, hoverflies visits, although not recorded in all populations or years, positively affected seed set in control treatment. Hoverflies can have a strong positive effect on seed set (Vance et al., 2004), and in some plant species characterised by open flowers, these insects are more efficient pollinators than bumblebees (Fontaine et al., 2006). Moreover, hoverflies can carry pollen over longer distances than bees (Lysenkov, 2009), which may increase the probability of cross-pollination with the pollen of unrelated individuals, resulting in the production of more seeds of better quality, in P. caeruleum. Additionally, our observations indicate that flies, including hoverflies, visited fewer flowers within an inflorescence than bees, which may limit the deposition of pollen from the same individual.

Unlike the previous results of Zych et al. (2013) that reported bumblebees as the most efficient pollinators of P. caeruleum, this study did not reveal that higher bumblebee visitation rate influence the reproductive success of this species. However, Zych et al. (2013) studied only a single population, where these insects accounted for most of the visits. In contrast, our results showed that insect assemblages vary spatially and temporally and may be shaped by beekeeping activities, to some degree. Apart from human activity, including beekeeping, also other factors, such as the immediate surroundings, type of habitat and population size of P. caeruleum, as well as the presence and abundance of other groups of insects, and meteorological conditions in a given year, probably influenced the spatial and temporal diversity of insects visiting P. caeruleum flowers (Somme et al., 2014; Paajanen and Cronk, 2020; Rohde and Pilliod, 2021).

We observed that bees and, to a lesser degree, flies were the frequent visitors of P. caeruleum flowers. Likewise, a previous study (Zych et al., 2013) showed that these insect groups carried considerable P. caeruleum pollen loads. As the mean number of ovules in P. caeruleum flowers is rather low (n = 29), all morphogroups are expected to serve as efficient pollinators.

In P. caeruleum, the lack of specialisation in terms of pollination biology and the spatiotemporal variations in the composition of pollinator assemblages may result in high divergence in flower traits related to pollination. Our previous study on the chemistry of P. caeruleum nectar showed high variability in these traits among populations (Ryniewicz et al., 2020). We also observed some divergence in terms of the anatomy of P. caeruleum flowers, such as differences in flower size or the arrangement of reproductive parts. We assume that the morphology of P. caeruleum flowers shows certain adaptations to specific groups of pollinators; however, because of the high instability of the environment, these adaptative traits are not fixed, and such a variability of floral traits can be interpreted as “adaptive wandering” (Zych et al., 2019). In addition, the selection of randomly appearing features in P. caeruleum populations may be disrupted due to the presence of honeybees, which depend on human activity; this may affect the presence and diversity of other pollinators. Unpredictable visitation of flowers is a common feature among generalist plant species (Ghazoul, 2005; Ollerton et al., 2007; Wolowski et al., 2014), and this may lead to the development of self-compatibility (Arista et al., 2017). Local adaptations and variable degree of self-compatibility in P. caeruleum populations warrant further research.

The present study revealed that meteorological factors significantly affect seed production in control treatment and the PL index. Seed set was negatively affected by average monthly temperature in September of the previous season, when plants start to enter a state of dormancy and developed cold hardiness. In addition, seed set was negatively affected by average monthly precipitation in August of the previous season but positively affected by average monthly precipitation in June of the same season, when the field experiments were conducted. In case of meteorological factors which were analysed for the impact on the PL index, we observed only the influence of average monthly temperature in October, which increased the PL index.

P. caeruleum is a boreal plant species that prefers moist soils and does not tolerate well high temperatures, and a seasonal drought, during the vegetative season. All these weather extremes are becoming increasingly frequent and common in Poland (Kundzewicz and Matczak, 2012), posing a severe threat to the development of P. caeruleum populations. Our results indicate that unfavourable meteorological conditions during late summer and autumn the strongest affected the reproductive success and pollen limitation in populations of P. caeruleum (in the next season). Higher temperatures during autumn delay the development of cold hardiness (McKenzie and McLean, 1980) and can render plants weak, predisposing them to damage during winter.

This result supports the notion that temperature is an important factor shaping P. caeruleum distribution. On the other hand, the positive influence of average monthly precipitation in June, during P. caeruleum flowering and seed set indicate that the accessibility of soil water is especially important in this crucial time for plant reproduction.

However, there is some uncertainty in the context of the effect of monthly meteorological conditions. For instance, we observed differences in phenology among populations and years. Occasionally, the distance of some populations from the nearest weather station was over 40 km. Moreover, in addition to reproduction, favourable external conditions are essential for other activities, such as vegetative growth or processes related to extending survival (Sletvold and Gren, 2011; Sletvold and Ågren, 2015).

Finally, population size affected some of the measured parameters. P. caeruleum populations were selected at the beginning of the project considering their sizes: five small, five medium, and five large; however, over the course of the study, the number of individuals in most populations changed significantly (Supplementary Table S1).

The effects of plant population size on seed production and pollen limitation through pollinator activity are well-known. For instance, this may be attributed to the lower attractiveness of small populations characterised by lower resources and/or deposition of higher amounts of heterospecific pollen on the stigma, particularly among generalist plant species (Jennersten and Nilsson, 1993; Ågren, 1996; Ashman, 2005; Knight et al., 2005; Aizen and Harder, 2007; Morales and Traveset, 2008). Although our results did not confirm the effect of population size on the visitation frequency of all insects, its effects on the seed set in control treatment and PL index were evident. Moreover, the visitation frequency of all insects negatively affected the PL index. Of note, honeybees and bumblebees, which were responsible for nearly 90% of all visits to P. caeruleum flowers, are characterised by flower constancy. The decrease in plant population and availability of limited food resources may be perceived by these insects as a severe threat, forcing the major pollinators to switch to different food sources. This may further worsen pollen limitation, ultimately leading to the complete disappearance of plant population.

Despite the wide array of factors tested in our study, still many of them, which are potentially important in determining degree of pollen limitation were not analysed in the present study, including uneven distribution of resources in individual parts of the inflorescence (Strelin and Aizen, 2018), species richness of co-flowering plants, interspecific competition to attract pollinators (Bell et al., 2005; Vamosi et al., 2006), resource availability (Haig and Westoby, 1988; Ashman et al., 2004), and habitat traits (Somme et al., 2014). Moreover, reproduction costs may differ between populations (Sletvold and Gren, 2011), which is also likely to be the case for P. caeruleum, as evidenced by the effect of populations on seed production as per our GLMM.

Long-term studies on pollen limitation are paramount for identifying trends in population dynamics and establishing the direction of change. Our results underscore the risk of further decline in small populations of P. caeruleum, which are subjected to more severe pollen limitation. This may decrease the genetic diversity in populations, ultimately leading to their complete disappearance. The present study may serve as the scientific basis for the conservation of P. caeruleum in its Polish range.