Climatic variation may affect plant–pollinator interaction potential by shaping intra-annual patterns of flowering and foraging, thereby altering the windows of temporal overlap between partners. We examined whether spatial variation in temperature and precipitation predicted plant–pollinator interaction potential across a dryland ecosystem elevational gradient spanning desert to pine forest. Five sites representing distinct communities along a 2230 m elevational gradient in the Santa Rosa Mountains, CA, USA, were studied from February to August in 2021 and 2022, including 82 flowering plant species and 386 pollinator species. Using data on flowering phenology and pollinator visits to flowers, we tested whether temperature and precipitation explained variation in plant–pollinator interaction potential as quantified in two ways: whether pairs of species flowered and foraged in the same year and site (binary temporal co-occurrence) and the magnitude of overlap in their phenological distributions (phenological synchrony). Plants and pollinators showed reduced temporal co-occurrence and synchrony under warmer temperatures, and phenological synchrony between pollinators and flowers was negatively related to precipitation. Shorter pollinator foraging periods were associated with increased synchrony with flowering plants, whereas shorter flowering periods were associated with decreased synchrony with pollinators. Overall, warmer temperatures were associated with a lower likelihood of temporal co-occurrence and reduced synchrony, and higher precipitation was also linked to lower synchrony. These findings indicate that climatic conditions influence the timing and duration of plant–pollinator interactions in arid ecosystems and may alter the temporal structure of these communities.

Study sites

This study was conducted at five sites representing different communities along a 2230 m elevational gradient in the Santa Rosa Mountains along the eastern edge of the Coachella Valley, CA, USA (Figure 1A). The three lowest sites were located in the Boyd Deep Canyon Desert Research Center (BDCDRC), and the two upper sites in the San Bernardino National Forest (Figure 1B). The lowest elevation site, site one (210 m a.s.l.), was located in a low desert wash habitat, site two (830 m a.s.l.) in a high desert habitat, site three (1300 m a.s.l.) in a high desert-pinyon pine habitat, site four (1960 m a.s.l.) in a coastal chaparral habitat, and site five (2440 m a.s.l.) in a Jeffrey pine forest habitat.

We quantified climatic conditions using monthly temperature and precipitation values, which are based on daily data, from each of the five sites. For the lower three sites, we extracted monthly average temperature (°C) and monthly total precipitation (mm) data from the BDCDRC weather stations located within 30–1200 m of each site (https://deepcanyon.ucnrs.org/weather-data/, accessed on 1 May 2023). Because permanent weather stations were not available for sites four and five, we used the same weather variables (monthly averages for temperature (°C) and monthly totals for precipitation (mm)) extracted from 800 m grid cells from the parameter-elevation regressions on independent slopes model (PRISM) climate group database (Oregon State University, https://prism.oregonstate.edu, data created 2023), which is particularly suitable for mountain ranges. Over the 32-year period from 1990 to 2022, the climate across the elevational gradient became significantly hotter and more arid, but precipitation levels did not change (Figure S1).

During the 2 years of the study, site one was the hottest and driest community (Table S1). Site one and two, along the western edge of the Sonoran Desert, received precipitation almost entirely during the winter. Site three received the bulk of precipitation in the winter (mainly rain), with some North American monsoonal precipitation during the summer. Site four received both winter (rain and snow) and monsoonal precipitation in the form of summer rain storms. Site five received considerable winter snowfall as well as the highest amount of monsoonal rainfall (https://deepcanyon.ucnrs.org/weather-data/, accessed on 1 May 2023; PRISM, Oregon State University, https://prism.oregonstate.edu, data created 2023).

Data Collection

In both 2021 and 2022, we collected data weekly at each site, starting with the onset of flowering and stopping 1 week after all plants had ceased flowering, within our transects. Each site comprised three fixed 100 m × 10 m transects. Each transect was separated from the others by approximately 20 m. Each week, we collected data on the community composition of plants that were flowering, floral abundances and plant-pollinator interactions. To document which species were flowering, we first inspected all transects to detect which plants were in bloom along their length. This exploratory survey enabled us to locate where to place floral abundance quadrats such that all species that were flowering on that date were represented in our flower counts. To account for diurnal variation in flowering, we performed floral abundance surveys twice per day, in the morning and afternoon. We semi-haphazardly placed 1 × 1 m quadrats approximately every 10 m per transect, 1–5 m from the centre of the transect and alternating left and right sides. Varying the quadrat locations between surveys reduced the likelihood of double counting the same flowers on the same date. We counted all flowers or inflorescences for each species in the quadrats, for a total of 60 m² surveyed per site per day. As long as all species in bloom on a given date were represented in our quadrats, quadrats could contain zero flowers.

To collect data on plant-pollinator interactions, we surveyed each of the three transects for 20 min (1 h per session), four times per sampling day (twice in the morning, twice in the afternoon). This resulted in a total of 4 h of survey effort per site per day. Upon seeing an insect contacting floral reproductive parts (stigmas or anthers), we caught the pollinator and noted the flowering species it was visiting. We identified some pollinators (Apis mellifera, Xylocopa spp., some Lepidoptera) by sight, and recorded their interactions without collecting them. We recorded all visible pollinators during transect surveys. Individuals were identified in the field when possible; uncollected individuals were included in abundance tallies. Collected pollinators were netted and identified following permit conditions.

Plant and Pollinator Identification

We identified plants at the species level in the field, and confirmed plant identification with research staff at BDCDRC (Table S2). We caught four major groups of pollinators (wild bees, n = 902; flies, n = 582; wasps, n = 134; butterflies and moths, n = 61), and identified them to the species level when possible, or to morphospecies otherwise (Table S3). We identified pollinators at the Entomology Research Museum at the University of California, Riverside, using reference collections, dichotomous keys and the help of expert taxonomists. Voucher specimens were deposited in the Entomology Research Museum at the University of California, Riverside.

Quantifying Interaction Potential

For plants and pollinators to be observed interacting in our surveys, flowers and foraging pollinators must be present at the same time and site (Figure 2). Thus, flowering and foraging within the same year (i.e., temporal co-occurrence) act as a first filter on interaction potential (Figure 2A), and the magnitude of phenological overlap (i.e., synchrony) in flowering and foraging acts as a second filter (Figure 2B,C). To assess temporal co-occurrence, that is whether a plant-pollinator pair that was observed interacting in a given site in 1 year of sampling also co-occurred in that site in the other year of sampling, we used a binomial metric based on our plant-pollinator interaction dataset. When a pair of known interaction partner species were both flowering/foraging at the same site in the same year, they were considered to temporally co-occur (1); otherwise, instances where one or both species did not occur at the same site in the same year were considered non-co-occurrence (0). This metric does not include species pairs that were not known to interact in one or both sampling years. Thus, the possible site-specific values of co-occurrence for plant-pollinator species pairs in our dataset are: 1 (2021) and 0 (2022); 1 (2021) and 1 (2022); 0 (2021) and 1 (2022). We did not include species pairs that were never observed to interact in the 2 years of our sampling (0 [2021] and 0 [2022]), as these forbidden links can arise from various factors that are independent of climate, such as morphological mismatches or behavioural constraints (Jordano et al. 2003; Olesen et al. 2011). Thus, we examined co-occurrence only for species pairs that had known potential to interact.

All analyses were conducted in R 4.5.2 (R Core Team 2025). To estimate phenological synchrony between temporally co-occurring plant-pollinator species pairs (i.e., species pairs that had a co-occurrence value of 1 for a given site and year), we first imputed daily-resolution phenological distributions for each species × site × year combination using von Mises kernel density from the ‘densityFit’ function in the ‘overlap’ package (Ridout and Linkie 2009). We selected distribution sensitivity (bandwidth) using the ‘vm.kde’ function in the ‘Directional’ package (Tsagris et al. 2023). For plant species, we imputed distributions using data from the floral abundance surveys, summing the total number of flowers across the four daily surveys. To determine pollinator phenological distributions, we used the number of specimens collected (or number of individuals identified by sight) during the weekly plant-pollinator surveys, summing the number of times specimens of the same species were caught or observed within a day for each date. We scaled each phenological distribution from 0 to 1, such that the peak abundance for a given species equaled 1 (Figure 2). We then calculated phenological onset and cessation as the day that flowering or pollinator foraging activity reached 0.25 or 0.75 of the cumulative activity for flowering or pollinator foraging, respectively. We then calculated phenological duration as the difference between the two dates.

To calculate phenological synchrony, we quantified the proportion of overlap between plant and pollinator phenological distributions. We performed pairwise comparisons of the distributions of all temporally co-occurring plant-pollinator species pairs, modifying methodology from Carter et al. (2018), Stemkovski et al. (2020) and Fisogni et al. (2022). For each species pair, we quantified the shared area under the phenological distributions of the two species using the ‘pmin’ function and the ‘integrate.xy’ function in the ‘sfsmisc’ package (Maechler 2025). We then divided the shared area by the area under the distribution of each individual species (Figure 2). This yielded, for each pollinator and plant species separately, the proportion of its distribution shared with its interaction partner, which ranged from 0 to 1. In cases where interaction partners co-occurred but had non-overlapping phenological distributions (as in Figure 2B), the proportion of overlap was 0, indicating complete asynchrony.

Quantifying Aridity

For analyses related to temporal co-occurrence, we calculated the mean monthly temperature (°C) and monthly total precipitation (mm) for the 12 months preceding phenophase onset (i.e., onset of flowering or foraging) for plant and pollinator species that flowered or foraged in a given site and year. When a species was not active in a given year, we used the phenophase onset day from the year a species was active to then calculate mean monthly temperature and mean total precipitation for the 12 months preceding that day in the inactive year (i.e., if a species was inactive in 2021, we used the 2022 onset day of year to gauge when the species was likely to begin flowering/foraging and applied it to the 2021 climate data). Species had to be active in at least one of the two study years to be included in the analysis, providing a species-specific measure of aridity aligned with the expected phenophase onset even in years when the species was inactive.

For analyses related to phenological synchrony, we calculated the mean monthly temperature and monthly total precipitation for the 12 months prior to phenophase onset (i.e., onset of flowering or foraging) for each temporally co-occurring plant-pollinator species pair. This produced a unique value for each plant and pollinator species that flowered/foraged in different months and years.

Statistical Modelling

Models were fitted using the package ‘glmmTMB’ (Brooks et al. 2017). Prior to fitting each model, we scaled the predictor variables from 0 to 1 inclusive.

To determine whether temporal co-occurrence varied with climate, we used a generalised linear mixed-effects model (GLMM) with a binomial error distribution and logit-link. The full model included temporal co-occurrence (0,1) as the response variable, mean monthly temperature and mean monthly total precipitation as fixed effects, and elevation (as an ordered factor) as a random effect. To integrate temperature and precipitation at the species pair level, we averaged the relevant values for each pair of potentially co-occurring species. We then fitted nested models (mean monthly temperature only; mean monthly total precipitation only) and selected the best fitting model using ‘lrtest’ in the ‘lmtest’ package (Zeileis and Hothorn 2002). To test for multicollinearity, we measured the variance inflation factor (VIF), where values < 5 indicate low to moderate collinearity and acceptable values (Burnham and Anderson 2002). We performed uniformity, dispersion and outlier tests on the selected model using the ‘DHARMa’ package (Hartig 2025). We reported the AICc using the ‘dredge’ function in the ‘MuMIn’ package (Bartoń 2025), confidence intervals using the ‘confint’ function in the ‘stats’ package (R Core Team 2025), and marginal and conditional R² values using the ‘r2’ function in the ‘performance’ package (Lüdecke et al. 2021).

To determine how phenological synchrony from the plant or pollinator perspective varied with climate, we fitted two GLMMs with an ordered beta error distribution (i.e., allowing for 0–1 boundaries; Kubinec 2023) using a continuous value between 0 and 1 inclusive for plant or pollinator phenological synchrony as the response variable, mean monthly temperature and mean monthly total precipitation values specific to each plant or pollinator as fixed effects, and elevation (as an ordered factor) as a random effect. Species pair identity was included as a random effect for the pollinator model but not the plant model as the latter was overspecified. We performed the same model selection and residual diagnostic tests for these models. Examining synchrony from each partner's perspective accounts for possible inherent differences in the duration of phenophases for plants vs. pollinators. For example, a synchrony value of 1 for a pollinator that foraged for only 4 weeks is more likely than a synchrony value of 1 for a plant that flowered for 8 weeks (Figure 2C).

To determine how phenological synchrony from the plant or pollinator perspective varied with flowering or foraging duration, respectively, we fitted GLMMs with an ordered beta error distribution (Kubinec 2023). Each model used a continuous value between 0 and 1 inclusive for phenological synchrony as the response variable, phenological duration as a fixed effect and elevation (as an ordered factor) as a random effect. We again fitted separate models for plant and pollinator perspectives, selected the best-fitting model and performed residual diagnostic tests on the selected models.