Ecosystem services directly affect human welfare, health and economy (Millennium Ecosystem Assessment 2005). Amongst them, pollination by animals, especially by insects, results in both ecological and economic outcomes that are beneficial for humans and is one of the most crucial ecosystem services (Potts et al. 2016; Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services 2019). This service enhances, amongst others, the maintenance of sexual reproduction of numerous wild plant species, the functioning of associated ecosystems, as well as the global crop yields (Klein et al. 2007; Bretagnolle and Gaba 2015 ; Potts et al. 2016). Bees are the most diverse, specialised and effective group of pollinators, with approximately 20,000 species worldwide and 2,000 species in Europe (Ascher and Pickering 2014; Nieto et al. 2014; Potts et al. 2016). The European honeybee (Apis mellifera L.) is the most widely managed pollinator species for crop pollination (Klein et al. 2007; Garibaldi et al. 2017; Aizen et al. 2020), as well as for the products provided by colonies. What is less well known is that more than 50 bee species are now domesticated and used for pollination of a range of crops, including the widespread use of Bombus terrestris L. for pollination of horticultural crops, such as courgettes and tomatoes (Velthuis and van Doorn 2006; Knapp et al. 2018) and that both managed and wild bee species contribute significantly to crop pollination globally (Garibaldi et al. 2013).

Several anthropogenic factors affect pollinator abundance and diversity and directly threaten the pollination services they provide (Potts et al. 2010; Vanbergen and the Insect Pollinator Initiative 2013; Klein et al. 2018; Dicks et al. 2020). Potts et al. (2016) identified five main drivers of pollinator decline: (i) landscape changes and the decrease in floral resources (quality and quantity), (ii) climate change, (iii) misuse/overuse of pesticides, (iv) worldwide trade of managed species and their associated pathogen spread and (v) invasive exotic species. The impacts of these factors on bee population trends (Duchenne et al. 2020) and bee physiology (Roger et al. 2017; vanEngelsdorp et al. 2017; Martinet et al. 2020) are increasingly well established. However, the impacts of these anthropogenic stressors on bee morphology and immunity have received less attention, especially under field conditions.

Yet, both morphology and immunity are, naturally, of great importance for bees. First, size and shape can be particularly crucial for insect fitness. Larger body size is often associated with larger foraging ranges, higher survival rate, higher fecundity and reproductive success (Greenleaf et al. 2007; Kingsolver and Huey 2008; Beukeboom 2018). At a continental scale, the body size of bee assemblages naturally increases with latitude (Gérard et al. 2018), but this trend seems less pronounced over the last century, where both body size increase (Gérard et al. 2019; Gérard et al. 2020) and decrease (Nooten and Rehan 2019) have been observed in relation to latitude, depending on the bee species. In addition to body size, the wing morphology of bees is also functionally essential for their flight performance (Wootton 1992; Wakeling and Ellington 1997) and, therefore, also their foraging and dispersal abilities (Bots et al. 2009; Johansson et al. 2009). Overall, two main buffering mechanisms are known to prevent the phenotype from undergoing deleterious changes: (i) canalisation, measured as the phenotypic consistency at the inter-individual level and (ii) developmental stability, assessed at the intra-individual level by measuring fluctuating asymmetry (FA). FA is the deviation from perfect bilateral symmetry due to random phenotypic deviations between right and left sides of an organism during its development, i.e. developmental stability decreases with increasing FA (Waddington 1957; Palmer and Strobeck 1986; Debat and David 2001). Both of these buffering mechanisms can be impacted by biotic and abiotic factors, resulting in phenotypic modifications (Hoffmann et al. 2002;Hoffmann et al. 2005), particularly under stressful conditions, such as elevated temperatures during development (Beasley et al. 2013).

In addition to their potential impacts on morphology, environmental stressors can also affect the immunological capacity of bees (Brandt et al. 2016, Di Prisco et al. 2013;Di Prisco et al. 2013) . The fat body of bees is the main tissue involved in the synthesis of immunoproteins, making this tissue a suitable proxy for bee immunocompetence (Hetru et al. 1998). Experimental and field studies showed that fat body mass is dependent on both physiological and environmental conditions and is a good indicator of immunocompetence in individual bees, defined as the capacity to mount an immune response (Keeley 1985; Korner and Schmid-Hempel 2005; Alaux et al. 2010; Vesterlund et al. 2014). In addition to this important role in the immune system, the fat body is also the central storage tissue for nutrient and energy reserves, i.e. carbohydrates, lipids and proteins (Arrese and Soulages 2010). While the impact of potential threats on the phenotype and the immunocompetence of bees are well studied under laboratory conditions (Gérard et al. 2018b for phenotype; Alaux et al. 2010; Vanderplanck et al. 2016 for fat body), there are, as yet, too few field studies (e.g. Dellicour et al. 2017) to assess the suitability of bee phenotype and fat body as indicators of impacts of global change under real-world, landscape-level conditions.

In this article, we compile both wing morphology and fat body data for two major pollinator species of European crops, the honeybee (A. mellifera) and the buff-tailed bumblebee (B. terrestris), following field-level exposure to stressors in two major entomophilous crops (perennial apple trees and annual oilseed rape plants) at 128 sampling sites in eight European countries and two sampling occasions: a first batch of specimens collected before crop bloom and a second batch after the crop bloom had ceased and the larval and imago stages of the bees had potentially been exposed to pesticides used in the crop and surrounding landscapes. The dataset also compiles the environmental factors associated with each sampling site, namely landscape metrics and local pesticide use. This dataset offers the opportunity to test whether phenotypic variability and immunocompetence of bees can be affected by a range of real-world, landscape-level environmental drivers in a context of global change.

Overview of the study

This study was carried out in the framework of the Horizon 2020 project PoshBee (http://poshbee.eu), which aims to support healthy bee populations, sustainable beekeeping and pollination in Europe. The goal of Work Package 1 (WP1) is described as “developing a site network for assessing exposure of bees to chemical, nutritional, and pathogen stressors” and led by Trinity College Dublin (Ireland). WP1 also includes 30 additional collaborators across 14 European countries (http://poshbee.eu/partners). The goal of PoshBee Work Package 2 (WP2) is described as “measuring chemical exposure, pathogens and aspects of nutrition in honey bees, bumble bees and solitary bees” and is led by the Agence Nationale de la Sécurité Sanitaire de l’alimentation, de l’environnement et du travail (ANSES, France). WP2 also includes seven additional beneficiaries across five countries.

Eight apple orchards sites and eight winter-sown oilseed rape sites were selected according to a gradient of land-use intensity in each of eight countries chosen to represent four major European biogeographical areas: Boreal (Sweden and Estonia), Atlantic (Ireland and United Kingdom), Continental (Germany and Switzerland) and Mediterranean (Spain and Italy), making a total of 128 different sampling sites. Three honeybee hives (A. mellifera) and three B. terrestris colonies were placed in each site and standardised following internal PoshBee protocols (Hodge and Stout 2019). Two sampling sessions were conducted at each of these sites during the local crop blooming period: a first sampling session (T0) when the hives/nests were installed on the sites (between late March and mid-May 2019, depending on the country) and a second sampling session (T1) after bloom of the focal crop and potential pesticide exposure (between mid-May and late July 2019, depending on the country).

At each site, several specimens of B. terrestris and A. mellifera were collected during each sampling session to perform the morphological and fat body analyses. For the wing morphological analyses, eight individuals of each species were collected per site and per sampling session, for a total of 5096 specimens (1024 individuals per species per sampling session). Several individuals had to be excluded from the analyses due to wing damage or because some teams were not able to collect enough individual bees (Specimens from T1 Ireland have been unfortunately lost somewhere during the delivery between Ireland and Belgium and bumblebees from T1 Estonia seem to have been sent to another university, without being able to find them). The number of individuals used to compute the morphological analyses is summarised in Table 1.

For the analyses of the fat body mass, five specimens of each species were analysed per site for the second sampling session (T1), for a total of 990 specimens. The first sampling session was not considered for these analyses since the T0 specimens were also used for other analyses that involved destructive methods. As with the morphological dataset, some teams were not able to collect enough specimens and some specimens were excluded from the analyses because they were compromised (see Table 2 for a summary of the number of specimens available).

The sampled specimens enable the evaluation and comparison of the phenotypic variability and fat body content in bees that were exposed to agrochemicals along a gradient of pesticide use intensity, including exposure during the larval stage, with similar metrics from non-exposed bees.

Morphometric measurements

Morphometric measurements were conducted using an Olympus SZH10 microscope coupled with a Nikon D200 camera to photograph each bee wing. After uploading pictures in the tpsUTIL 1.69 software (Rohlf 2013a), we digitised left and right forewings of each specimen with a set of 18 two-dimensional landmarks (Fig. 1; tps DIG 2.27, Rohlf 2013b). Then, using the function readland.tps from the geomorph package, each landmark coordinate is multiplied by its scale factor provided for each specimen in the TPS file created by the software tpsUTIL (Adams and Otárola-Castillo 2013; R Core Team 2017; Suppl. material 1). We then used the Generalized Procrustes Analysis superimposition method that remove all the non-shape components, i.e. by translating specimens to the origin, scaling and rotating each landmark configuration to minimise the distance between each corresponding landmark of each landmark configuration. To do this, we used the “gpagen” function of geomorph package (Adams and Otárola-Castillo 2013; R Core Team 2017). We used centroid size (CS) to estimate the wing size, which is calculated as the square root of the sum of squared distance between all landmarks and the wing centroid (e.g. Gérard et al. 2015). To test for measurement error, a sub-dataset of 128 wings was digitized twice by the same experimenter (MG). The statistical procedure was undertaken on R version 3.6.4.

Fat body

The abdominal fat body content of specimens was measured following Ellers 1996. Isolated abdomens were weighed after drying at 70°C for 3 days. Dried abdomens were then immersed in 2 ml of diethyl ether for 24 h to extract fat, rinsed twice and weighed again after drying at 70°C for 7 days. The fat body content was defined as the abdominal weight loss during this process, standardised by abdomen weight before extraction to avoid biases linked to specimen size (Suppl. material 2).

Environmental predictors

In addition to the type of crop, we collected information regarding different types of environmental predictors:

• Landscape structure: Within a 1-km radius centred at the location of sentinel colonies, all land cover features were manually identified using high-resolution images provided by World Imagery (ESRI) and digitised using Geographical Information Systems Software (ArcGIS Pro, 2.4.1, ESRI). Identified land cover features were then classified into ten final categories: surface running waters, waterbodies, wetlands, grasslands, woodlands and heathlands, bare areas, orchards, crops, roads and urban areas. Landscape structure was quantified by calculating five independent metrics of landscape composition and configuration using the software FRAGSTATS 3.3 (McGarigal et al. 2012). First, the proportion (%) of all ten land cover features was calculated for each landscape. A landscape intensity gradient (LIG) was then defined using the sum of proportion of crops and orchards per landscape as a proxy. As measures of landscape composition, the Shannon’s Diversity Index (SHDI) and the total area (CA) of both grasslands and urban areas was calculated. As a measure of landscape configuration, edge densities of semi-natural habitats (ED) were measured by dividing the edge length of semi-natural habitats by the total area of the corresponding habitat map.

• Pesticide use: Each site experienced different levels of pesticide application in focal crop field. We used the sum of pesticide applications (measured as the sum of kg/ha and l/ha of plant protection products applied) as a proxy for the intensity of pesticide application. Pesticide data (all organic and synthetic herbicides, fungicides and insecticides applied to the field) for each field was acquired through directly questioning farmers who own or lease the sites from which bees were sampled. For each pesticide and each application, farmers provided application rates (in l/ha or kg/ha depending on the pesticide) by date (ranging from 1-5 applications per pesticide per site). Only applications between October 2018 and June 2019 (period preceding specimen collection) were considered. Pesticide use intensity is the sum of all applications of all pesticides over that period. This is only a general proxy and does not account for the relative toxicity of different pesticides, the volume of active ingredients or the impacts of applications at different times during bee life cycles. We also included the active(s) ingredient(s) (AI) contained in the pesticides used. Not all farmers provided all the requested information: we obtained this information for 83 out of the 128 sites.

Bee wing morphology

An excel table with 7238 rows (without column headings) and 80 columns. Each row represents a bee wing.

Column headings: Specimen ID, Side, Individual, Session, Species, Replicate, Country, Crop, Latitude, Landscape intensity gradient, SNH edge density, CA Grassland, CA Urban, Shannon Diversity Index, Pesticide use, centroid size, x and y coordinates of the 18 landmarks on the bee wings, AI1 – AI28.

Geographical coverage: Eight European countries (Estonia, Germany, Ireland, Italy, Spain, Sweden, Switzerland and United Kingdom).

Spatial resolution of the landscape structure: 1 km radius around the centroid of the site.

Input data:

Data: Morphometric measurements (i.e. centroid size as well as 18 landmarks to quantify wing shape

· Specimen ID: a code containing the type of crop, the number of the site, the individual code and the side

· Side: right (R) or left (L) forewing

· Individual: the individual code of a specimen

· Session: the session to which the specimen has been collected (i.e. T0 or T1)

· Species: the species to which the specimen belongs (i.e. A. mellifera or B. terrestris)

· Replicate: first or second session of landmark digitalisation (1 or 2)

· Country: the country to which the specimen belongs (i.e. Estonia, Germany, Ireland, Italy, Spain, Sweden, Switzerland, UK)

· Centroid size (continuous quantitative data; mm): square root of the sum of squared distance between all landmarks and the wing centroid

· X and y coordinates of the 18 landmarks on the bee wings in a Cartesian coordinate system of origin (0,0) (continuous quantitative data; mm): 18 columns of the x coordinate of each landmark and 18 columns of the y coordinate of each landmark

Environmental conditions:

· Different variables characterise the landscape structure (continuous quantitative data): Shannon’s Diversity Index (“no units”), total area of grasslands (CA Grassland; hectares), total area of urban areas (CA Urban; hectares), landscape intensity gradient (“no units”), semi-natural habitat edge density (SNH edge density; metres per hectare)

· Pesticide use (continuous quantitative data; l/ha)

· AI (1-28): names of the Active Ingredients contained in the pesticides used

· Type of crop (qualitative data with 2 levels)

· Latitude (continuous quantitative data, decimal degrees)

Object name

TableS1_Morpho.xlsx

Creation date

June 2020

Dataset creator

Maxence Gérard

Dataset contributors

See list of co-authors

Repository location

This paper (Suppl. Material S1)

Bee fat body

An excel table with 990 rows (without columns heading) and 41 columns. Each row represents an individual.

Columns heading: Individual, Field label, Species, Country, Type of Crop, Latitude, Landscape intensity gradient, SNH edge density, CA Grassland, CA Urban, Shannon Diversity Index, Pesticide use, Fat Body.

Geographical coverage: Seven European countries (Estonia, Germany, Italy, Spain, Sweden, Switzerland and United Kingdom).

Spatial resolution of the landscape structure: 1 km radius around the centroid of the site.

Input data:

Data:

· Individual: the individual code of a specimen

· Field label: a code containing the country, the type of crop, the species and the number of the site

· Species: the species to which the specimen belongs (i.e. A. mellifera or B. terrestris)

· Country: the country to which the specimen belongs (i.e. Estonia, Germany, Ireland, Italy, Spain, Sweden, Switzerland, UK)

· Fat body: Fat body content (quantitative data; %)

Environmental predictors:

· Different variables characterise the landscape structure (continuous quantitative data): Shannon’s Diversity Index (“no units”), total area of grasslands (CA Grassland; hectares), total area of urban areas (CA Urban; hectares), landscape intensity gradient (“no units”), semi-natural habitat edge density (SNH edge density; metres per hectare)

· Pesticide use (continuous quantitative data; l/ha)

· AI (1-28): names of the Active Ingredients contained in the pesticides used

· Type of crop (qualitative data with 2 levels)

· Latitude (continuous quantitative data, decimal degrees)

Object name

TableS2_FatBody.xlsx

Creation date

November 2020

Dataset creators

Victor Lefebvre and Maryse Vanderplanck

Dataset contributors

See list of co-authors

Repository location

This paper (Suppl. Material S2)

Re-Use potential

Any re-use of these data must cite this source. The authors may be contacted in case of doubts with the use and the interpretation of the data.

The views expressed in this article are those of the authors and do not necessarily reflect an official position of the European Union.

This research has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement No 773921 for the POSHBEE project. We acknowledge all the farmers who participated in the project, allowing us to perform our study in their crops and kindly answering our questionnaire.

The authors declare that they have no competing interests.