Signatures of increasing environmental stress in bumblebee wings over the past century: Insights from museum specimens

Abstract Determining when animal populations have experienced stress in the past is fundamental to understanding how risk factors drive contemporary and future species' responses to environmental change. For insects, quantifying stress and associating it with environmental factors has been challenging due to a paucity of time‐series data and because detectable population‐level responses can show varying lag effects. One solution is to leverage historic entomological specimens to detect morphological proxies of stress experienced at the time stressors emerged, allowing us to more accurately determine population responses. Here we studied specimens of four bumblebee species, an invaluable group of insect pollinators, from five museums collected across Britain over the 20th century. We calculated the degree of fluctuating asymmetry (FA; random deviations from bilateral symmetry) between the right and left forewings as a potential proxy of developmental stress. We: (a) investigated whether baseline FA levels vary between species, and how this compares between the first and second half of the century; (b) determined the extent of FA change over the century in the four bumblebee species, and whether this followed a linear or nonlinear trend; (c) tested which annual climatic conditions correlated with increased FA in bumblebees. Species differed in their baseline FA, with FA being higher in the two species that have recently expanded their ranges in Britain. Overall, FA significantly increased over the century but followed a nonlinear trend, with the increase starting c. 1925. We found relatively warm and wet years were associated with higher FA. Collectively our findings show that FA in bumblebees increased over the 20th century and under weather conditions that will likely increase in frequency with climate change. By plotting FA trends and quantifying the contribution of annual climate conditions on past populations, we provide an important step towards improving our understanding of how environmental factors could impact future populations of wild beneficial insects.


| INTRODUC TI ON
Environmental change over the past century has placed multiple pressures on biodiversity, threatening important functional groups of organisms (Carmona et al., 2021;De Palma et al., 2016;Grab et al., 2019;Lenzen et al., 2012;Tilman et al., 2017;Trisos et al., 2020). Especially worrying have been the reports of declining insect pollinators Nieto et al., 2014;Powney et al., 2019;Zattara & Aizen, 2021), as the vast majority of angiosperms (Ollerton et al., 2011) including >75% of the world's leading food crops (Klein et al., 2007) are to some degree dependent on their pollination service. Ultimately, if pollinators disappeared, half of the c. 350,000 species of flowering plants would lose an estimated ≥80% of their seed production, with a third of flowering plants estimated to fail to produce any (Rodger et al., 2021). Thus, distinguishing which factors are contributing to insect-pollinator losses remains a research priority, and will help inform safeguarding strategies and predict future pollination services. Several global factors are associated with insect-pollinator losses, including climate change, agricultural intensification and spread of invasive species (Ollerton, 2021;Potts et al., 2016). However, we are still limited in our ability to quantify the degree to which such factors have contributed to stress being placed on wild populations (Gill et al., 2016;Ollerton et al., 2014); defined here as a reduction in an individual's energy allocation towards development and reproduction (Beasley et al., 2013). We have typically been constrained to making inferences about past responses of wild populations based on occupancy modelling (Casey et al., 2015;Schleuning et al., 2016), laboratory studies (e.g. Kenna et al., 2021;Martinet et al., 2021;Smith et al., 2020) or snapshot comparisons of populations between past and contemporary landscapes (Carvell et al., 2017;De Palma et al., 2016).
But this often does not account for the fact that the magnitude of these environmental stressors has varied historically in space and time. Part of the challenge of quantifying past stress is that time-series data spanning the past century and beyond are scarce, stifling progress in understanding when historic insect-pollinator populations actually experienced stress. Moreover, without comprehensive temporal data and a baseline understanding of stress experienced by insect-pollinator populations, it is difficult to place any change into an appropriate historical context.
A solution to help fill this data gap is to leverage museum specimens, which represent an underutilised biological data repository, to investigate how historic populations have responded to stressors (Freedman et al., 2020;Heberling, 2020;Oliveira et al., 2016;Owens et al., 2020;Scheper et al., 2014). By combining collection data across multiple museums, we can also overcome spatiotemporal biases in sampling (Wildman et al., 2022).
Mapping when and where stress has occurred is complicated, given that different stressors can impose varying lag effects on population sizes and geographic distributions (Wearn et al., 2012).
This makes associating population declines or range contractions with historic environmental change challenging. An alternative approach is to use a biomarker or morphological signature as a proxy for levels of stress at the time environmental pressures were experienced, thereby enabling a more accurate assessment of the effects of such stressors. One such biomarker and potential proxy of stress is fluctuating asymmetry (FA), which are deviations from the bilateral symmetry of body plans (Palmer, 1994), reflecting developmental stability. Developmental stability refers to the ability of an organism to produce a consistent phenotype under different conditions (Palmer, 1994). As the left and right sides of bilaterally symmetrical traits are determined by the same genes, deviations from perfect symmetry are thought to be due to exogenous factors affecting the developmental process, for example the production of stress-induced hormones (Benderlioglu, 2010). Although random deviations from asymmetry are expected in traits under normal conditions, high FA or increasing levels of FA are thought to reflect an inability to buffer development against suboptimal conditions, such as those induced by environmental change (Freeman et al., 1993;Palmer, 1994;Palmer & Strobeck, 1986). Directional changes in the levels of FA over time could therefore indicate that the underlying drivers of FA (that we interpret here to be 'stressors') are also changing. 4. Species differed in their baseline FA, with FA being higher in the two species that have recently expanded their ranges in Britain. Overall, FA significantly increased over the century but followed a nonlinear trend, with the increase starting c. 1925. We found relatively warm and wet years were associated with higher FA.
5. Collectively our findings show that FA in bumblebees increased over the 20th century and under weather conditions that will likely increase in frequency with climate change. By plotting FA trends and quantifying the contribution of annual climate conditions on past populations, we provide an important step towards improving our understanding of how environmental factors could impact future populations of wild beneficial insects.
For example, bumblebee Bombus terrestris wings have been shown to become less symmetrical under temperature (33°C) or parasitic stress (Gerard et al., 2018), as do honeybee Apis mellifera wings when individuals are exposed to pesticides (Friedli et al., 2020). Although some studies have not detected a link between FA and stress (e.g. Bjorksten et al., 2001;Floate & Fox, 2000), a meta-analysis found stressors explained 36% of variation in FA across insect studies and concluded FA is a sensitive biomarker (Beasley et al., 2013). Indeed, FA has been ob- but see Clarke, 1995); yet, to our knowledge, long-term studies investigating changes in FA are lacking.
Here we present a study using historical collections from a network of British museums to investigate changing FA as a signature of stress experienced by four bumblebee species (Bombus hortorum, B. lapidarius, B. muscorum and B. pascuorum) within Britain throughout the 20th century. Bumblebees are globally important pollinators (Goulson, 2010), and considered sensitive to environmental change (e.g. Settele et al., 2016), making them an important functional group to study. These four species were chosen as their distributions across Britain have changed differently over the past century: B. hortorum and B. lapidarius have maintained a relatively stable geographic distribution since the 1960s; B. muscorum has declined; and B. pascuorum expanded its distribution (Casey et al., 2015). To quantify levels of FA in these historic bumblebee populations, we landmarked photographs of forewings and collected label metadata for >3,300 specimens from five museum collections. After conducting a set of rigorous filtering steps (see Section 2) to consider only high-quality specimens, we investigated: (a) if the amount of FA varies between species, and whether higher or lower FA is associated with distribution change over the century (described above); (b) if FA is higher in the second half of the 20th century given anthropogenic impacts on the environment have increased over the century, and whether change has shown a steady and consistent trend (linear relationship) or a more dynamic trend (nonlinear); and (c) if there is a correlation between FA and the climatic variables, mean annual temperature and precipitation, given both have been implicated in impacting bumblebee biology.

| MATERIAL S AND ME THODS
We first photographed the museum specimens of the four bumblebee species and extracted corresponding metadata (such as collection location) from each label. We summarised wing shape by manually landmarking the left and right forewings, and used the difference in the left and right wing shapes to estimate FA. Next, we took a conservative approach to specimen inclusion by applying several data-filtering steps (detailed below), to exclude specimens that might have artificially high FA (e.g. through having tilted wings as a consequence of the curatorial process). Additionally, to ensure an even distribution of data across the 20th century, we rarified our data (also detailed below). We did not require ethical approval for this study. Lite 5500K fluorescent ring bulbs. From these, we extracted sample metadata from the label images, by transcribing collector name, collection date, location and caste, where present. All specimens had been collected throughout Britain and geolocating of specimens was carried out using Google Maps' Geocoding application programming interface, which returned a formatted address, country of origin and latitude and longitude coordinates (of the location's centroid using Mercator projection) for each specimen, using R code adapted from Lynn (2014). Location precision ranged from country to post code ( Figure S1), with 66.6% of specimens located to city level ('locality') or higher resolution.

| Initial dataset
To associate each specimen with climatic data, we grouped locations into UK Met Office (2020) 'climate regions' (East Anglia; east and northeast England; northwest England and north Wales; southeast and central England; Midlands; south Wales and southwest England; east Scotland; north Scotland; west Scotland).
Climatic data in the UK Met Office's UK and Regional Series are available from 1910 for each climate region as monthly, seasonal and yearly values. We used mean annual temperature and precipitation, as annual precipitation correlated with the precipitation for each season and the number of days with ≥1 mm rain in each season (Pearson's correlation: r ≥ 0.50, p ≤ 0.001), and mean annual temperature correlated highly with the yearly and seasonal maximum, minimum and mean temperatures (r ≥ 0.57, p ≤ 0.001; Figure S2; but we include, in the Supporting Information, additional analyses based on maximum and minimum annual temperatures). Although annual climatic data encompass months in which bumblebees are not developing, they provide an overview of the climatic conditions throughout each year, and they remove the need to make assumptions about the age of each bumblebee when collected, to infer when it had developed. However, we also provide analyses using climatic data from just the spring and summer months in the Supporting Information.
We focused on specimens collected between 1900 and 2000, which were represented by greater availability of data, with a recorded collection location and date, and with wings positioned in a manner in which venation could be assessed from a dorsal image (Bombus hortorum = 656; B. lapidarius = 818; B. muscorum = 511; and B. pascuorum = 1,351 individual specimens; Table S2).
Additionally, we only considered queens collected later in the year. For example, for B. hortorum, we excluded queens collected before Julian day 150, as these likely represent queens that had developed in the year preceding collection and had overwintered (Julian day for B. lapidarius = 164; B. muscorum = 164; B. pascuorum = 159; Methods S1).

| Wing landmarking
To capture forewing shape, we used the software tpsDig (v. 2.31; Rohlf, 2015) to manually plot 13 landmarks, which are homologous across bumblebees and have been used previously to investigate FA (Klingenberg et al., 2001; Figure 1). This resulted in two sets of 13 2D Cartesian coordinates per bumblebee (one set per forewing). Specimens were not re-mounted or re-positioned in any way before taking photographs to avoid damaging these irreplaceable and fragile specimens. Two people performed the landmarking, and to ensure consistent positioning of landmarks, we performed a repeatability analysis, where a subset of images (20 B. hortorum drones) were landmarked every 3 weeks over the landmarking period by both data collectors. For each specimen, there was no significant difference between wing shape measurements taken by the two data collectors over the landmarking period (linear mixedeffects models; F = 0.203, p = 0.654; Table S3; Figure S3). Only images where all 13 landmarks were clearly visible were landmarked (Table S2; Methods S2).

| Estimating wing shape asymmetry
Differences in wing scale, location and rotation were first removed (Kendall, 1977), by performing a Generalised Procrustes Alignment using the 'geomorph' package (Adams et al., 2020). To assess measurement error, all specimens were landmarked a second time, and the repeat measurements were subjected to a Procrustes analysis of variance (ANOVA; Side × Individual) for each species and sex (see Methods S3). We found no evidence of significant measurement error (Procrustes ANOVA; mean squares for measurement error <0.0002 relative to ≥0.00156 for Side × Individual across all species; Table S4; following Friedli et al., 2020).

| Dealing with low-quality images, angled wings and a skewed distribution of specimens
Tilt of the specimens relative to the position of the camera, image pixelation, differences in specimen illumination and specimens with warped wings can all cause error in measuring FA (Webster & Sheets, 2010).  We therefore undertook two filtering steps: (a) we compared the wing shape of all specimens (of both their left and right forewings) to the mean wing shape across all specimens (using Procrustes distance) and removed all specimens that had a Procrustes distance above the upper quartile of this distribution; (b) we removed specimens with a wingangle differential (i.e. absolute difference between the left and right forewing angle) that was larger than the upper quartile plus the interquartile range (Methods S4; Table S5; Figure S4). This was because a large difference between the positioning (i.e. angles) of a specimen's wings (before Procrustes alignment) correlated significantly with higher wing shape FA in all species (all r > 0.18, p < 0.001; Figure S5).
Therefore, by filtering out specimens with high FA, our findings and effect sizes can be considered conservative.
We then implemented a third filtering step where we rarefied the data to ensure an even distribution of specimens across the 20th century. First, we reduced the dataset to contain a single randomly chosen specimen for each species (n = 4), month (n = 9), year (n = 93), climate region (n = 9) and caste (n = 3). Second, we limited each year to a maximum of 20 randomly selected specimens (as some years contained a disproportionately large number of specimens), to prevent heavily weighted years biasing trends over time, leaving 590 specimens for our analyses (Table S2; Figure 2).

| Statistical analyses
We investigated how FA (dependent variable): (1) differed between the four bumblebee species when pooling all specimens across the century; (2a) responded between the two halves of the 20th century and (2b) whether the trend over the 20th century was nonlinear as a function of time; and (3) has responded nonlinearly with the combined effect of mean annual temperature and annual precipitation. For each of these three questions, we first fitted 'maximal' models containing all biologically relevant variables, and then simplified these maximal models by removing terms and comparing the reduced models using ANOVAs (with model variances compared against the F distribution); we additionally used these F-tests to assess the significance of model terms and report them in Section 3.
We additionally explored comparing reduced models using AIC values, and the results were qualitatively similar.
To examine how FA differed between the four species (1) (South, 2017), 'ggspatial' (Dunnington & Thorne, 2020) and 'ggplot2' (Wickham et al., 2019) were used to produce this image. function, with a separate smooth for each species (controlling for the specimens being collected from different locations across Britain); and year (treating each year as a separate level within a factor) as a random effect (to control for the effects of FA changing over time). After model selection, the fitted GAMM included: caste; latitude and longitude, as a global smooth; and year as a random effect.
To examine how FA changed between the two halves of the century (2a), we divided the century into two periods (1900-1949 and 1950-2000), and modelled the interaction between century half and species identity using a generalised additive model (GAM).
The maximal model additionally included caste as a covariate (and an interaction between species and caste), and latitude and longitude year as a separate level within a factor) as a random effect. The fitted GAMM included: caste; mean annual temperature and precipitation combined in a global 2D tensor product smooth; latitude and longitude, as a global smooth; and year as a random effect.
Model diagnostic plots were checked for normally distributed residuals and homogeneity of variance. FA was log-transformed to ensure model residuals were normally distributed. The fitted models were visualised using the 'predict.gam' function (Wood, 2012) and 'ggplot2' (Wickham et al., 2019).

| Which bumblebee species show the highest levels of fluctuating asymmetry?
Wing shape FA indicated differences between the four bumblebee

| Change in fluctuating asymmetry over the 20th century
When investigating how FA changed between the first half of the century

| Changes in fluctuating asymmetry with environmental variables
When investigating how FA changes in response to the combined effect of mean annual temperature and annual precipitation, we found the trends (smooth estimates) did not vary between species (F-test: model DF for GAMMs with and without species-specific smooths = 60.8 and 48.0, F = 1.09, p = 0.354; Tables S16-S19), indicating the four species responded similarly. Indeed, we found a significant interaction between the two variables for bumblebees as a whole (GAMM: EDF = 8.99, F = 2.33, p = 0.00806; adjusted R 2 = 17.7%; Table S19; Figure 5

| Bumblebee species likely differ in their baseline levels of fluctuating asymmetry
Our results indicate that the degree of baseline FA can vary even between relatively closely related species, with lower levels of FA in B.

F I G U R E 4
Nonlinear trend in wing shape fluctuating asymmetry (FA) over the 20th century. Points represent the raw FA data for all species, with different shapes for the different bumblebee species. The solid black line shows the generalised additive model (GAM) prediction for all species (based on a worker bumblebee at mean latitude and longitude, as the model prediction accounts for variation in FA across bumblebee castes and sample collection locations). Shading shows 2 × standard error of the GAM estimates. Note that the y-axis is on a log scale. Sample size of B. hortorum = 133, B. lapidarius = 179, B. muscorum = 75 and B. pascuorum = 203 specimens.
Intriguingly, the lower levels of FA in B. muscorum are seemingly at odds with its status as a vulnerable species (Rasmont et al., 2015) and a relative 'loser' from environmental change, as evidenced by large range contractions in Britain since the 1960s (Casey et al., 2015).
However, given the rarity of B. muscorum in the latter half of the century, the low levels of FA could be explained by B. muscorum being primarily confined to its favoured habitat (mainly coastal regions, with tall, flower-rich grasslands; Edwards, 2012). Consequently, low FA may be an artefact of natural history collectors being unable to sample specimens from locations where B. muscorum would be rarer and possibly subject to greater environmental stressors (e.g. Adamski & Witkowski, 2002). Using FA as a possible signature of stress may therefore be more useful for detecting early-warning signs for generalist species experiencing a wider range of conditions (i.e. a larger geographic distribution) than for more specialist species with narrower habitat requirements. The higher baseline levels of FA in B. pascuorum and B. lapidarius are also interesting, as these are the two more common species and both have recently expanded their ranges northwards (Casey et al., 2015). Higher FA could therefore be a potential signature of stress when colonising new environmental and climatic conditions (Krause et al., 2016;Liebl & Martin, 2013).
Furthermore, this result suggests that, although FA may be an indicator of stress in populations (Beasley et al., 2013;Clarke, 1995;Leary & Allendorf, 1989;Parsons, 1992;Van Dongen, 2006), high FA levels in a population may not always equate to negative population outcomes (reviewed in Clarke, 1995), and, indeed, not all studies have detected a link between high FA and stress (e.g. Bjorksten et al., 2001;Floate & Fox, 2000;Servia et al., 2004;Ward et al., 1998). Previously, this has been suggested to be due to studies having small sample sizes (Babbitt et al., 2006), large measurement error (Palmer, 1994) or measuring traits less susceptible to FA (Arundell et al., 2019;Lazić et al., 2013). However, the lack of association between FA and fitness here may be as result of the most resilient individuals being able to bear higher stress levels (similar to Zahavi's (1975) handicap principle).

| Signatures of higher fluctuating asymmetry in the latter half of the 20th century
When investigating trends over the century, FA showed a relative increase for all four species after c. 1925, with FA levels at the end of F I G U R E 5 Bumblebee wing shape fluctuating asymmetry (FA) changes nonlinearly with the combined effect of mean annual temperature and total annual precipitation. Heatmap colours represent GAMM predictions of FA for workers sampled at mean latitude and longitude, as the GAMM used to predict FA accounted for variation across bumblebee castes and sample collection locations. Darker colours represent higher predicted ('Pred.') FA. Points indicate the raw ('Obs.') FA values, with the size indicating the magnitude and different shapes representing the four species. To visualise what wings look like with relatively low (Procrustes distance = 0.010) and high (0.065) FA, two sets of wing outlines are placed either side of the scale bar, representing the lower and upper predicted asymmetry ranges respectively. Grey lines with blue points represent right forewings, and black dashed lines with black points represent left forewings (superimposed onto the right wings using Procrustes alignment). Wing outlines were created using a modified version of the 'plotRefToTarget' function from the 'geomorph' package (Adams et al., 2020). the century higher than pre-1925 levels. Although drivers of stress are multifactorial and complex (Zaragoza-Trello et al., 2021), this trend of increasing FA coincides with large changes in agricultural practices, driven by government policies, that occurred in Britain following the First World War and during and after World War II, including expanding arable land and increasing use of pesticides and herbicides (Ollerton et al., 2014;Robinson & Sutherland, 2002). Such practices have individually been shown to increase FA in laboratory and field experiments in a wide range of invertebrates (Abaga et al., 2011;Chang et al., 2007;Costa & Nomura, 2016;Friedli et al., 2020;Mpho et al., 2001) and vertebrates (Coda et al., 2016).
Alternatively, other global change drivers, such as emergent diseases (Fürst et al., 2014;Yordanova et al., 2022), spread of invasive species (Meeus et al., 2011) and climate change (Forister et al., 2018;Halsch et al., 2021;Kerr et al., 2015), may also explain increasing FA over the latter part of the century. Additionally, the past century has seen possible increased competition between managed honeybees and wild bumblebee species (Thomson, 2016). Although to our knowledge no studies have assessed the link between bumblebee interspecific competition and FA, competition between larvae is associated with higher levels of FA in butterflies (e.g. Parage aegeria; Gibbs & Breuker, 2006).
Indeed, a growing body of evidence supports the idea that environmental stressors can impact individuals not just as adults but from the start of their life cycle, that is, during larval development (Pellegroms et al., 2009;Smith et al., 2020). But regardless of the exact mechanism, quantifying FA has potentially given insights into historic stress, and specifically provided a more temporally accurate assessment of when this may have been experienced in wild bumblebee populations.
Despite the consistent trend across species over time, we note that unexplained variance in FA remains. This could be due to differences in local habitat or microclimate (Henriques & Cornelissen, 2019;Kark, 2001;Schmeller et al., 2011) or other stressful conditions, which could interact in non-additive ways (Zaragoza-Trello et al., 2021).
Alternatively, variance in FA could also be influenced by the positioning of wings on museum specimens. However, as FA was not higher in older specimens, we can be confident that the increase in FA over time is not a result of damage to specimens within individual collections.
Ideally, FA would be measured on wings that have been detached from the specimen and slide mounted. However, as we are dealing with irreplaceable biological artefacts, a challenge for future studies would be to develop high-throughput scanning methods to more accurately measure FA from wings of pinned museum specimens (e.g. Perrard et al., 2012;Plum & Labonte, 2021).

| Warm and wet years are associated with high fluctuating asymmetry
FA was consistently higher in warmer and wetter years. Importantly, this finding seems to support other studies finding a strong signal of climate in driving insect population trends (Forister et al., 2018;Halsch et al., 2021;Kerr et al., 2015;Román-Palacios & Wiens, 2020; but see Guzman et al., 2021). It is also consistent with most bumblebee populations being better adapted to colder conditions, with different populations-regardless of whether they are from warm or cold environments-having similar tolerances to high temperatures (Martinet et al., 2021;Pimsler et al., 2020). Moreover, temperature stress has been linked to increased FA across a range of invertebrate species (Gerard et al., 2018;Imasheva et al., 1997;Nishizaki et al., 2015;Vishalakshi & Singh, 2008). Additionally, warmer, wetter years could promote pathogen spread and infection (Lafferty, 2009;Neidel et al., 2017), which are also associated with increased levels of FA (Arundell et al., 2019). Given that mean annual temperature for many regions will likely increase under climate change, our findings could reveal an early-warning sign, with some species potentially being unable to sustain continued increases in FA. Alternatively, these climatic conditions could correlate with increased nutritional stress, as declines in floral abundance, and consequent loss of insect pollinators, have been reported with climate extremes (Høye et al., 2013;Iserbyt & Rasmont, 2012;Thomson, 2016). Expanding this analysis, and network of museums, to include more bumblebee species and regions representing the southern and northern limits of species' ranges that may experience more severe climate change (such as the Mediterranean and Arctic regions in Europe) is therefore a research priority, to understand the extent of impact that natural populations will face in the future.

| CON CLUS IONS
By leveraging morphological data in museum specimens, this study has, for the first time, tracked FA over an extended period and at a temporal resolution that potentially allows more accurate assessment of when during the century bumblebee populations may have experienced stressors. By using wing shape FA (rather than other measures of a population's response, such as geographic distribution), we avoided potential lag effects of stressors that often obscure how populations are responding at a given time point. Although the relationship between FA and stress is not solely determined by environmental stressors, our results have potentially revealed periods of historical signatures of developmental stress worthy of further investigation. Moreover, across our studied species, increases in FA over the 20th century were found in warmer, wetter years. Moving forward, measuring FA in contemporary specimens, and especially relating it to land-use and climate change hotspots, may provide a mechanistic insight into the causes of stress and provide earlywarning indicators to inform safeguarding strategies to protect beneficial insects. Ultimately, our study highlights the crucial role long-term collections of museum specimens can play in understanding past and future pollinator responses to environmental change.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available from the Environmental Information Data Centre (EIDC) repository https://doi.org/10.5285/26965 35e-564a-4c6a-877e-51599 6fa97a1 (Arce et al., 2022).