Impact of landscape fragmentation and climate change on body size variation of bumblebees during the last century

Body size is a key parameter of organism fitness. While the impact of climate change on body size has received increasing attention, the long-term consequences of landscape fragmentation are still poorly known. These two major global threats may potentially induce opposite trends: the decrease of body size in warmer environments (e.g. individuals developing faster) or the selection of larger individuals in fragmented habitats (e.g. large individuals more capable of reaching distant patches). We assessed the relationship between temperature and landscape fragmentation with mean body size during the last century, within four European regions (Austria, Belgium, England and above the Arctic circle in Scandinavia) and among queens of five bumblebee species. At the regional scale, we first analysed the variation over time of body size and the two hypothesised drivers, temperature and landscape fragmentation. Then, at the local landscape scale, we tested whether body size varied according to these drivers irrespective of the region. At the regional level, we observed a statistically clear increase of queen body size corresponding to an increase of landscape fragmentation (i.e. in Belgium and England). There was no increase of size when fragmentation did not increase (i.e. in Austria and above the Arctic Circle). Temperature also increased through time in all regions. At the local landscape scale, we found that all species were impacted by changes in both climate and landscape fragmentation but show different trends. The body size of the two largest species significantly increased at landscape level with higher fragmentation while body size of the two smallest species decreased with higher fragmentation. We highlight that, in a context of global changes, landscape fragmentation can also be a major driver of body size clines. Depending on the dispersal abilities of species, larger species could be positively selected for and overcome landscape fragmentation.


Introduction
Body size is a fundamental trait to understand fitness, from metabolic processes (Kleiber 1932, Brown et al. 2004 to intra-specific interactions and ecosystem functioning (Woodward et al. 2005). Shifts of body size can thus alter ecosystem stability, affecting both primary production and biodiversity (Woodward et al. 2005, Yvon-Durocher andAllen 2012). Variation in body size is caused by a range of factors; the most studied being temperature and resource availability (e.g. through habitat fragmentation) in a context of latitudinal or altitudinal clines (Atkinson andSibly 1997, Meiri et al. 2007). In this context, Bergmann's rule is probably the most widely known relationship (Bergmann 1847). This rule states that higher ambient temperatures are expected to lead to smaller body sizes in homeothermic vertebrates; later it was postulated that increased heat loss in smaller endotherms and higher metabolic rates in smaller ectotherms could explain these trends (Blackburn et al. 2008). While there is a disagreement if Bergmann's intention was to state this rule at the intra-or inter-specific level, it is now commonly used 'sensu lato', without concern of the level and the mechanism (Shelomi 2012). Lower resource availability can lead to smaller body size too, especially in higher latitudes where the resources can be spatially and temporally scarce (Borthagaray et al. 2012).
Variation of body size across long-term temporal clines has received less attention than latitudinal or altitudinal clines. When it comes to global changes, climate change is by far the most extensively studied factor that has been described as affecting body size (Sheridan and Bickford 2011). By analogy with latitudinal clines, a decrease of body size consequent to climate change was predicted among a wide range of taxa and it has even been proposed as a universal response to climate change (Millien et al. 2006, Gardner et al. 2011. By contrast, the impacts of habitat loss and landscape fragmentation on body size have received far less attention even though they are among the major drivers of biodiversity loss (Hanski 1998, Dirzo andRaven 2003) and can also alter ecosystem processes (Valladares et al. 2006) and filter species traits (Gamez-Virues et al. 2015). Larger individuals seem to be selected because of their abilities to cross unsuitable habitat to reach distant habitat patches (Tscharntke and Brandl 2004, Stevens et al. 2014, Keinath et al. 2017, Merckx et al. 2018 and are also often more resistant to starvation (Peters 1983, Gergs andJager 2014). Recent theoretical frameworks encapsulate the positive relationship between body size and landscape fragmentation (Hillaert et al. 2018a, b), but convincing empirical evidence is lacking, particularly in invertebrates. Moreover, all these drivers of body size clines can also counteract each other, leading to very slight or nonsignificant body size clines. For example, while landscape fragmentation during the last century can select for larger individuals, lower resource availability or climate change can lead to smaller ones so the trade-off between the energy allocated for movement and the energy allocated for growth are dependent on the degree of habitat discontinuities, resources access, climatic conditions and the dispersal ability of the species (Peters 1983, Bonte et al. 2012. Bees are an ideal model to study the long-term impacts of global changes including landscape fragmentation and climate change on body size. Indeed, many entomologists have collected them throughout the world during the last century and several species have long time-series in museum collections. Moreover, bees are experiencing declines in species richness and abundance in certain areas (Goulson et al. 2015); both climate change and habitat loss/isolation are involved as major causes and will potentially provide selection pressures leading to phenotypic changes. While it is predicted that climate change will severely impact some bee species during the next decades (Rasmont et al. 2015, Ollerton 2017, loss and fragmentation of habitat are already identified among the main drivers of decline (Potts et al. 2010, Goulson et al. 2015, but see Cane 2001, Smith and Mayfield 2018 for contrasting results) and are expected to continue to impact bee abundance and distributions in the future alongside climate change (Marshall et al. 2018). Bee body size can decrease as a result of the increase in temperature (Scriven et al. 2016, Gérard et al. 2018a, b, Nooten and Rehan 2019 even if the pattern is not consistent (Dellicour et al. 2017, Gérard et al. 2018b). However, bee body size can also increase in fragmented landscapes (Warzecha et al. 2016, Gérard et al. 2019). Among bees, it is commonly observed that larger species can fly longer distances even if this is not an absolute rule as honeybees can fly much further than bumblebees although they are smaller (Osborne et al. 2008, Pahl et al. 2011. Nonetheless, many studies suggest that larger individuals can be favoured over smaller ones in fragmented landscapes because of the disconnection between nesting and foraging locations, (Gathmann and Tscharntke 2002, Greenleaf et al. 2007, Wright et al. 2015 even if alternative parameters like nest density and diet breadth also influence foraging range (Knight et al. 2005).
In the context of an expected continued increase of landscape fragmentation and climate change, we urgently need to understand how key species like wild bees phenotypically respond to these global changes. Even if a wide array of factors may affect body size (e.g. pesticide), we focused on these two drivers because, based on the literature, they are expected to have been the main factors affecting bumblebee body size during the last century, and because of the sufficient amount of data available across this time. Based on museum collections from the last century, recent bumblebee specimens and historical environmental data, we first assessed the trend of queen body size within five European bumblebee species (Bombus balteatus, B. hortorum, B. lapidarius, B. pascuorum, B. pratorum) as well as landscape fragmentation and temperature across four regions (Belgium, England, Austria and Scandinavia above the Arctic circle) over time. We expected that both landscape fragmentation and temperature will vary among the four regions. Secondly, we tested the effects of climate change and landscape fragmentation at the local landscape scale, i.e. the sites where specimens were collected. We considered the surrounding conditions of the landscapes where specimens were collected, and analysed the trends in body size taking into account the inherent variability between regions. We hypothesize that 1) if climate change is the main driver of queen body size variation, we will observe a decrease of body size if landscape fragmentation does not counteract the effects of climate change. However, if 2) landscape fragmentation is the main driver of queen body size, we will observe an increase of body size in regions with increasing landscape fragmentation. Moreover, we expect that if body size only increases in regions with increasing fragmentation, we would observe an increase in size along a fragmentation gradient in the analysis at local scale, irrespective of the region.

Sampling and body size measurement
We focused our study on the bee genus Bombus (i.e. bumblebees). Bumblebees represent an interesting model to explore changes in body size: 1) they are heterotherms (i.e. their body temperature vary depending on self-regulating mechanisms and surrounding environment temperature), 2) they can be found in both arctic and temperate regions and 3) they have been collected by many naturalists over the last century (Heinrich 1993). We only considered regions (i.e. central, western and northern Europe) and bumblebee species that were represented in the museum collections over long time series. Specimens of Bombus balteatus were only available in northern Scandinavia, B. hortorum, B. lapidarius and B. pratorum specimens were collected in Austria, Belgium and England and finally B. pascuorum specimens were available from the four regions (Austria, Belgium, England and Scandinavia).
We evaluated bumblebee body size using inter-tegular distance (i.e. ITD, distance in millimetres between the two insertion points of the wing) as a proxy, this distance is measured on the thorax which contains flight muscles and is also strongly correlated with other size parameters like dry body mass and wing size (Cane 1987, Greenleaf et al. 2007, Gérard et al. 2018a. The same researcher measured the ITD of all specimens using a Facom 1500 mm digital calliper (France, Morangis). We selected a total of 2457 bumblebee queens collected between 1900 and 2017: 80 queens of Bombus balteatus, 509 queens of B. hortorum, 632 queens of B. lapidarius, 728 queens of B. pascuorum and 508 queens of B. pratorum (Supplementary material Appendix 1 Table  A1, Fig. A1 for the details by region). For each year, we sampled between three and ten specimens per species and per region and calculated the mean ITD (mean number of specimens per year = 5.432, SD = 2.499). We selected only spring queens (i.e. queens collected from February to July) to try and select only first-generation queens, because they are strongly linked to the fitness and the success of the colony during the early stages of colony establishment (i.e. nest searching and early collection of food resources; Pyke 1978). All sampled specimens are conserved in museum collections at Gembloux Agro-bio Tech (Belgium), Royal Belgian Inst. of Natural Sciences (IRSNB, Belgium), Univ. of Mons (Belgium), Univ. Museum of Bergen (Norway), Natural History Museum of London (UK) and Biologiezentrum Linz (Austria).

Temperature data
Historical interpolations of mean monthly minimum and maximum temperatures for the period 1901-2016 were downloaded from the CHELSAcruts database at a resolution of 30 arc-sec, approximately 1 km (Karger et al. 2017). Due to the low density of weather stations before 1950, the climate data from the early part of the 20th century is not as accurate as more modern data. However, the CHELSAcruts dataset represents the best source to provide comparable climate trends across the spatio-temporal resolution of our study.
For each year we calculated the mean annual minimum and maximum temperatures per month (minT, maxT). For the analysis at the regional scale we extracted the mean value per decade of minT and maxT for each of the four regions. For the analysis at the landscape scale we created a 5 km buffer around each unique specimen observation and extracted the mean minT and maxT values for the year the specimen was collected and averaged these values over the entire buffer. For the specimens collected in 1900 and 2017, outside the range of the CHELSAcruts database, we used 1901 and 2016 as proxies. The two variables minT and maxT were positively correlated (Pearson's correlation coefficient [r] = 0.7), therefore we selected minT to use in the models as it was less correlated to the chosen landscape fragmentation variable (r = 0.04 and 0.2 respectively; Supplementary material Appendix 1 Fig. A3).

Landscape fragmentation data
Bumblebee species included in these analyses are overall openhabitat generalists (Rasmont et al. 2015). Unfortunately, detailed data on the habitat requirements for each of the bumblebee species were not available for all the sampled localities across the last century. Therefore, we considered landscape fragmentation as a substitute for habitat fragmentation. In order to examine comparable landscape fragmentation in each of the four regions across the last century, we required long-term historical land use/land cover (LULC) maps at a European scale for the 20th century. Spatially explicit, comparable maps of LULC were not available at the spatial and temporal scales of the study. Therefore, we established a proxy for landscape fragmentation over time using historical reconstruction maps for Europe, which are available per decade from 1900 until 2010 (Fuchs et al. 2015a, b). These datasets combine historical maps of LULC alongside a modelled reconstruction using Historic Land Dynamics Assessment or HILDA as the modelling approach (Fuchs et al. 2015a, b). The maps represent a thematic resolution of six LULC classes at a 1 × 1 km resolution, and include forests, grasslands, cultivated land, human settlements, water and other. The HILDA maps do not include Norway, therefore only the LULC from Finland and Sweden, above the arctic circle, were included in the regional scale analysis and the specimens from above the Arctic Circle were excluded for the analysis at the local landscape scale.
No single metric can adequately measure landscape fragmentation (Wang et al. 2014). We included metrics of dispersion (spatial distribution of individual land use classes) and inter-dispersion (the interactions between different land use classes) indices, as well as measures of edge and patch density. We calculated 11 separate landscape fragmentation metrics (Supplementary material Appendix 1 Table A2) using the 'landscape metrics' R package (ver. 1.1; Hesselbarth et al. 2019). For the regional scale analysis, we extracted the landscape fragmentation metrics as a single value for the whole region per decade. For the landscape analysis we extracted the landscape fragmentation metrics from a 5 km buffer around the coordinates of each specimen for the decade when the specimen was collected. All indices were strongly correlated (Supplementary material Appendix 1 Fig. A3). Therefore, we chose to use a single metric, the contagion index, because it incorporates measures of both dispersion and interdispersion of different land use types, giving a more complete picture of how the landscape is arranged (McGarigal and Marks 1995). Furthermore, the contagion index was strongly correlated to all other landscape metrics (r < −0.7 or > 0.7) (Supplementary material Appendix 1 Fig. A3). The contagion index is measured as a proportion/percentage and incorporates the proportions of the landscape occupied by each land use class, the number of shared edges between each land use class and the total number of classes present in the landscape (for full details on the calculation, see McGarigal 2015). To aid the interpretation of the contagion index, i.e. show increasing landscape fragmentation, we used the inverse of the contagion index. An inverse contagion index of 0 means there is only a single land use type in the surrounding landscape (lowest fragmentation) and when the values is 100 the land use types are maximally disaggregated (highest fragmentation; McGarigal 2015).

Body size and environmental change trends at regional scale
To test whether there were similarities in the temporal trends of body size, temperature and landscape fragmentation for each of the four regions, we used a multiple linear regression approach. To have reliable trends in body size and to compare with the temporal resolution of the landscape data, we took decadal means (from the eleven last decades) of each of the three dependent variables: body size, temperature, and landscape fragmentation. We calculated multiple linear regression models to predict each of the dependent variables based on decade, region and their interaction as independent variables. For the mean decadal body size, we also include species as an extra independent variable, as the different species vary in their average sizes. For landscape fragmentation the model did not meet the assumptions of linear regression. We used the betareg R package (v3.1.3; Cribari-Neto and Zeileis 2010) to fit a beta-distribution to deal with non-normality of the residuals and heteroscedasticity, and because landscape fragmentation is measured as a proportion. We also conducted post-hoc pairwise comparisons to observe the difference in trends between regions using the emmeans R package (v1.4.2; Lenth 2019).

Body size trends at local landscape scale
To understand the trends and drivers of queen body size at local landscape scale we analysed the individual specimen data with a linear mixed effects model with a Gaussian distribution (LMM) using the lme4 R package (v1.1-21, Bates et al. 2015). We extracted the mean annual minimum temperature and the inverse contagion index within a 5 km buffer surrounding the coordinate value of each specimen. As specimens above the Arctic Circle were excluded from this analysis, we fit the model with 2292 individual specimen records. Due to aforementioned collinearity between predictors, we used minT and the contagion index (landscape fragmentation) alongside year of collection and species as the fixed effects (predictors) for the model. We fitted a full LMM with body size as the response variable and included all two-way interactions between the chosen predictor variables, with a Gaussian distribution. Multiple specimens were collected with same the coordinates from locations we define as sites, and these sites were found throughout our three regions (Austria, Belgium and England). Therefore, we included site nested within region as a random effect and not as a fixed effect, to account for spatial dependency between specimens from the same region, to account for spatial autocorrelation of body size from nearby regions and because we were interested in the overall trends of body size and not the inherent variation between sites and regions. We selected the best model as the model with the lowest AIC (Burnham and Anderson 2004) after testing all possible model combinations. Model visualisations were then produced using the sjPlot R package (v2.7.2; Lüdecke 2019). All analyses were conducted in RStatistics (R Core Team).

Queen body size and environmental trends at regional scale
Preliminary results using single linear models for each species in each region show that every species in Belgium and England increased significantly in mean body size over time Multiple linear regressions were then fit to predict temporal trends in body size, temperature and landscape fragmentation in each region. We observed statistically clear linear trends for all three dependent variables over time. Both body size and landscape fragmentation were predicted to have different trends over time depending on the region. In Belgium and England queen bumblebees demonstrated statistically clear increases in body size per decade (0.03 and 0.02 mm per decade, Table 1), these trends were significantly (p = 0.03 and p < 0.001 respectively) different from the trends observed in Austria where no clear increases were demonstrated (−0.002 mm per decade; Table 1, Fig. 1A). The models predicted comparable trends for increases in landscape fragmentation and increases in body size of queen bumblebees in Belgium (Table 1 see 'Scaled trend'). In Belgium and England, the model predicted statistically clear trends of increased landscape fragmentation per decade (0.78% and 0.94% per decade), significantly different from Austria (p < 0.001 and p < 0.001 respectively), which showed no clear trend (−0.01% per decade) and above the Arctic Circle (p < 0.001 and p < 0.001 respectively), where landscape fragmentation was predicted to have significantly decreased per decade (−0.43% per decade, Fig. 1B). Minimum temperature on the other hand was predicted to have no clear difference between the trends over time in the four regions. Each region shows an approximate 0.1°C increase per decade (Table 1, Fig. 1C).

Queen body size trend at local landscape scale
The best model (next best model: ΔAIC 1.4, for other models with ΔAIC < 2, Supplementary material Appendix 1 Table  A4-A6) to explain body size patterns over time included year, minT, landscape fragmentation and species, as well as the interactions of species with minT and species with landscape fragmentation. The interaction between species and both minT and landscape fragmentation illustrated that the four species included in the analysis do not show a consistent relationship with the environmental drivers. Conversely, we did not observe a clear statistical interaction between date and species suggesting that all species do show a similar trend in increased size over time. Indeed, we observed an average increase of body size over time for each of the four species (0.002 mm yr −1 ; Fig. 2A Table A3). Finally, there is also a statistical clear and strong effect of landscape fragmentation. Overall, as landscape fragmentation increases at sites the model predicts larger specimens. Again, this is highly dependent on species. The structure of each model consists of y = decade × region, y being each of the three response variables. All data were averaged for each decade from 1900 until 2010. Significance is measure as p < 0.05 = *, p < 0.01 = **, p < 0.001 = ***. The scaled model trend is based on the response variable being transformed to have a mean of 0 and standard deviation of 1. The scaled trends are to aid in the comparison between models.  Table A2).

Discussion
The results show that queen body sizes have, on average, increased during the last century, whereas we would expect an overall decrease of body size if climate change was the only dominant driver of body size changes. In parallel, landscape fragmentation increased in two countries (Belgium and England), remained stable in Austria and decreased above the Arctic Circle. In this first analysis, no body size decrease was observed over time as we would expect if increasing temperature was the main driver. Interestingly, the results for the drivers of body size trends at the local landscape scale showed a more complex relationship between body size and the hypothesized drivers. Overall, all species showed an increase in size over time when analysed across the whole study area, but showed varying relationships with climate and landscape fragmentation. This indicates that there is not a consistent role of climate and landscape fragmentation on body size, even within a single genus. Indeed, we see that the two drivers have a contrasting relationship with body size which overall leads to the small increases in size over time. While we observed a positive relationship between body size and landscape fragmentation, and a negative relationship with climate for two species (i.e. Bombus hortorum and B. lapidarius), we observed the inverse relationship for the two other species (i.e. B. pascuorum and B. pratorum) with smaller bodies observed in more fragmented landscape and larger specimens in warmer regions. This suggests that the overall trends observed at the regional level can be refined with local scale data.
Most of the literature in recent decades has focused on the relationship between climate change and body size. However, in the analyses focusing on each region separately, bumblebee body sizes did not decrease alongside an increase in temperature, which is a good indicator that temperature per se did not explain the overall body size clines. This absence of a clear systematic impact of temperature could be due to the isolation of the nest from ambient temperature (i.e. underground), which could buffer body size variation against temperature fluctuations to a certain extent (Heinrich 1993). Habitat fragmentation is far less studied as a driver of body size over time. In insects, a general trend towards larger individuals along a positive spatial gradient of fragmentation has been described, for example, in butterflies (Hill et al. 1999), damselflies (Taylor and Merriam 1995), bush crickets (Berggren 2005), and in one genus of solitary bees (Warzecha et al. 2016). Among bumblebees, this pattern could be explained by the selection of larger queens in spring. Larger queens could be positively selected for because of their greater dispersal abilities in more fragmented habitats to find nesting locations and collect the first resources needed for colony initiation (Greenleaf et al. 2007, Wright et al. 2015. Interestingly, all these clades have high dispersal ability; our findings thus corroborate the recent review of Merckx et al. (2018). Indeed, in urbanized habitats, both decreases and increases of body size have been observed, depending on the taxa and their dispersal ability. When a taxon has dispersal abilities that allowed them to compensate for the increase of landscape fragmentation (e.g. larger moths, butterflies, orthopterans), body size increased with increasing urbanisation, while in less mobile taxa (e.g. spiders, weevils, cladocerans), they observed a decrease of body size with increasing landscape fragmentation (Merckx et al. 2018). Bumblebees have relatively large body sizes (approximately 4-8 mm of thorax width for queens, Gérard et al. 2018b) and high dispersal abilities (> 10 km, Greenleaf et al. 2007). For the larger species, the body size trend during the last century is thus probably dependent on a trade-off between the degree of landscape fragmentation and the dispersal capacities of model species. We could expect an increase of body size only in species with high dispersal capacities (i.e. larger bumblebees) while small-ranging species (i.e. smallest solitary bees) may mostly display the opposite trend (Merckx et al. 2018). Interestingly, when we computed the analysis on the raw specimen data (Fig. 2C), the body size of the two largest species (i.e. with higher dispersal abilities) increased with higher landscape fragmentation while the body size of the two smaller species (i.e. with lower dispersal abilities) decreased with higher landscape fragmentation, in congruence with the hypothesis of Merckx et al. (2018). The selection of larger body size in spring queens consequently to higher landscape fragmentation could also be mitigated by climate change. Indeed, warmer wintering temperatures could limit the need to use fat reserves during diapause and thus increase dispersal abilities at the end of hibernation. When considering the effect of temperature on the raw specimen data at the local scale (Fig. 2B), this driver seems indeed to lead to intermediate body size for each of the four species. In parallel with our results, a recent paper assessed the evolution of body size of four bumblebee species during the last century, in USA (Nooten and Rehan 2019). They found the opposite trend: the body size of the four species had declined during the last century. They do not assess the drivers of this variation and assessed this trend on workers. One of their main explanations is that the scarcity of food resources resulted in a decrease of body size over time, as bee body size is primarily linked to the quantity of food resources. As their study has been conducted on workers, the selection pressure on body size could differ from queens, which are in charge of the first stages of colony foundation. Furthermore, their dataset is less homogeneously distributed throughout the century making it difficult to attribute any driver to the observed changes in body size.
Alternatively, the increase in body size could have been caused by agricultural intensification and the use of pesticides, which are correlated with landscape fragmentation (Tscharntke et al. 2005). Larger body size can improve physiological resistance to pesticides (Tahori et al. 1969, Thompson andHunt 1999). Thus, with the increase of pesticide spread during the last decades, larger bumblebees could have again been positively selected for. Unfortunately, statistically and quantitatively testing pesticide use during the last century to assess its relationship with body size modification remains elusive because no standardized dataset is available to measure this variation over time across different European countries. Bumblebee body size is also strongly dependent on diet quality and quantity. It seems that this parameter remained quite stable during the last decades for generalist bumblebees and could not explain body size variation (Kleijn andRaemakers 2008, Roger et al. 2017). These claims need however to be qualified: warmer temperatures can increase the flowering period and thus the resources availability for larvae, which could also lead to an overall increase of body size through time as it has been observed at the local scale for each species. Gardner et al. (2011) pointed out that food availability or quality is indeed one of the main drivers of body size changes and it cannot be excluded. A possible consequence of body size mismatches between plants and pollinators would be the disruption of their mutualistic network as body size strongly constrains this interaction and its morphological mechanistic fit (Bartomeus et al. 2016, Kendall et al. 2019. Body size can impact, for example, the activity period of a bee (i.e. smaller bees are less prone to forage in dim light; Streinzer et al. 2016) but has also a strong influence on pollination success which can decrease with decreasing body size (Armbruster et al. 1989, Stang et al. 2009, Solis-Montero and Vallejo-Marin 2016. A morphological mismatch in this interaction could thus decrease pollination effectiveness and the delivery of ecosystem services (Jauker et al. 2016, Gérard et al. 2020a.
Our study highlights the importance of taking into account several environmental parameters when assessing the modification of a morphological trait. In the future we recommend assessing these shifts of body size along a long time series of species with different levels of dispersal ability to better understand how landscape fragmentation can induce body size modifications and their consequences on ecosystem functioning. We also encourage experimental researches in mesocosms to distinguish the effect of temperature and fragmentation on body size.