Climate warming changes the phenology of many species. When interacting organisms respond differently, climate change may disrupt their interactions and affect the stability of ecosystems. Here, we used GBIF occurrence records to examine phenology trends in plants and their associated insect pollinators in Germany since the 1980s. We found strong phenological advances in plants, but differences in the extent of shifts among pollinator groups. The temporal trends in plant and insect phenologies were generally associated with interannual temperature variation, and thus likely driven by climate change. When examining the synchrony of species-level plant-pollinator interactions, their temporal trends differed among pollinator groups. Overall, plant-pollinator interactions become more synchronized, mainly because the phenology of plants, which historically lagged behind that of the pollinators, responded more strongly to climate change. However, if the observed trends continue, many interactions may become more asynchronous again in the future. Our study suggests that climate change affects the phenologies of both plants and insects, and that it also influences the synchrony of plant-pollinator interactions.

The given datasets contain estimated species phenology shifts for time and temperature. The way the slopes were estimated is given below.

We worked with occurrence records of plants and insects available from the GBIF database (GBIF.org 2021c, 2021a, 2021b, 2021d). For the plants, we restricted ourselves to species covered by the BiolFlor database of plant traits (Klotz et al. 2002), because we originally intended to classify plants by their level of pollinator dependence – an idea we later abandoned. For the insects we restricted ourselves to beetles (Coleoptera), flies (Diptera), bees (Hymenoptera) as well as butterflies and moths (Lepidoptera), as these groups contain most insect pollinators (Kevan and Baker 1983). We used the R package rgbif  (Chamberlain and Boettiger 2017) to download all available records of the above taxa from GBIF. Our basic criteria for including records were that they originated from Germany, and that their basis of record (as defined in GBIF) was either a living specimen (e.g., a captured insect), a human observation (i.e., an observation of a species made without collecting it), just an observation (i.e., when the exact type of observation was not clear), or a preserved specimen (e.g., an herbarium record or a collected specimen). If names of plant species were not accepted names, we used the R package taxsize (Chamberlain and Szöcs 2013) to check the names against the GBIF backbone taxonomy and determine the actual accepted name.

Prior to the data analyses, we subjected the data to several steps of quality control. First, we removed all records from before 1980 as these turned out to be too inconsistent, with few records per year and large gaps due to consecutive years without records. We also removed the records from 2021 as the year had not been complete at the time of our analysis. Second, we removed all records from the days of year (DOY) 1, 95, 121, 163, 164, 166 and 181, and in particular DOY 365 and 366 from the National Museum of Natural History in Luxembourg because the high numbers of records on these days indicated that either records without collecting date had been assigned these by default, or the dates were used by BioBlitz events where very large numbers of records are taken on a specific day of the year. Including these data would have strongly biased the intra-annual distributions of our records. Finally, we removed some records for which no elevation or temperature data could be obtained (see below).

To ensure reasonable coverage of the studied time interval, we then restricted the records to species which had at least 10 occurrence records in every decade covered (with the year 2020 included in the last decade).

After these data curation steps, we maintained just above 12 million occurrence records that covered altogether 1,764 species, with 11.4 million records of 1,438 plant species, around 590,000 records of 207 species of butterflies and moths, some 76,000 records of 20 bee species and 30,000 records of 22 fly species, and almost 25,000 records of 77 species of beetles (Table 1). There were large differences between plants and insects not only in the numbers of records but also in their temporal distribution across the studied period (Figure S 1). While plants had relatively even record numbers across years, the insect groups, in particular flies and bees, were strongly underrepresented in the earlier decades, and record numbers increased rapidly in the last 20 years, probably due to the advent of platforms like iNaturalist.org and naturgucker.de, which allow recording of species occurrences by citizen naturalists, and which made up most of the insect occurrence data for Germany in GBIF.

Temperature and elevation data, and individual interactions

Besides the main phenology data from GBIF records, we obtained several other data sets required for our analyses. To test for associations with climate, we used temperature data from the Climate Research Unit (CRU, https://crudata.uea.ac.uk) at the University of East Anglia, specifically the Time-Series dataset version 4.05 (University of East Anglia Climatic Research Unit et al. 2021), which contains gridded temperature data from 1901-2020 at a resolution of 0.5° longitude by 0.5° latitude. From this dataset we extracted the monthly mean temperatures and averaged them to obtain the annual mean temperatures at the sites of occurrence records. To be able to control for elevation at the locations of occurrence records, we used elevation data at a 90 m resolution from the NASA Shuttle Radar Topography Mission (SRTM), obtained from the SRTM 90m DEM Digital Elevation Database (Jarvis et al. 2021) and accessed through the raster package (Hijmans 2020) in R.

Data analysis

All data wrangling and analysis was done in R (R Core Team 2008). Before analysing phenology data, we examined patterns of climate change in Germany through a linear model that regressed the annual mean temperature values at the collection sites (= the corresponding 0.5° × 0.5° grid cells) over time.

To understand phenology changes in plants and insects, we first estimated the phenological shifts in each taxonomic group (i.e., plants, beetles, flies, bees, and butterflies/moths) as the slopes of linear regressions linking activity, i.e., the DOY of a record to its year of observation. We estimated taxonomic group-specific phenological shifts in two linear mixed-effect models: one that estimated shifts over time and one that related phenology variation to temperature. Both models included the DOY of a record as the response variable, and the latitude, longitude, and elevation of a record as fixed effects. The temporal-change model additionally included the year of a record as a fixed effect and as a random effect across species; the temperature-change model instead included the annual mean temperature at the site of a record as a fixed effect and as a random effect across species. We used the lme4 package (Bates et al. 2015) in R to fit these models, and assessed model fits by visually inspecting the relationships between residuals and fitted values, and between residuals and covariates (Supplementary diagnostic plots). As the random effects from the models agreed well enough with more complex generalized additive models, we considered our linear model reasonably robust.

Literature

Literature

Bates, Douglas; Mächler, Martin; Bolker, Ben; Walker, Steve (2015): Fitting Linear Mixed-Effects Models Using lme4. In J. Stat. Soft. 67 (1). DOI: 10.18637/jss.v067.i01.

Chamberlain, Scott A.; Boettiger, Carl (2017): R Python, and Ruby clients for GBIF species occurrence data. DOI: 10.7287/peerj.preprints.3304v1.

Chamberlain, Scott A.; Szöcs, Eduard (2013): taxize: taxonomic search and retrieval in R. In F1000Research 2, p. 191. DOI: 10.12688/f1000research.2-191.v2.

GBIF.org (2021a): Occurrence Download. DOI: 10.15468/dl.z9wz76.

GBIF.org (2021b): Occurrence Download. DOI: 10.15468/dl.3anj3a.

GBIF.org (2021c): Occurrence Download. DOI: 10.15468/dl.t2x9cj.

GBIF.org (2021d): Occurrence Download. DOI: 10.15468/dl.hb2fgs.

Hijmans, Robert J. (2020): raster: Geographic Data Analysis and Modeling. Available online at https://CRAN.R-project.org/package=raster.

Jarvis, Andy; Reuter, Hannes I.; Nelson, Andy; Guevara, E. (2021): Hole-filled seamless SRTM data V4. International Centre for Tropical Agriculture (CIAT). Available online at https://srtm.csi.cgiar.org, updated on 9/7/2021.

Kevan, P. G.; Baker, H G (1983): Insects as Flower Visitors and Pollinators. In Annu. Rev. Entomol. 28 (1), pp. 407–453. DOI: 10.1146/annurev.en.28.010183.002203.

Klotz, Stefan; Kühn, Ingolf; Durka, Walter (2002): BIOLFLOR – Eine Datenbank zu Biologisch-Ökologischen Merkmalen der Gefäßpflanzen in Deutschland. In : Schriftenreihe für Vegetationskunde, vol. 38, pp. 1–333.

R Core Team (2008): R: A language and environment for statistical computing. Vienna, Austria. Available online at http://www.r-project.org./.

University of East Anglia Climatic Research Unit; Harris, I. C.; Jones, P. D.; Osborn, T. (2021): CRU TS4.05. Climatic Research Unit (CRU) Time-Series (TS) version 4.05 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901- Dec. 2020). NERC EDS Centre for Environmental Data Analysis. Available online at https://catalogue.ceda.ac.uk/uuid/c26a65020a5e4b80b20018f148556681, checked on 9/1/2021.

Data description

The data consists of two files of comma separated values with header rows, one containing the random effects of the time shifts model and one containing the random effects of the temperature shifts model.

The column names are the same across both files:

species = Genus and species epithet for each species for which phenology shift was estimated
group = name of the random variable the parameters were estimate for (all "species")
intercept = estimated random intercept
slope = estimated random slope
intercept_std_err = standard error of the intercept estimate
slope_std_err = standard error of the slope estimate
main_var = main fixed variable of the model, either "year" or "temp"