Impact of crop exposure and agricultural intensification on the phenotypic variation of bees

In a context of rapid global change, understanding how environmental stressors can impact phenotypic variation, and which phenotypic traits are predominantly affected can be particularly relevant. Indeed, potential phenotypic modifications could affect the functionality of traits from taxa that are in decline but that are keystone species in many ecosystems. In this study, we assessed the impact of environmental drivers and agricultural intensification on two crucial pollinator species: the honeybee ( Apis mellifera ) and the buff-tailed bumblebee ( Bombus terrestris ). Among eight countries representing four major European biogeographical regions [i.e., Boreal (Sweden and Estonia), Atlantic (Ireland and United Kingdom), Continental (Germany and Switzerland) and Mediterranean (Spain and Italy)] and two type of crops (i.e., apple orchards and oilseed rape) we assessed how landscape structure, latitude and pesticide management could impact their wing morphology. Two sampling sessions were conducted: the first one when the hives/nests were settled on the field (T0) and a second sampling session after the potential effect of agricultural intensification (T1). Using a dataset of more than 7238 wings, we measured the wing size, shape and asymmetry. We observed that, in several countries, a shift in most of the morphological traits occurred between T0 and T1. When focusing on the drivers of phenotypic variation in T1, the levels of significance for some potential drivers were sometimes high, while most of the variation remained unexplained. The latitude and, more rarely, grassland cover were found to partly explain the wing modifications. In light of these results, we conclude that phenotypic shifts can occur in a very short period, after encountering new field conditions. Further studies should be conducted to better understand which alternative drivers could explain morphological changes in the agro-ecosystem after crop exposition, as well as the potential consequences of these changes on foraging performance or pollination efficiency.


Introduction
High intraspecific phenotypic variability can be beneficial and can increase survival odds in fluctuating environments (Forsman and Wennersten, 2016;Maynard et al., 2019).Phenotypic variation can be induced by many factors at different levels from molecular, genetic, cellular, individual to populational levels (Willemore et al., 2007).In fluctuating environmental conditions, plasticity is a common response of an organism and is, by definition, a non-evolutionary response since it involves one generation of an individual (Gibert et al., 2019).While mechanisms such as phenotypic plasticity can increase phenotypic variability, others can help to produce the same phenotypes despite stresses outside the normal environmental bounds usually experienced by a population, such as canalisation and developmental stability (Waddington, 1942;Debat & Le Rouzic, 2018).The concept of canalisation is closely related to the concept of phenotypic plasticity, as canalization describes the ability to produce a same phenotype despite different environmental/genetic conditions, while phenotypic plasticity describes the ability to vary and express different phenotypes depending in changing conditions.Developmental stability also describes a developmental process that ensure phenotypic consistency.However, conversely to canalisation, it is measured under similar environmental and genetic conditions, despite developmental noise (Waddington, 1957;Debat and David, 2001).To ensure similar environmental and genetic conditions, developmental stability is thus measured using both sides of symmetrical organisms; e.g., bilateral animals: deviations from perfect symmetry between each side reflect random variation during development (Palmer and Strobeck, 1986).Deviation from perfect symmetry is called fluctuating asymmetry (FA); an increase in FA indicates a decrease in developmental stability.
There is strong interest for measuring morphological changes in response to environmental factors.It is often assumed that stressful conditions will lead to an increase of phenotypic variance, for example, the disruption of canalisation mechanisms by stressors can lead to an increased variance of the observed phenotypes, also known as the "decanalization" of a trait (Hoffmann and Parsons, 1997;Lack et al., 2016).Increases of FA have often been tested as a potential indicator of stressful conditions although the results are contrasted, even for the same species facing the same stressor (Kruuk, Slate, Pemberton & Clutton-Brock, 2003;Rasmuson, 2002;Beasley et al., 2013).Many debates emerged about FA sensitivity as a stress indicator leading to the conclusion that, at least, it cannot be generalised (Palmer, 1996;Houle, 1998;Debat, 2016).Moreover, shift of overall shape and size could also be induced by genetic or environmental drivers even if the use of shape variation as an indicator of stress has been less studied (Hoffmann et al., 2002;2005;Pinto et al., 2015, Dellicour et al., 2017;Gérard et al., 2021).
Insect wings are a particularly accurate means to assess phenotypic variation because they are flat, two-dimensional, and present many homologous landmarks.They are highly diversified across taxonomic orders, and sometimes even within populations (Grimaldi and Engels, 2005;Dehon et al., 2019).Wings are considered as crucial functional traits for flight efficiency (Wooton, 1992;Wakeling and Ellington, 1997), foraging, and dispersal capacities (Bots et al., 2009;Johansson et al., 2009).Their morphologies experience adaptive selective pressures (Kölliker-Ott et al., 2003), and can be driven by both genetic and environmental factors such as temperature (Morales et al., 2010), habitat fragmentation (Outomuro et al., 2013) or urbanization (Banaszak-Cibicka and Zmihorski, 2012).Indeed, long and slender wings are more efficient for long distance flight and may be more adapted in more fragmented habitats (Betts and Wootton, 1988;DeVries, Penz & Hill, 2010).Among insects, several studies highlighted that stressful conditions like high developmental temperatures or toxins can lead to changes in wing shape, size, or FA (e.g.Dingha et al., 2004;Beasley et al., 2013;Gérard et al., 2018).Consequently, these changes may affect foraging behaviour.For example, changes in wing morphology can affect the time that bees are spending in each foraging trip (Foster & Carter, 2011), which may ultimately affect pollination efficiency.
In this study, we gathered an extensive dataset from 128 sites across eight countries (Estonia, Germany, Ireland, Italy, Spain, Sweden, Switzerland and UK), two type of crops (apple and oilseed rapes), along a gradient of landscape structure (edge density, total grassland area, total urban area, etc.) and pesticide intensity during two sampling sessions.We studied the impact of these potential drivers of wing morphology in two bee species, the honeybee (Apis mellifera) and the buff-tailed bumblebee (Bombus terrestris), gathering a dataset of 3619 individuals.
We assessed the impact of pesticide management, crop type and environmental predictors (latitude, landscape structure) on phenotypic variation.First, we hypothesized that there would be significant differences in morphological traits between specimens collected in the first sampling session (when the hives/nests were placed in the field) and the second sampling session (i.e., once crop blooming occurred and after potential exposure to agricultural intensification).Indeed, specimens collected during the first sampling session developed before being exposed to any potential stressors, while specimens collected during the second sampling session could have experienced pesticide exposure, different landscape structure among different climates as well as different type of crops.Thus, we also hypothesized that, in fields with more intense agricultural practices, there would be a decrease of wing size, an increase of FA, a modification of wing shape and an overall increase of trait variance because of the increase of stressful conditions.In the context of global change and potential negative population trend, detecting the drivers of phenotypic changes is of great importance.Indeed, if we could identify a phenotypic trait that is consistently impacted by harsher anthropic pressures, it may be used as a stress indicator in future studies on bees.Phenotypic changes may also fundamentally affect the functionality of an organism, which is even more crucial for pollinating species.It is difficult to directly infer that changes observed under controlled conditions are equivalent to those that would appear in field conditions (Gibert et al., 2019), hence the need to lead large-scale studies testing for the impact of various conditions on phenotype.

Biological models and sampling
The species used in the study were the honeybee Apis mellifera and the buff-tailed bumblebee Bombus terrestris, which are both part of the Apidae family (Hymenoptera) and social species that produce workers, males and queens (Michez et al., 2019).We focussed on measuring the wing morphology of workers because they collect most of the resources for the colony and because wing morphology can differ between sexes and castes (Gérard et al., 2015).Specimens of A. mellifera and B. terrestris were sampled across eight countries representing four major European biogeographical regions: Atlantic (Ireland and United Kingdom), Boreal (Estonia and Sweden), Continental (Germany and Switzerland) and Mediterranean (Italy and Spain).Three colonies of each species were placed at each sampling site, each of the colonies was at least 5 m apart to avoid interference among them.All bumblebee colonies were purchased from local suppliers, provided by the company Biobest (Westerlo, Belgium), but they were derived from the same commercial population representing the local wild subspecies.For 6 of the 8 countries this was a single continental B. terrestris population, but for Ireland and the UK (where the use of non-native subspecies is restricted) it was a single B. terrestris audax population, to match the subspecies found in these islands.Before they were set up in the field, every bumblebee colony contained approximately the same number of workers and were supplied with a sugar syrup (BIOGLUC®, biobest) which was removed once the bees were on the field (see Hodge and Stout, 2019 for further details about the protocol before nests were settled on the crops).Apis mellifera colonies were provided by local suppliers, representing local races (Hodge et al., in prep; Table S1 for details about the honeybee subspecies).Every hive has been checked before the start of the experiments so that they were free of diseases, a normal strength for the season and the location, were similar in size, etc (see Hodge and Stout, 2019 for further details about the protocol before nests were settled on the crops).Every hives and colonies were set up on the crops one to two weeks before the blooming.In each country, two crops with contrasting growing systems and pollination biology were sampled: apple (Malus domestica), a perennial, self-incompatible crop grown in long standing (5-10 year) orchard plantations and winter-sown oilseed rape (Brassica napus), an annual, self-compatible crop planted irregularly as a break crop in arable rotations.There were eight sampling sites for each type of crop in each country.Thus, a total of 128 sites were selected along a gradient of land use intensity see Hodge and Stout (2019) for further details about the site selection).All sites were surveyed using the same methods as part of the European Union POSHBEE project (http://poshbee.eu),standardised following the PoshBee protocols (Hodge and Stout, 2019).
Two standardised sampling sessions were carried out across all sites: the first session (T0 sampling) occurred when the hives/nests were settled on the sites (from late March to mid-May 2019, according to the countries) and another session (T1 sampling) once the blooming occurred and after potential effect of agricultural intensification (from mid-May to late June 2019, according to the country).Every site has been sampled once, on a same day for T0, then on a same day for T1.On average, there were 45.1 days between the T0 sampling and the T1 sampling (mean value -see table S2 for details country per country).Eight workers from each species were selected for each sampling session (T0 and T1) for further analyses, leading to an overall data set of 2048 specimens of each species, coming from 384 colonies (1024 specimens collected in T0 and 1024 specimens collected in T1).Some of these individuals were excluded from further analyses due to wing damage, or because not enough specimens were collected (see Table 1 for the final dataset and Vanderplanck et al., 2021 for a detailed description).

Potential drivers of wing morphology
Landscape structure was quantified within a 1-km radius around each site by calculating landscape metrics derived from manually digitized habitat maps.A total of ten land cover categories including bare areas, crops, grasslands, heathlands, orchards, roads, surface running waters, urban areas, waterbodies, and wetlands were identified and mapped using Geographical Information Systems Software (ArcGIS Pro, 2.4.1.ESRI) and high-resolution images provided by World Imagery (ESRI).First, we calculated the proportion (%) of all ten land cover features for each landscape.Then we defined the landscape intensity gradient (LIG) as the sum of proportion of crops and orchards per landscape.Landscape structure was assessed by calculating (i) the Shannon's Diversity Index (SHDI) at the landscape level with all 10 land cover categories, (ii) the edge density of semi-natural habitats (ED) by dividing the edge length of semi-natural habitats by the total area of the corresponding habitat map and (iii) the total area (CA) of two land cover categories: grasslands and urban areas.All landscape metrics were calculated using the software FRAGSTATS 3.3.(McGarigal et al., 2012).
Pesticide data were collected directly from participating farmers using a standardised questionnaire hosted on Qualtrics and translated into local languages.Alongside demographic questions, the questionnaire asked, (i) which plant protection product (PPP) they applied, (ii) for each product which dates they applied and (iii) for each date of application, what was the total volume applied (l/ha).Because some PPPs are mineral, they were applied as kg/ha but, in our analysis, this is equivalent to volume.The total volume of each PPP individually and all PPPs in total was then calculated for each site to give a measure of pesticide pressure.All participants were asked if this regime was representative of the management of their other fields and replied that it was, indicating that these responses are likely to be representative of the wider landscapes around the survey sites.The questionnaire was checked by representatives of farming organizations for clarity of language prior to dissemination and approved by the University of Reading Ethics Committee.Table S3 summarizes the response rate.Note that not all respondents supplied complete data.

Wing size and shape data acquisition
Specimens from A. mellifera and B. terrestris were analysed separately due to substantial differences in wing morphology.The right and the left forewings of each individual were cut off at the base, flattened and photographed using an Olympus SZH10 microscope coupled with a Nikon D70 Camera (Shinjuku, Japan).Each wing was digitized using a set of two-dimensional coordinates of 18 landmarks (Fig. 1) from right and mirror-reflected left wings photographs using tps-DIG v2.32 (Rohlf, 2013).While both wings of each individual were used for fluctuating asymmetry (FA) analysis, only the right wing was used for wing size and shape analyses.Using the Generalized Procrustes Analysis superimposition, we scaled, translated, and rotated the landmark configurations against the consensus configuration to remove all non-shape components for further analyses (Bookstein, 1991;geomorph package, Adams and Otarola-Castillo, 2013).Centroid size was used to estimate wing size: it corresponds to the square root of the sum of squared distance between each landmark and the centroid of the landmark configuration and is also a good proxy for body size in insects (Bai et al., 2016;Gérard et al., 2018).To account for measurement error, a subset of 100 wings was digitized twice.We calculated individual vectors of size FA by subtracting the centroid size of the right and left wings of each individual and selecting the absolute value of this subtraction.We then assessed individual vectors of shape FA by calculating the square root of the sum of each squared value of each landmark (x and y).

Statistical analyses
First, we assessed collinearity between the potential drivers of morphological variation and removed variables that had a Pearson's correlation coefficient higher than r = 0.7 (Dormann et al., 2013).Six final variables were selected for inclusion in the models: latitude, the area of grassland, the cover of urban areas, pesticide use, crop type and country (Supplementary material Fig. S1).
All the following analyses were carried out first for A. mellifera then for B. terrestris.We assessed if any significant differences were observed between specimens collected at T0 and T1 among the different countries based on their (i) wing size, (ii) wing shape, (iii) size FA and (iv) shape FA.After checking for residuals normality, we used three separate linear mixed models (LMM, lmer4 R package) for wing size, size FA and shape FA as well as Post-Hoc comparisons (emmeans function, emmeans R package), and a Procrustes ANOVA (procD.lmfunction; package geomorph, Adams and Otarola-Castillo, 2013) for wing shape.The procD.lm function allowed us to use a multidimensional response variable in a linear model, which is crucial when using landmark configurations.In  these models, the sites were included as random variable.We also assessed if any differences in the variance of size, size FA and shape FA occurred in each country between the T0 and T1 sampling using F-tests (var.testfunction; package stats).
We then assessed the effect of environmental drivers on phenotype.As a preliminary analysis, we checked if any differences existed among the morphological traits within the different countries in T0.As no significant differences were observed (except for the wing size of A. mellifera from Ireland, which are not used in T1), we directly assessed the differences among specimens collected in T1.As different sites necessarily used different subspecies, we also assessed if any significant morphological differences existed in T0 among A. mellifera or B. terrestris subspecies.The method was similar to the analyses described below and described in greater detail in the Supplementary Material (Appendix 1-2).To better understand wing size changes of specimens collected at T1, the effects of environmental drivers on centroid size were analysed using linear models (LM) and linear mixed effects models with Gaussian distribution (LMM) using the lmer4 R package.After checking for residuals normality (Shapiro test), we fitted a full mixed linear model (LMM) with wing size as the response variable, country and site (i.e. to account for pseudoreplication) as random factors and each predictor as fixed factors, as well as all two-way interactions.We selected the best model using the lowest AIC (Burnham and Anderson, 2004) after testing for all LM and LMM model combinations, and depending on the amount of variance explained by the random factors (i.e., country and site).We then tested the impact of the potential drivers on size and shape FA using the same statistical procedure as described above for centroid size data.
To analyse the impact of the selected drivers on wing shape of specimens collected at T1, we first assessed the potential allometric component of the wing shape variation.Specifically, we conducted a Procrustes ANOVA (Goodall, 1991) with permutation procedures to assess if centroid size (response variable) significantly explained wing shape variation (explanatory variable) using the procD.lmfunction (package geomorph, Adams and Otarola-Castillo, 2013).To understand the effects of environmental drivers on wing shape, we then analysed the landmark configurations by fitting a linear model using the function lm.rrpp (package RRPP, Collyer and Adams, 2018).This function fits a linear model using an ordinary least squares (OLS) estimation of coefficients on multidimensional data (i.e., landmark configurations), and uses a randomized residual permutation method.This method is also appropriate for including a random effect in models using multidimensional data, which was crucial for including country and site as random factors.We fitted a full linear model with landmark configurations as the response variable, each predictor as fixed factors as well as all two-way interactions, and country as a random factor.We then selected the best model based on the function model.comparaison(package RRPP, Collyer and Adams, 2018), which calculates the log-likelihood, the trace of the residual covariance matrix and AIC of each tested model.When the random factors were not included in the model to assess wing shape variation, we used the procD.lmfunction mentioned previously.Finally, to test the potential effects of subspecies on the morphological changes observed in T1, we used the same statistical method presented above, and described the results in Supplementary Material (Appendix 1 & 3).Each analysis was conducted in the R software version 4.1.1.(R Core Team, 2020).
Data on landmark coordinates and the potential drivers of wing morphology are compiled in Vanderplanck et al. (2021).

Honeybee wing analyses between T0 and T1
Honeybee specimens collected at T1 had significantly larger wings than specimens from T0 (p = 0.03) and the random factor site explained 7% of the variance that was left in the residuals after accounting for the variance explained by the fixed effects.More specifically, we observed that specimens from Estonia (p = 0.003), Sweden (p < 0.001) and UK (p = 0.008) had significantly larger wings in T1 than in T0, but we found no significant differences in other countries (all p-values > 0.2) (Fig. S2A).No significant increase of wing size variance was observed between T0 and T1 in any country (all p-values > 0.1).There was a significant difference in wing shape between specimens collected in T0 and T1 (p < 0.001).We also observed significant increase of size FA in T1 (p < 0.001) as well as different shape FA (p < 0.001).The random factor site explained respectively 8.6% and 10% of the variance that was left in the residuals after accounting for the variance explained by the fixed effects.More specifically, when looking within each country, the only significant difference in size FA was in Estonia, where it significantly increased in T1 (p < 0.001) (Fig. S3A).Between T0 and T1, size FA variance significantly increased in Estonia (p < 0.001), Germany (p = 0.04) and UK (p < 0.001) and significantly decreased in Italy (p < 0.001) but was not significantly different in other countries (all pvalues > 0.1) (Fig. S3A).Shape FA increased in each country (all pvalues < 0.009), except in Germany (p = 0.82) and in Switzerland (p = 0.98) (Fig. S4A).Finally, shape FA variance significantly increased in Estonia (p = 0.003), Germany (p = 0.003), Spain (p < 0.001), Switzerland (p = 0.01) and UK (p < 0.001) but was not significantly different in other countries (p > 0.3 for both) (Fig. S4).

Drivers of honeybee wing morphology in T1
The best model (next best model ΔAIC 6.3, Table S4-5) that explained wing size of honeybees included latitude, crop type, and their interaction (adjusted R-squared = 0.04).The interaction of latitude with crop showed that honeybees collected in apple orchards from higher latitudes had larger wings (p < 0.001; Table S4; Fig. 2A).We observed a significant increase in wing size with latitude (p < 0.001; Table S4; Fig. 2A), and the wing size of specimens collected on apple orchards was significantly larger than for oilseed rape specimens (p = 0.004).
The best model (next best model ΔAIC 74.96, Table S6-7) that explained honeybee wing shape included latitude as fixed factor and country as random factor (R-squared = 0.15; p < 0.001).We observed a significant impact of country (p = 0.001, R-squared = 0.11; Table S6) and a slightly significant impact of latitude (p = 0.035, R-squared = 0.04; Table S6) on wing shape.The best model (next best model ΔAIC 16.79, Supplementary material, Table S8-9) to explain size FA variation of honeybees included latitude, crop and their interaction (Adjusted R-squared = 0.08).We observed a significant increase of size FA with latitude (p < 0.001), and among specimens collected in apple orchards (p < 0.001).The interaction of latitude with crop showed that honeybees collected in apple orchards from higher latitudes had higher size FA (Table S8; p < 0.001).The best model (next best model ΔAIC 26.13, Table S10-11) that explained shape FA variation of honeybees included latitude, crop and their interaction (Adjusted R-squared: 0.06).Shape FA increased significantly with latitude (p < 0.001) and was higher specimens collected in apple orchards (p < 0.001).The latitude/crop interaction also significantly impacted wing shape (p < 0.001): specimens collected in higher latitudes from apple orchards had significantly larger shape FA.

Bumblebee wing analyses between T0 and T1
Bumblebee specimens collected in T1 had, overall, significantly smaller wings than specimens from T0 (p < 0.001) and the random factor site explained 8% of the variance that was left in the residuals after accounting for the variance explained by the fixed effects.More specifically, we observed that specimens from Germany, Sweden and UK (all p-values < 0.001) had significantly smaller wings in T1 than in T0, but no significant differences were found in the other countries (p > 0.5) (Fig. S2B).Between T0 and T1, wing size variance significantly increased in Germany (p < 0.001), Spain (p < 0.001), Sweden (p = 0.03) and Switzerland (p = 0.003) (Fig. S2B) but was not significantly different in Italy (p = 0.1) or the UK (p = 0.18).We observed a significant difference in bumblebee wing shape (p = 0.001).Finally, between the two sampling sessions, we observed no significant differences of size FA (p = 0.207) but significant differences in shape FA (p < 0.001).The random factor site was explaining respectively < 0.001% and 5% of the variance that was left in the residuals after accounting for the variance explained by the fixed effects.There were significant differences in shape FA in each country (p < 0.004) except in Switzerland (p = 0.99) (Fig. S3B).Size FA and shape FA variances did significantly change in all countries (all p-values > 0.1) (Fig. S3B, S4B).

Drivers of bumblebee wing morphology in T1
The best model (next best model ΔAIC 2.03, Table S12-13) that explained wing size variation of bumblebees included latitude, crop type, the area of grassland, and every interaction between these three fixed factors (Adjusted R-squared = 0.18).Wing size decreased slightly with latitude but the trend was not significant (p = 0.092; Fig. 2B).Wing size was slightly larger in specimens collected in apple orchards, although the relationship was not significant neither (p = 0.059; Fig. 2B).Wing size increased significantly with the area of grassland (p = 0.039).We also detected significant effects for the interaction between the area of grassland and crop type (p = 0.019): the wing size of specimens collected in apple orchards with a larger area of grassland was larger.The interactions between the other fixed factors were not significant (all p-values > 0.05).
The best model (Table S16) that explained bumblebee wing shape included country only (r-squared = 0.06; p = 0.001).All countries were significantly different from each other (all p-values < 0.01).
The best model (next best model ΔAIC 0.5, Table S15-16) that explained size FA variation of bumblebees only included latitude, but the model was not significant (p = 0.37).The best model (next best model ΔAIC = 1.99,Table S17-18) that explained shape FA variation in bumblebees included latitude, crop type, the area of grassland, and every interaction between these three fixed factors (Adjusted r-squared = 0.07).Shape FA significantly decreased as the area of grassland decreased (p < 0.001).Shape FA was also significantly higher for specimens collected in smaller grassland area from lower latitudes as shown by the significant interaction between the area of grassland and the latitude (p < 0.001).

Discussion
Using an extensive dataset collected across a European landscape gradient, we explored and identified environment-related changes in phenotypic variation for both honeybees and bumblebees.We demonstrated that shifts in morphological traits could occur quickly when experiencing new environmental conditions (e.g., crop exposure during the full blooming).Wing size, wing shape, size fluctuating asymmetry (FA) and shape FA of honeybees and bumblebees were shown to mainly vary between the two sampling sessions (T0 versus T1).Most of the observed changes were in accordance with our hypotheses: size decreased in the bumblebees in the second session (T1), size FA increased in the honeybees in T1, and shape FA increased both in the honeybees and bumblebees in T1.However, no changes occurred in the size FA of bumblebees, and size FA only increased in one country in the honeybees in T1.We also observed that trait variance increased regularly between T0 and T1 among the different countries.However, even if the levels of significance for potential morphological drivers were sometimes high, most of the variation in shape and size remained unexplained by the selected drivers in most of the models.Variation in morphological traits was often significantly impacted by latitude, the type of crop (in honeybees) or the interaction of latitude and the type of crop but pesticides management and landscape composition largely failed to explain the observed modifications (summarized in Table 2).The absence of any effect of pesticide use on phenotype may be due to how pesticide management was measured in our study.While we used the volume of plant protection product (PPP) applied as a proxy for pesticide management, their specific active ingredients, real occurrence in pollen and/or nectar, and interaction effect were not measured, and may well be more important to assess the impact of pesticides on phenotypes.Similarly, our habitat maps did not distinguish between different types of grasslands, which may explain the weak relationships between landscape composition and phenotype.With no distinction made between intensive grasslands and natural grasslands, our dataset could not discriminate grassland habitats that are favourable for pollinators.Some unexplored drivers could also explain the differences between T0 and T1 within a same country.For example, crop exposure could lead to pathogen spillover, notably from honeybees to bumblebees (Nanetti et al., 2021).Indeed, honeybees may act as super-spreader through plant-pollinator-pathogen networks (Proesmans et al., 2021).The exposure with pathogens could explain the morphological differences between the two sampling sessions, as pathogen infection can lead to phenotypic changes (e.g.Gérard et al., 2018).
Concerning stressful conditions, it is also possible that the conditions experienced by the sampled bees were not sufficiently stressful to induce significant morphological changes.The factors that constitute environmental stress have been debated for decades, leading to different and multiple definitions depending on the field of application (Selye, 1950;Schulte, 2014).For many ecologists, a stressor can be defined as "an environmental factor that, when first applied, impairs Darwinian fitness" (Sibly and Calow, 1989).In this case, stressful factors potentially affect the individual's survival by creating a shift in the allocation of energy, weakening the homeostasis, and ultimately impacting the phenotype (Møller, 2006).In bees, fitness measures can include the number of dead workers, the number of larval cells produced, and the number of ejected larvae or colony growth (e.g., Tasei and Aupinel, 2008;Vanderplanck et al., 2019).While this information has not been collected in this study, future studies could help us to assess if stressful conditions have indeed been suffered by collected specimens.The bumblebees collected in T1 potentially experienced environmental conditions that are outside of the range they usually experience, whereas the specimens sampled in T0 grew in controlled conditions.These new environmental conditions could lead to changes that range from no observable modification of the phenotype if they are completely buffered, to selection of novel morphological traits or extinction of a population, depending on the frequency and the intensity of stressors (Badyaev, 2005).Furthermore, we used two domesticated species that may be particularly capable of buffering different environmental conditions.For instance, in controlled conditions, B. terrestris is especially resistant to hyperthermic stresses, having an extremely high time before heat stupor (Martinet et al., 2021).Moreover, plastic traits also allow species to cope with fluctuating environments; the differences between specimens collected in T0 and T1 could reflect this plasticity.
Previous work has shown that insect size and shape can be strongly impacted by changes in environmental conditions and stressors (e.g., Hoffmann and Shirriffs, 2002), particularly in controlled conditions (Gérard et al., 2018), and along latitudinal clines (Shelomi, 2012).While we did not specifically test for this relationship, body size and wing centroid size are strongly correlated (e.g., Dujardin, 2008;Bai et al., 2016).Although bee size generally increases with latitude in Europe when pooling species of all families together, the decrease of size with latitude can also be observed on several genera like bumblebees (Gérard et al., 2018).Similarly, our results showed that honeybees from higher latitudes had a larger wing size than those from lower latitudes.These trends in size are mostly explained by thermoregulatory mechanisms: being smaller can generally help dissipate temperature in warm environments (Meiri, 2011).However, even if wing centroid size is strongly correlated with body size, it would be interesting to assess if wing size reduction is indeed correlated with body size reduction, or if the allometry of the trait is impacted (i.e. the rate of wing size changes and body size changes are different).In insects, many alternative drivers explain these clines (e.g., food availability, physiological processes linked to metabolism and growth rate, etc; Blanckenhorn and Demont, 2004), but our study did not allow us to specify if one of these drivers could better explain the observed trend.We are not aware of studies focusing on latitudinal clines of honeybees.However, we should not assess the impact of latitude on honeybee wing size alone, because the interaction between crop type and latitude was also significant.Indeed, honeybee wing size was larger at higher latitude but more specifically among apple orchards.This might suggest that honeybees are larger when foraging on apple orchards than on oilseed rape crops, and that apple trees might provide floral resources with higher nutritional quality, or in higher amounts.Indeed, food quality and quantity are among the primary drivers explaining changes in bee body size (Vaudo et al., 2015;Sauthier et al., 2017), though specific analyses would be needed to know if it is truly the case between oilseed rape and apple trees.
Trait variance increases in T1 for both honeybees and bumblebees across several countries, which supports our initial hypothesis.Indeed, stressful conditions can reveal higher phenotypic variability, which may be advantageous, and allow for the selection of new phenotypes and novel adaptations (Bradshaw and Hardwick, 1989;Badyaev, 2005).Phenotypes that were previously canalized (i.e., low variability) can express hidden variation which potentially allows them to deal with new conditions, sometimes caused by the phenomenon called phenotypic decanalization.The clearest evidence of this loss of canalization in wild populations was recently provided for the first time by Lack et al. (2016) and our results corroborate these findings.However, an increase in trait variance can also have potential deleterious effects if the morphological trait moves further away from its functional/favoured optimum, leading to non-adaptive plasticity (Ghalambor et al., 2007).The impact of environmental conditions on wing shape has received less attention (Hoffmann et al., 2002(Hoffmann et al., , 2005;;Gérard et al., 2018). Yet, Hoffmann et al. (2005) suggested that changes in wing shapewhich has also been observed consequently to stressful conditions in laboratory experiment on bumblebees (Gérard et al., 2018) could be an indicator of stress.Both bumblebee and honeybee wing shape differed between T0 and T1; in both cases, the driver that best explained this variation was the country.
Furthermore, the existing literature on the relationships between stressful conditions and FA over the last few decades shows contradicting results and publication biases (e.g., the lack of publication of negative results), which makes it difficult to evaluate (Van Dongen, 2006, 2011).Even if FA has been suggested as a tool to evaluate developmental instability in changing conditions such as found in agroecosystems (Benitez et al., 2020), the relationship is often weak and dependant on the trait, the population, and the species (Lens, Van Dongen, Kark & Matthysen, 2002;De Coster et al., 2013).Despite these discrepancies in the literature, our study showed relatively similar trends both in bumblebees and honeybees: an increase of FA in T1, mostly in shape FA.Our results contradict a previous study where nutritional stress, heat stress, parasitic stress and inbreeding rarely impacted the FA of B. terrestris under controlled conditions (Gérard et al., 2018).However, higher levels of FA are often detected under laboratory conditions, most likely due to positively selected buffering effects in wild populations, or to diluting factors in natural conditions (Beasley et al., 2013).One potential explanation could be that, in contrast to laboratory experiments, pollinators sampled in real

Table 2
Overview of trend between T0 and T1, and the significant impact of the fixed factors and their interactions on the different morphological traits in T1.Only factors being significant at least once are indicated.NA indicates that the driver was not in the model.↑ indicates an increase from T0 to T1, ↓ indicates a decrease from T0 to T1 and X indicates a shift from T0 to T1.Finally, throughout the foraging season, honeybee and bumblebee worker size can increase (Knee and Medler, 1965;Sauthier et al., 2017), though it has not been systematically observed (Couvillon et al., 2010).The drivers of this size increase could be genetically driven, or a side effect of the accumulation of resources during the first weeks of the colony (Lawson et al., 2016;Withee and Rehan, 2016).These phenotypic changes, which can be naturally observed for size, could also occur for the other morphological traits we studied, though this has not to our knowledge been assessed previously.In our study, honeybee worker size increased between T0 and T1 in several countries, thus we cannot exclude that this increase is, for example, genetically driven.However, as we found smaller bumblebees at T1, these previous studies could not explain our observed bumblebee size trend.Moreover, the phenotype responses to the environmental drivers may be related to the subspecies characteristics.Although subspecies was not a powerful predictor of the morphological changes in T1 (Supplementary material; Appendix 3), some differences among both honeybee and bumblebee subspecies in T0 were observed.However, subspecies and country were sometimes confounded in our models, since some subspecies were used only in one country.Because each subspecies may not have the same level of phenotypic plasticity environmental drivers may have stronger effects on specific subspecies (Manfredini et al., 2019).In order to fully disentangle the effects of innate phenotypic plasticity, subspecies and country, future research should also perform genetic analyses.
In conclusion, this study highlighted that shift of morphological traits occurs over a very short period of time.Particularly, trait variance and asymmetry in shape were often higher when exposed to the field conditions.However, while our study highlighted that there is probably a strong environmental basis for the observed phenotypic variation, we were not able to highlight a major driver of this variation in T1.Thus, future studies should aim to identify and investigate the presence of alternative drivers that could explain these morphological changes (e.g.exposure to parasites and pathogens) and assess the consequences of these changes on foraging performance or pollination efficiency (Gérard et al., 2020).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Impact of latitude and crop type on honeybee (A) and bumblebee (B) wing size.Red dots and red line represent the workers sampled in apple orchard, blue dots and blue line represent the workers sampled on oilseed rape.

Table 1
Total dataset used for morphological analyses.It contains 7238 analysed wings sampled across eight countries, within two bee species (A.mellifera and B. terrestris) and two sampling sessions (T0 and T1).