Historical and citizen‑reported data show shifts in bumblebee phenology over the last century in Sweden

Bumblebees are a key taxon contributing to the provision of crop pollination and ecosystem functioning. However, land use and climate change are two of the main factors caus-ing bee decline across the world. In this study, we investigated how the flight period of bumblebee spring queens has shifted over the last century in Sweden, and to what extent such shifts depended on climate change, landscape context, latitude, and the phenology of bumblebee species. We studied ten species of bumblebees and used observations from museum specimens covering 117 years from the southernmost region in Sweden (Scania), combined with citizen-reported observations during the past 20 years across Sweden. We found that the flight period of bumblebees has advanced by 5 days on average during the last 20 years across Sweden. In the agriculture-dominated region of Scania, we found that in the late 2010s bumblebee spring queen activity in simplified landscapes had advanced by on average 14 days, compared to 100 years ago. In addition, in simplified landscapes the flight period of early species was significantly earlier compared to in complex landscapes. Our results provide knowledge on the intraspecific variation of phenological traits, indicating that early species (often common species) exhibit a higher plastic response to the environment, which may facilitate adaptation to both climate and landscape changes, compared to the late species of which many are declining.


Introduction
Insect pollinators such as bumblebees (genus Bombus) contribute critically to both crop pollination (Klein et al. 2007) and ecosystem functioning (Ollerton et al. 2011). Bumblebees are in decline in many parts of the world at local or regional scales (Biesmeijer et al. 2006;Potts et al. 2010;Zattara and Aizen 2021) which could threaten the provision of pollination services and ecosystem health (Potts et al. 2010, Vanbergen andthe Insect Pollinators Initiative 2013). The persistence of populations of central-place foragers like bumblebees is determined by the availability of food (flower) resources in the vicinity of nesting sites during their period of activity (Osborne et al. 2008;Smith et al. 2014). The availability of food resources is, in turn, affected by both land use and climate change, which are two of the main drivers of bee decline (IPBES 2016;Soroye et al. 2020;Kammerer et al. 2021;Prestele et al. 2021).

Land-use change and bumblebees
Since the 1960s and the Green Revolution, agricultural management has shifted towards more intensive production to fulfil global needs of food, feed, and fibre. This intensification of agriculture has led to landscape scale simplification of land use, e.g. by increased field sizes, more monocultures, and specialization of farming into either crop or animal production, with concomitant reduction of crop diversity and loss of semi-natural habitats such as pastures and meadows, resulting in a general biodiversity loss (Donald et al. 2001;Stoate et al. 2001). Bumblebees show trait-dependent vulnerability to simplification of agricultural landscapes, such that species with above-ground nests, few workers per colony, a late phenology and long colony cycle are rarer in simple landscapes dominated by arable and annual crops, compared to complex landscapes with smaller sized arable fields, permanent pastures and other semi-natural habitats (Goulson et al. 2005;Persson et al. 2015).
As bumblebees are eusocial insects, the selection of nest sites by founder queens is crucial for colony establishment and development. Nest site selection may be influenced by landscape context and availability of resources when queens emerge (Lye et al. 2009). Landscape simplification has the potential to alter the availability of both nesting and floral resources for bumblebees (Goulson et al. 2008;Roulston and Goodell 2011;Persson and Smith 2013). The combination of lack of semi-natural habitats and a simplified crop monoculture can lead to seasonal floral resources gaps, e.g. lack of late-flowering resources, which can affect the development of the colony at different stages (Persson and Smith 2013;Rundlöf et al. 2014;Timberlake et al. 2019). Although less is known about the landscape scale availability of nesting sites, one can assume that the general lack of preferred nesting habitats, such as semi-natural grassland, non-crop field borders and edges of woodland (Osborne et al. 2008, O'Connor et al. 2017, reduce nest site availability in simplified agricultural landscapes (Persson et al. 2010).

Climate change and bumblebees
Climate change consists of gradual changes in average conditions and an increased frequency of extreme events such as droughts and heatwaves (IPCC 2021). Gradually changing conditions can lead to reduction of species distribution ranges (Kerr et al. 2015;Sirois-Delisle and Kerr 2018;Marshall et al. 2020), and disrupted plant-pollinator interactions (Memmott et al. 2007) due to phenological mismatches (Scaven and Rafferty 2013;Gérard 1 3 et al. 2020a, b;Freimuth et al. 2021). Bumblebee species respond differently to changes in climatic conditions (e.g., Williams et al. 2007; Bartomeus et al. 2013;Pradervand et al. 2014;Rasmont et al. 2015;Herbertsson et al. 2021), and studies suggest there are limits to their adaptive capacity to track climate change (Kerr et al. 2015, Martinet et al. 2020a. Importantly, land-use-driven changes to landscape structure coupled with climate change will affect the availability of suitable habitats for bumblebees, with differences among species (Prestele et al. 2021). In particular, declines are expected for habitat specialists with poor dispersal capacity (Bartomeus et al. 2013;Wood et al. 2019), while warm adapted species may benefit if suitable habitat is available (Bartomeus et al 2013;Herbertsson et al. 2021). This highlights the need for studies combining land use and climate change, to assess possibilities to simultaneously mitigate/alleviate the negative effects of both factors (Titeux et al. 2016;Prestele et al. 2021).
Bumblebee queens can use spring air or soil temperature as a cue to emerge from hibernation (Alford 1969;Goodwin 1995;Kudo and Cooper 2019). Increased (spring) temperatures can therefore directly lead to advanced phenology of bumblebee populations (Herbertsson et al. 2021). In addition, bumblebee populations may track flower resources, and can thus respond indirectly to climate warming via phenological shifts in resource availability (Ogilvie et al. 2017). Phenological mismatches in plant-pollinator networks can occur when the response of bee activity and flowering phenology to climate warming are not equally strong or use different cues (Renner and Zohner 2018;Kudo and Cooper 2019;Gérard et al. 2020a, b). Such mismatches have been observed at the species level (e.g. Osmia spp. in Germany (Kehrberger and Holzschuh 2019)) and at the community level (e.g. wild bees in Belgium (Duchenne et al. 2020a) and in Europe (Duchenne et al. 2020b). There is, however, variation in the phenological responses among bee species: while some species advance their activity period, others delay it, and some show no change (Duchenne et al. 2020b). A better mechanistic understanding of this heterogeneity and its consequences for pollinator communities could be obtained by, for example, analysing how species' traits affect their response to temperature and land use changes (Bartomeus et al. 2013;De Palma et al. 2015). For example, one trait related to bumblebee phenology where different strategies exist is the date of queen emergence after hibernation: while some species emerge early during spring, others emerge later in the summer (Benton et al. 2006). Phenology of bumblebee queen emergence could indicate how temporally well-matched the colony development and flower resources are (Memmott et al. 2007).

Phenology of Swedish bumblebees
In Sweden, bumblebee flight period ranges from March to September, with the main period of queen emergence coinciding with the flowering of Sallow (Salix caprea) in spring (Mossberg and Cederberg 2012). After emergence, each queen searches for a nesting site, lays and incubates eggs, forages and provides for the first brood. When the first workers emerge, the queen will remain in the nest to lay eggs while workers start provisioning the nest and care for the brood. Later in the growing season, as the colony reaches a critical size, the queen will produce sexual individuals (males and new queens) (Mossberg and Cederberg 2012). Sweden is a large country with a large gradient of climate conditions across latitudes (SMHI 2021). This means that the bumblebee community structure differs along the latitudinal gradient, but also that the emergence of the same species is expected to differ across latitudes. Relatively to the pre-industrial era, temperature has increased more at northern latitudes compared to the south (Christensen et al 2019). One could therefore expect a stronger response and a more rapid phenological shift towards earlier queen flight periods at northern latitudes. The large climate gradient in Sweden thus makes a perfect case study to explore the response of bumblebee emergence to difference regimes of climate warming. On the other hand, the region of Scania in the very south of Sweden has large variation in agricultural intensity (from specialized cereal production to extensive mixed farming systems, Persson et al 2010), but a relatively homogeneous climate, which allows the study of phenological shifts due to agricultural-induced landscape change along similar trajectories of climate warming.

Long-term monitoring using museum specimens and citizen science observations
Long-term monitoring of pollinators is essential to assess both population trends and the efficiency of interventions to safeguard pollinators (IPBES 2016). Insect natural history collections (NHC) provide a historical baseline, both to identify pollinator declines and to study global change impacts (Kharouba et al. 2019;Bartomeus and Dicks 2019). For example, NHC have become an important source to track morphological shifts in response to climate (MacLean et al. 2019) and land-use changes (Gérard et al. 2020a, b). In addition, the increasing interest of citizens to report wildlife observations and the creation of transdisciplinary online platforms for wildlife reporting have allowed the compilation of large data sets to monitor current pollinator trends (Birkin and Goulson 2015;Domroese and Johnson 2017;Bartomeus and Dicks 2019;Garratt et al. 2019;Zattara and Aizen 2021). The Swedish Species Observation System (hereafter Artportalen) is one such system. Citizen science data arising from systematic sampling can be of very high quality while data arising from opportunistic reporting has limits, e.g. concerning measures of abundance and over-sampling of rare species (e.g. Bird et al. 2014).
Here, we aimed to investigate long-term shifts in the phenology of bumblebee spring queen flight period in Sweden and to what extent any such shifts depend on climate change, latitude, and agricultural-induced landscape simplification. In addition, we explored if bumblebee phenology (early vs. late queen emergence) may explain any shifts in flight period. To this end, we used observations of bumblebee queens from museum collections from the southernmost region in Sweden (Scania) covering 117 years, and observations from Artportalen covering the whole country during the past 20 years.
We ask the following questions and hypothesize the following answers: Hypothesis 1 How has the flight period of bumblebee queens changed over time during the last two decades across Sweden, and how does this depend on latitude? We expect a general rapid phenological shift towards earlier queen flight periods, and a stronger response at northern latitudes, where the relative increase in temperature is higher compared to the south.
Hypothesis 2 How has the timing of the flight period of bumblebee spring queens shifted over the past 117 years in the southernmost region of Sweden (Scania), and are any such trends similar in simple and complex agricultural landscapes? We expect that bumblebee queen activity has advanced and that this trend is more pronounced in simple than in complex landscapes, due to competition for nesting and flower resources and a flower resource gap later in the season in simple landscapes.

3
Hypothesis 3 Has the flight period of bumblebee spring queens across Sweden (national dataset) and in Scania (regional dataset) shifted over time and with mean annual temperature, as a proxy for climate change? We expect that the flight period has advanced both over time and with increased temperatures, and that the effects of time and temperature are indeed interchangeable due to temperature being the driving force.
Hypothesis 4 Are shifts in queen flight period explained by the species-specific bumblebee phenological trait queen emergence time (early vs. late) in Sweden and Scania, respectively? We expect a general shift in both early and late-emerging species towards earlier flight period in response to climate change and observed phenology shift in plants (earlier flowering). In addition, we expect queens of late-emerging species to have shifted towards earlier dates particularly in simple agricultural landscapes, where the availability of lateflowering resources (e.g. clover-rich leys, semi-natural grasslands) has declined over time.
Hypothesis 5 Is there a difference in queen activity period between species that have declined versus those that increased in relative abundance over the past 120 years? We expect that declining species are generally late emerging ones, while stable or increasing species have an earlier activity period.

Dataset
We selected 10 species of bumblebees to cover species with both early and late phenology ( Table 1). Most of these species are readily identified in the field, although females of B. terrestris and B. lucorum can be difficult to separate. Freshly emerging spring queens are likely easier to properly identify as they are large and their fur is not yet worn. However, B. lucorum cannot be separated from two other species (B. cryptarum and B. magnus) without the use of molecular methods, and these three species are therefore referred to as the B. lucorum-complex (Williams et al. 2012). Thus, we expect records of B. lucorum to also contain the other two cryptic species. The data used in this study consisted both of specimens from museum collections and citizen science observations: 1088 museum specimens collected between 1899 and 2016 in Scania, Sweden, of which 836 were assessed as spring queens (see below), and 9279 citizen science observations from Artportalen covering all of Sweden from 2000 to 2019, of which 4384 were spring queens. Both datasets contained data on all 10 bumblebee species. See Fig. 1 for the distribution of observations of the two datasets and Appendix Table S1 for the number of records per species.
The museum specimens were obtained from the entomological collection of the Biological Museum, Lund University, acronym MZLU (Fägerström 2020). From each specimen, we gathered information on the site, collection date, and name of the collector. The collection site (usually the name of a village and/or a specific part of the village) was used to obtain geographic coordinates for the analysis. We only included museum specimens for which we had information about both site and day of collection. The museum data mainly contained specimens collected in Scania between 1899 and 2016.
Citizen observations were downloaded from the Artportalen website, the Swedish Species Information Centre (SLU Artdatabanken 2020). Artportalen is an infrastructure hosted by the Swedish Agricultural University (SLU) which contains both opportunistic species observations and observations from systematic sampling, reported by the public (citizens)  and by conservation officers at public authorities and from the private sector (SLU Artdatabanken 2020). Reports are regularly checked and verified by experts, especially concerning rare species and observations outside known distribution ranges (pers. comm. Niklas Johansson, ArtDatabanken, 01/10/2021). On the 10th and 11th of March 2020, we searched the Artportalen database of occurrence records. We downloaded data for all 10 species selected (Table 1), for all possible years in Artportalen, and for observations during the flight period of bumblebees in Sweden. We thus used the following filters: [Species = (the scientific name for each bumblebee species); in Period, Months = "March to September" for "all possible years"; in Sighting properties, Gender = "Female" (to select queens, and distinguish from workers)]. We obtained 10 296 observations for all Sweden downloaded as an excel file and filtered for records where there was complete information of day of observation and location. In some cases, observations were reported within a time period of several days. In those cases, we calculated the mean date for the reported period. If the reported time period exceeded a month, we excluded the observation from the analysis.

Agricultural production areas
To analyse if shifts in bumblebee phenology in the southernmost region of Sweden relate to agricultural landscape context, we used the large agricultural production areas (Stora Produktionsområde in Swedish) defined by The Swedish Board of Agriculture (Statistics Sweden 2014). These consider the large-scale conditions for farming, including soils, topography, and climate, and are thus suitable as a proxy for general agricultural intensity and landscape heterogeneity. The production areas largely correlate with latitude in Sweden: from more productive and intensive agricultural land in the south towards lower productivity in the north. In Scania, three production areas are represented: arable plains with simplified landscape structure highly dominated by arable fields (area 1), mixed landscapes (area 2), and forested mixed landscapes (area 5). The latter two contain smaller fields, more grass and clover leys, permanent pastures, and other semi-natural habitats (Persson et al. 2010;Persson and Smith 2013). Following Herbertsson et al. (2021), we used these production areas as a proxy for historical land-use change in Scania, because since the 1950s, the proportion of arable land and field size have continued to increase in the simplified landscapes (area 1), whilst agriculture has partly been abandoned and replaced by forest in the mixed landscapes (area 2 and 5) (Ihse 1995). We pooled records from areas 2 and 5 (hereafter complex landscapes) as only a few specimens were sampled in area 5 and explored the difference in bumblebee phenological shifts between simplified landscapes (area 1) and complex landscapes (areas 2 and 5) (Fig. 1). For each data point collected, we matched the site coordinates of queen observations with the agricultural production areas using QGIS 3.16.5.

Temperature data
As a proxy for the extent of climate change over the sampling period of the citizen observations dataset we obtained the mean annual temperature and the mean spring annual temperature (March, April, May) at the national scale in Sweden from 1899 to 2019 from the Swedish Meteorological and Hydrological Institute (SMHI 2021). For the museum dataset covering only Scania, the southernmost part of Sweden, we instead used a more detailed resolution for regional temperature data from the Climatic Research Unit (CRU) 1 3 time-series temperature data (Jones and Harris 2012), covering the years 1901-2011 (the available time-scale on the dataGuru data repository). For the museum dataset this meant a slight reduction of the sampling period and thus sample size (46 specimens lost from before 1901, and 1 specimen lost from after 2011). Since no specimens were recorded in 1901, 2010 or 2011, the final dataset for the museum data covers the years 1902-2009.

Separation between spring queens and daughter queens
To study shifts in the flight period of emerging queens during spring, it was necessary to exclude any daughter queens produced during summer from analyses. This was done by using an algorithm to identify activity peaks of the two groups of queens in the datasets, corresponding to the peak of a normal distribution, and thus obtain a breakpoint between peaks. Because of the large geographic area covered by the Citizen observations, the date of the breakpoint separating the flight periods of spring queens and daughter queens vary greatly, especially along the latitudinal gradient. We therefore first defined climatic regions with homogeneous annual temperatures based on spatially resolved temperature data for all sampling points. We calculated the mean annual temperature for 1980-2011 (the available time-scale) for each sampling point, based on the spatial interpolation (using the kriging method; Cressie & Wikle 2015) of the Climatic Research Unit (CRU) time-series temperature data (Jones and Harris 2012). We performed the spatial interpolation using the variogram() function (gstat R package) (Pebesma 2004), with an exponential variogram. We then clustered all the sampling points based on their mean annual temperature using one-dimensional kernel-density estimation, with the pdfCluster() function (pdfCluster R package) (Azzalini and Menardi 2013). We obtained seven groups of sampling points with homogeneous annual temperatures. Because the northern regions of the country contained fewer observations than the southern ones, we grouped the four most northern regions, to form a total of four groups for further analyses. This division between the most northern climatic region and the three southern ones matched the biogeographical border called the 'Limes Norrlandicus' (Biologiska Norrlandsgränsen, in Swedish). In each climatic zone, we then clustered sampling points for all species based on their collection day (calendar day), using the pdfCluster() function. We found at least two queen flight periods in each climatic region, and we selected the date between the first and second flight periods as the breakpoint between the spring (nest-founding) and summer (daughter) queen flight periods (Fig. S1). Only specimens assigned as spring bumblebee queens were used in the analyses on phenological shifts.
For the Museum specimens, we used the same procedure to calculate the breakpoint date between the spring-emerging queens and summer daughter queens, but it was applied to the whole Museum dataset as Scania is one homogeneous climatic region (Fig. S2).

Statistical analyses
We performed our analysis in R, version 4.2.1 (R Core Team 2022), using linear mixedeffects models. To analyse how the flight period of bumblebee queens shifted over time, we assumed that the specimen observation date (citizen science dataset) or collection date (museum dataset) were proxies for queen flight period. Both the Museum and Artportalen datasets cover large spatial extents and extended periods, combining data from different observers who did not follow a standardized sampling procedure. As sampling effort was not homogeneous in space and time, we had to account for sampling bias in both datasets, 1 3 a recurring difficulty when dealing with historical data (Rollin et al 2020;Isaac & Pocock 2015;Carvalheiro et al 2013). We adapted the methodology developed by Carvalheiro et al. (2013), by drawing a 10 × 10 km grid over Sweden. In each grid cell, we estimated the yearly observation date of bumblebee species by calculating the median observation dates of all specimens of a given species for each year. We conducted all statistical analyses at the grid cell scale, with median observation date per species as a response variable, and observation year and geographical location as main explanatory variables. To correct for sampling bias among grid cells, we weighted each grid cells by its yearly sampling effort (total number of collected or observed bumblebee queens per species) using the weight argument from the lmer() function (lme4 R package), so that grid cells with more samples a given year had a higher weight in the analysis than poorly sampled locations. All statistical models also included a grid cell identifier as a random effect, to consider spatial and temporal dependency in the data (grid cells were revisited across years and neighbouring grid cells are spatially dependant) (Fig. S3). For the Scania Museum data, we followed the same procedure, but analyses were conducted within a grid of 5 × 5 km. We used this finer spatial resolution in order to accurately qualify the landscape context of each sampling site (simple or complex landscapes).
To evaluate shifts in flight periods in the Scania Museum data, we used linear mixedeffect model (LMM) with median observation date as a response variable, and all possible two-way interactions between landscape type, year (collection year) and species as explanatory variables. For the Sweden Citizen-collected data, we ran another LMM including all possible two-way interactions between latitude, year (collection year) and species. All continuous variables were scaled prior to analysis. Non-significant interactions were removed following stepwise backward elimination to improve model parsimony.
To explore whether shifts in flight periods differed between early and late emerging species, we used the same structure as previous models but replaced the species variable by species-specific queen emergence trait (categorical variable: early vs. late). We included all possible two-way interactions between emergence trait, landscape type and year (for museum data); all the two-way between emergence trait, latitude, year (for the Citizencollected data).
To assess whether species-specific trends showed similar responses to time and changing temperatures we ran an additional model, where we substituted the explanatory variable collection year by temperature. This model thus included all possible two-way interactions between mean spring temperature, landscape type and species for the museum dataset; and mean spring temperature, latitude and species for the citizen-collected dataset. Model comparisons are found in the Appendix.
We excluded Bombus distinguendus from the analyses at the national scale (citizencollected dataset) because of the very small sample size of spring queens (N = 6, Fig. S4). We report model results using an analysis of variance ANOVA (type II) Satterthwaite's method, and used post-hoc Tukey tests with the function lsmeans() (emmeans R package) (Lenth et al. 2018) to test for differences in collection day between the different levels of categorical explanatory variables. Model predictions and prediction errors were calculated using the function effect() (effects R package) (Fox and Weisberg 2018), and we illustrated the model fit results with the function plot_model() (sjPlot R package) (Lüdecke 2018).
We checked homoscedasticity and normality of model residuals by visual inspection of diagnostic plots and did not detect any major issues. The correlation coefficients between the explanatory variables were |r|< 0.7, meaning no substantial collinearity between the variables (Dormann et al. 2013). All model residuals were checked for spatial autocorrelation with a Mantel test using the function mantel.rtest() (ade4 R package) (Dray and 1 3 Dufour 2007) (see Table S2), and none of the models showed spatial autocorrelation. In models where mean spring temperature was substituted for year, we checked for temporal autocorrelation of model residuals. For the regional scale, museum data model, we found no signs of temporal autocorrelation (see Aucorrelation function (ACF) graphs Fig. S5a). However, for the national scale, citizen-collected data model, the ACF graphs showed signs of moderate temporal autocorrelation (Fig. S5b), which we handled by including year as a random factor, nested in the grid cell identifier (Fig. S5c).
To investigate whether species activity period or intraspecific variation in activity period were associated with species susceptibility to decline over the past 120 years, we combined both datasets of spring queens (museum and Citizen-collected data), covering a time frame from 1899 to 2019. We measured change in species relative abundance compared to historical relative abundance, by dividing each species' current relative abundance by their historic relative abundance, following Austin and Dunlap (2019) and Hatfield et al. (2014), based on the IUCN criteria for red list assessments (IUCN 2022). To estimate differences of the current population size relative to historical levels, we split the time range in three periods, to cover a historic period (P1) from 1899 to 1940, a mid-historic (P2) from 1941 to 1980, and a contemporary period (P3) from 1980 to 2019. We calculated the relative change between the first and second period (relative abundance P2/relative abundance P1), between the second and the third, and averaged it over the two periods. Following Austin & Dunlap (2019), values below 100% indicate a decrease in relative abundance of the species, values above 100% indicate an increase while values of 100% indicate no change in relative abundance. We also calculated species mean observation date and intraspecific variation in observation date. We used the coefficient of variation of observation date (CV = standard deviation/mean) to estimate intraspecific variation in observation date. We calculated the correlation between species activity period (or intraspecific variation in activity period) and species trend in relative abundance by using linear models, with species trend in relative abundance as the response variable and mean activity period (or CV of activity period) as the explanatory variable. Models' residuals were visually checked and showed no sign of heteroscedasticity or non-normality.

Results
For the national scale, citizen-collected data, we found that the flight period of spring bumblebee queens has shifted significantly, advancing on average by 5 days during the past 20 years (Table S3). We found that bumblebee queens were generally observed later in northern, compared to southern parts of Sweden, with a difference of approximately 10 days between the two most extreme latitudes (model 1 in Table 2, Fig. 2a). A significant interaction between latitude and year indicated that, contrary to our expectation, the flight periods of queens at southern latitudes has advanced more compared to northern latitudes, by ~ 8 versus ~ 5 days respectively (Table 2, Fig. 2a, Table S3). At the species level, B. terrestris and B. lucorum were the two species with the earliest observation dates, while B. ruderarius, B. subterraneus, and B. sylvarum were the species with the latest ones (Fig. 2b).
We tested if the effects of time (year) and temperature were interchangeable, that is, if the flight period of bumblebee spring queens shifted with mean annual temperature in the same way that they shifted over time. Results show that when time (year) was substituted for mean spring temperature in the general model (model 1 in Table 2), both the interaction between latitude and temperature, and between temperature and species had significant effects on the flight period of spring queens (model 2 in Table 2, Fig. 3a). Again, the most pronounced advancement of flight period in response to increasing spring temperatures occurred at southern latitudes, compared to northern latitudes (Fig.  S6). The significant interaction between temperature and species suggests differences in sensitivity of species to changes in temperature. Both B. terrestris, B. lucorum and B. pratorum responded to increasing temperatures by advanced queen flight periods, while the other species showed a less marked response (Fig. 3a). The response is thus similar to the one suggested by the interaction between species and year (Fig. 3b, models 1 and 2 in Table 2).
In the analyses including bumblebee phenological trait (early or late queen emergence), we found significant interactive effects between trait and year, and between trait and latitude (model 3 in Table 2) on flight period date. We found that early-emerging species advanced towards even earlier dates over the course of 20 years, while late-emerging species did not shift (Table 2, Fig. 2c, Table S3). For latitude, the significant interaction with emergence date showed that, although both early-and late-emerging species were observed

observation of bumblebee queens in Sweden
The response variable is calendar day. Results of the analysis of variance anova type II (Satterthwaite's method) are reported. Sum squares (Sum Sq), degrees of freedom (Df), F-value and p-value are shown for each model predictor. Conditional R 2 (R 2 (cond)) and marginal R 2 (R 2 (marginal)) are shown per each model. Significant p-values are shown in bold. Non-significant interactions were removed (results not shown) and models rerun to obtain the results. See Appendix for the Tukey tests when applicable. Number of observations: 3012, grouped by grid_id, 685. Number of observations per species: B. hortorum N = 135; B. lapidarius N = 477; B. lucorum N = 361; B. pascuorum N = 452; B. pratorum; N = 466; B. ruderarius N = 89; B. subterraneus N = 44; B. sylvarum N = 107; B. terrestris Fig. 2d, Table S3). For the regional scale, Scania museum dataset, we found that the collection date of spring bumblebee queens was significantly affected by the interaction between landscape type and year (model 4 in Table 3). The time-gap between queens collected in simple versus complex landscapes increased on average by 30 days during the past 20 years, while this gap was almost non-existing in early 1900s (Table 3 Table S4 for model fit values).
When we ran the same model substituting time (year) for temperature, the interaction effect (between landscape type and temperature) was no longer significant (model 5 in Table 3), which is contrary to the results for the citizen-collected data. Instead, we found that increasing temperatures were related to an advancement of the flight period in both landscape types and for all species (Fig. S7). The significant effect of the species factor indicates that species in the museum dataset had different spring queen flight periods (Fig. 4b). In the trait-based analysis, we found that the interaction between species' phenological trait (early or late) and landscape type, had significant effect on queen flight period (model 6 in Table 3, Fig. 4c). Early emerging species were collected significantly earlier (by approximately 13 days) in simple than in complex landscapes (Fig. 4c, Table S4, Table S5), while the collection date of late-emerging species did not differ between complex and simple landscapes (Fig. 4c).
The analysis on the extent to which species decline can be related to species emergence date (Fig. 5a), showed a strong significant correlation between long-term trends in relative abundance and mean calendar day of observation (F 1,8 = 89.2, p < 0.001, r 2 = 0.92; Fig. 5b), indicating that the relative abundances of early species have increased in time, while late emerging ones have decreased.
Moreover, additional analyses showed a significant association between the longterm trends in species relative abundance and the variation (coefficient of variation, CV) in calendar day of observation (p = 0.02, r 2 = 0.47, Fig. S8a), such that species with low variation in observation date were more susceptible to decline than species with higher variation in observation date. This was explained by the negative relation the mean observation date and the CV of observation date across species (Fig. S8b), where early observed species tended to have more varied observation dates compared to late species (p = 0.005, r 2 = 0.65). Table 3 Linear mixed model results describing the effects of species or traits, land-use, and time or temperature on the collection of spring bumblebee queens in Scania The response variable is calendar day. Results of the analysis of variance anova type II (Satterthwaite's method) are reported. Sum squares (Sum Sq), degrees of freedom (Df), F-value and p-value are shown for each model predictor. Conditional R 2 (R 2 (cond)) and marginal R 2 (R 2 (marginal)) are shown per each model.   . 4 Predicted values (marginal effects) for each specific model terms interactions for bumblebee spring queen observations for Scania and the museum collection from 1899 to 2016 for the 10 species included in the study (shaded areas and error bars indicate the 95% CI). a Bumblebee observations across time according to landscape type. b Average first flight day (± 95% CI) in the ten studied species. c Differences between early and late emerging species according to landscape type. *Post-hoc test Tukey test, P = < .0001. See Table 3 for model results

Discussion
In this study, we assessed the phenological shifts in the flight period of the spring-emerging queen bumblebees in Sweden, using museum specimens covering the past 100 years in southernmost Sweden, and citizen observations at the scale of the whole country covering the last 20 years. We found that, in both datasets, the flight period of spring bumblebees has shifted markedly over time and in relation to temperature. In addition, we found that bumblebee queens in southernmost Sweden (Scania) show long-term shifts in phenology related to land-use changes. There are challenges in teasing apart the effects of climate and land-use changes, as they have occurred simultaneously. While some studies have tried to statistically control for both, using historical maps and temperature records (Gérard et al. 2021;Herbertsson et al. 2021), it is still difficult to assess their relative effects. Here we show that both landuse and climate change have caused bumblebee phenological shifts in an area that went through drastic landscape simplification during the twentieth century.

Potential effects of climate change on bumblebee phenology
We show that the flight period across the country has advanced markedly during the past 20 years for several, but not all, bumblebee species investigated. Importantly, these shifts are trait-dependent, such that the group of species previously classified as early emerging has advanced its flight period while the late-emerging group has not -or rather, we do not detect a trend in late-species because some have advanced and some other have delayed their flight period. Queen emergence time in spring is likely triggered by temperature thresholds in the air or soil (Alford 1969;Goodwin 1995;Kudo and Cooper 2019) and therefore increasing temperatures due to climate change are expected to trigger the advancement of emergence in bumblebee queens over time. Differences in duration of spring snow cover across the country may further affect soil temperature and hence queen Relative abundance for each bumblebee species averaged for each sampling periods a, and relationship between the susceptibility to decline and mean observation dates of the ten bumblebee species included in the study b emergence (Kudo and Cooper 2019). As expected, bumblebee queens at northern latitudes in Sweden, and thus exposed to a colder climate and longer periods of snow cover, have a later flight period than bumblebees at southern latitudes.
Temperatures in Sweden have increased during the past decades, with spring temperatures increasing by approximately 1.5 degrees Celsius over the past 100 years (SMHI 2020). In response, an advancement of phenological stages of plants in spring have been observed (Rosenzweig et al. 2007), including onset of flowering and leafing (Menzel et al. 2006). Thus, temperature shifts affect both the phenology of bumblebees and the flowering plants they forage on (Ogilvie et al. 2017;Kudo and Cooper 2019). Due to a higher increase in spring relative to summer temperatures, plant species that flower early have a stronger tendency for an advanced first flowering date compared to species that flower later in the season (Fitter et al. 1995;Lesica & Kittelson 2010). Temperature can therefore also affect bumblebees indirectly by affecting the flowering of plants they forage on. Both direct and indirect temperature effects, separately and combined, may explain why early-emerging bumblebee species show a stronger advancement of flight period over time, compared to late-emerging ones. In addition, in the citizen-collected data we found that, for early emerging species, the advancement of flight period over time is similar to the advancement in response to increasing spring temperatures (Fig. 3). This indicates that during the past 20 years, climate warming is one of the main driving forces behind the shift, and that these species (B. terrestris, B. lucorum, and B. pratorum) may be better able to track optimal temperatures compared to late-emerging ones.

Effects of land-use change on bumblebee phenology
Another possible explanation for the advanced flight period of bumblebee queens is landuse change, which we explored in the museum dataset (focused only on Scania). Scania is a region largely dominated by arable agriculture (Persson et al. 2010). Semi-natural habitats such as permanent grasslands, meadows, and other non-crop habitats provide nesting and floral resources and are therefore essential habitats for bumblebees in agricultural landscapes (Lye et al. 2009;Rundlöf et al. 2014;Svensson et al., 2000). Over the past 60 years, farmland habitats in Sweden, including Scania, have changed considerably. The amount of semi-natural grasslands has decreased, mainly replaced by arable fields or productive forests (Ihse 1995;Auffret et al. 2018). This has led to drastic changes in flora over time, with dramatic declines of arable weeds and grassland species (Tyler et al. 2020). As expected, we found that landscape simplification influenced the flight period of bumblebee queens in Scania. Bumblebee queens sampled in simple landscapes flew on average 14 days earlier than those sampled in complex landscapes. Interestingly, the effect of landscape type on the advanced bumblebee flight period was detected only for early emerging species, and only in simple landscapes, while late-emerging species did not differ in flight period between landscape types.
These results suggest that intensive agricultural management poses a particular pressure on early emerging bumblebee species to adapt their phenology. A possible reason is that floral resources are more restricted in time in simple landscapes, such that midto late season resources may be lacking (Persson and Smith 2013;Timberlake et al. 2019), potentially leading to a higher evolutionary pressure on bumblebees to start colonies early and manage reproduction before flower resources decline. The documented lower diversity of flowering plant species in simple landscapes (Baessler and Klotz 2006;Clough et al. 2014) could also explain the differences in flight period observed 1 3 between simple and complex landscapes. For example, farm specialization indirectly affects the timing of spring floral resources by removing habitats with more temporally dispersed flowering (e.g. extensive leys, semi-natural grasslands, flowering trees) and replacing them with more homogenous and narrow flowering onsets (e.g. oilseed rape, leys with multiple harvests per season). The pulse of resources provided by oilseed rape fields benefits early bumblebee colony growth (Westphal et al. 2009), while positive effects on reproduction may require complementary late-season resources, e.g. in the form of semi-natural habitats or red clover fields (Rundlöf et al. 2014). Because of a more diverse flowering plant community, complex landscapes may likely pose less of a constraint on bumblebees to track resources, since there is a lower risk of resource discontinuity.
In addition, although not specifically tested here, competition for limited nesting sites (as suggested by e.g. McFrederick and LeBuhn 2006) may further induce an advanced flight period in simple landscapes with few non-crop habitats suitable as nesting sites. The advanced emergence of early species in simple landscapes could thus potentially be explained by occupying free nest sites early in the season to avoid competition with lateemerging species (Higginson 2017). This could in turn partly also explain the generally lower abundance of late-emerging species in simple landscapes (Persson et al 2015).
Another possible explanation for the advancement of flight period in simple landscapes may be that soil surfaces are more exposed to atmospheric conditions compared to complex ones, because of fewer trees and less complex vegetation, so that changes in temperatures are more likely to influence emergence date (Pawlikowski et al. 2020). This could similarly affect the onset of flowering in plants.
The interaction between time and emergence trait was not significant in the model of the museum dataset. However, the robust interactions both between landscape and trait, and landscape and time, could suggest that early emerging species are driving the detected trend of advanced flight period, while later species move towards the opposite direction and emerge later. Four of the five early emerging bumblebee species included in this study (B. lapidarius, B. lucorum, B. pratorum, B. terrestris) have increasing or stable populations, according to studies monitoring bumblebee abundance and community composition in red clover fields in southern Scandinavia from 1950 until the early 2010s (Dupont et al. 2011;Bommarco et al. 2011). In contrast, the other early-emerging species, B. pascuorum, has been shown to decline in clover fields according to Bommarco et al., (2011) but increasing according to Dupont et al. (2011). In fact, although considered an early species in this study, B. pascuorum's flight period is intermediate between early and late species (Fig. 2b). In comparison, the late-emerging bumblebee species (B. sylvarum, B. hortorum, B. ruderarius, B. distinguendus) are all classified as declining species (Dupont et al. 2011;Bommarco et al. 2011). According to our assessment of population trends of the ten bumblebee species, late-active species are more susceptible to be declining while earlier-active are more likely to have stable (B. lucorum) or increasing populations (B. lapidarius, B. terrestris, B. pratorum, Fig. 5, Fig. S8).
The fact that early emerging species are active earlier in simple landscapes, and that bumblebees in simple landscapes have an advanced queen flight period dates over the past 110 years (Fig. 3), could be associated with a plastic response to both climate and landuse change. Such plasticity could allow these species to cope better with environmental pressures and global change (McCarty 2001) by shifting their colony cycle to track both resources and optimal temperatures. Correspondingly, the absence of a shift in flight period for late-emerging species could be the result of an inability to shift their colony cycles and track optimal conditions in increasingly warm conditions (McCarty 2001), which could explain their geographical range contraction documented at southerly latitudes (see Nooten and Rehan 2020).

The use of museum collections and citizen science data in detecting long-term trends
Museum specimens and long-term records of citizen collected data are key resources to document and evaluate long-term pollinator declines and their ecological consequences (Bartomeus et al. 2013;Scheper et al. 2014;MacPhail et al. 2019;Zattara and Aizen 2021). We acknowledge the uncertainties involved with both opportunistic and targeted sampling from our museum data, and non-systematic citizen science observations, as these may not consistently reflect the occurrence or abundance of specific species or species groups (Kremen et al. 2011). However, although both our databases covered different spatial and temporal ranges in Sweden, and therefore are not totally comparable, they showed consistent trends.
Datasets based on opportunistic reporting can indeed provide a complementary source of information compared to systematic sampling, especially when deployed at large spatial scales (Henckel et al. 2020;Zattara and Aizen 2021). The most abundant species are usually the easiest to identify and report and they will therefore also drive community trends. In contrast, there are many uncertainties surrounding rare species trends, since these are less abundant and perhaps not detected by observers. This is one of the downsides of citizen collected data and may affect the capacity to capture robust trends for less common species. Nevertheless, we believe that our dataset contains a good sample size to represent all species, traits, and landscape types included (Table S1).
It may be worth mentioning that the span of queen flight periods in the two datasets do not match: the flight period in the national-citizen collected data ranges from 110 to 130 (day of the year) (Fig. 2b) while the regional museum dataset ranges between 155 and 180 (day of the year) (Fig. 4b). This is likely due to the fact that the former dataset is more recent, while museum records span periods where species were on average observed later in the season.

Implications for bumblebee conservation
Climate and land-use change will continue to threaten bumblebees during this century (Soroye et al. 2020). However, the sensitivity of bumblebee species to these global threats depends on a combination of ecological and morphological traits (Bartomeus et al. 2011;Persson et al. 2015). Habitat loss leads to shifts in wild bee species composition, generally benefiting highly mobile and generalist species in agricultural landscapes (Bommarco et al. 2011). As an example, B. distinguendus (late-emerging, long-tongued, and specialized in deep corolla flowers such as red clover), which is red listed in Europe (Nieto et al. 2014), is extinct in Scania since the late 1960s (Cederberg et al. 2021, Fig. S1), mainly due to an extreme decrease in red clover cultivation in southern Sweden (Bommarco et al. 2011, Dupont et al. 2011. Conservation measures are needed in agricultural landscapes to reverse such declines, e.g. by promoting more complex landscapes to provide nesting sites and ensure a continuous provision of floral resources to counteract seasonal resource gaps (Schellhorn et al. 2015;Hass et al. 2018;Timberlake et al. 2019). This could entail conservation of semi-natural grasslands, well designed flower-rich set-asides (flower strips, hedgerows), sowing of clover-rich leys, and the promotion of late mowing (Persson and 1 3 Smith 2013; Klatt et al. 2020;Rundlöf et al. 2014;Piqueray et al. 2019). Our study provides evidence that measures to increase landscape complexity are important conservation actions, irrespective of the effects of climate change on bumblebees.
Phenology studies are key to monitor species trends and provide effective conservation measures in the light of climate warming, and while museum collections provide past records, we here show that engaging the public to track species occurrences, and invest in pollinator monitoring schemes, can indeed provide high quality complementary information (Domroese and Johnson 2017;Breeze et al. 2021).

Conclusion
Several bumblebee species in Sweden are showing a marked advancement in their flight period over the past decades and century. Our results suggest that the temperature increase induced by climate change has a clear and strong effect on the advancement of queen flight period. However, the land-use changes that simplified agricultural landscapes have gone through during the 20th century have also driven species to emerge earlier, likely because of a combination of resource gaps through the summer season (e.g. Persson & Smith 2013;Hass et al 2018) and microclimatic conditions that fluctuate more than in complex landscapes (Schmidt et al. 2019).
Landscape simplification is known to select for several bumblebee traits, for example, increased queen size (Gérard et al. 2020a, b), and species with larger colonies are more common in simple landscapes (Persson et al. 2015). In addition, we here show that spring queens are emerging earlier in simple landscapes. As previous studies focused on speciesspecific traits, they could only document variation in community-traits values resulting from a change in species composition (filtering effect). In addition to these communitylevel processes, our study suggests that the intraspecific variation in traits such as flight period may also drive communities' response to climate and landscape changes. We show that some species, especially the ones with a late activity period, might not be able to track shifts in temperatures and/or resources in simple landscapes, possibly leading to a higher likelihood of local extinctions. Our study therefore calls for the implementation of management interventions to counteract the effects of climate warming on pollinators and their flower resources, and thus ensure that shifts in the activity periods of some dominant species do not harm the integrity of rare ones (Herbertsson et al. 2021).