# Building composite indices in the age of big data – Application to honey bee exposure to infectious and parasitic agents

Minh Huyen Ton Nu Nguyen, Stéphanie Bougeard, Aurélie Babin, Eric

Dubois, Christine Druesne, Marie-Pierre Rivière, Marion Laurent,
Marie-Pierre Chauzat

► To cite this version:

Minh Huyen Ton Nu Nguyen, Stéphanie Bougeard, Aurélie Babin, Eric Dubois, Christine Druesne, et al. Building composite indices in the age of big data – Application to honey bee exposure to infectious and parasitic agents. Heliyon, 2023, 9 (4), pp.e15244. 10.1016/j.heliyon.2023.e15244. anses-04068240

HAL Id: anses-04068240

https://anses.hal.science/anses-04068240v1

Submitted on 9 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Building composite indices in the age of big data - Application to honey bee exposure to infectious and parasitic agents

Huyen Ton Nu Nguyen¹, M., ²Bougeard S., ³Babin A., ³Dubois E., ⁴Druesne C., ³Rivière M.P., ³Laurent M., ¹,³Chauzat, M.P.

## Affiliations.

¹Paris-Est University, ANSES, Laboratory for Animal Health, Maisons-Alfort, France.

²ANSES, Ploufragan-Plouzané-Niort Laboratory, Epidemiology and welfare of pork, France.

³ANSES, Sophia Antipolis laboratory, Unit of Honey bee Pathology, France.

⁴ANSES, Research funding & scientific watch department, Maisons-Alfort, France.

## Corresponding author.
marie-pierre.chauzat@anses.fr. ANSES. Animal Health Laboratory (LSAn). 14, rue Pierre et Marie Curie F-94701 Maisons-Alfort Cedex. France. https://orcid.org/0000-0001-8842-741X

## Abstract.

Pollinator insects play a crucial role in maintaining biodiversity and agricultural production worldwide. Yet they are subject to various infectious and parasitic agents (IPAs). To better assess their exposure to IPAs, discriminative and quantitative molecular methods have been developed. These tools produce large datasets that need to be summarized so as to be interpreted. In this paper, we described the calculation of three types of composite indices (numerical, ordinal, nominal) to characterize the honey bee exposure to IPAs in 128 European sites. Our summarizing methods are based on component-based factorial analyses. The indices summarized the dataset of eight IPAs quantified at two sampling times, into synthetic values providing different yet complementary information. Because our dataset included two sampling times, we used Multiple Factor Analysis (MFA) to synthesize the information. More precisely, the numerical and ordinal indices were generated from the first component of MFA, whereas the nominal index used the first main components of MFA combined with a clustering analysis (Hierarchical Clustering on components). The numerical index was easy to calculate and to be used in further statistical analyses. However, it contained only about 20% of the original information. Containing the same amount of original information, the ordinal index was much easier to interpret. These two indices summarized information in a unidimensional manner. Instead, the nominal index summarized information in a multidimensional manner, which retained much more information (94%). In the practical example, the three indices showed an antagonistic relationship between *N. ceranae* and DWV-B. These indices represented a toolbox where scientists could pick one composite index according to the aim pursued. Indices could be used in further statistical analyses but could also be used by policy makers and public instances to characterize a given sanitary situation at a site level for instance.
**Keywords.** Summarizing big datasets. Stress exposure. Honey bee health. Factor analysis. Synthetic indices.

**Highlights.**

- The dataset of honey bee exposure to IPAs was summed up using three type of composite indices

- Numerical, ordinal and nominal indices relied on factorial statistical methods

- Pros and cons of the three composite indices were discussed using applications

- The three indices showed an antagonistic relationship between *N. ceranae* and DWV-B

**Graphical abstract**

<figure><img src="image_2.png" /></figure>

**1. Introduction**

Pollinating insects play a crucial role in maintaining biodiversity and agricultural production worldwide. Bees are the most important group of pollinators, visiting more than 90% of the leading 107 global crop types ¹. Honey bees (*Apis mellifera*) are the most used pollinators in agriculture. They are essential to 35% of the world's agricultural production, animal pollination services representing an annual economic contribution of up to 200 to 500 billion US dollars ². Under normal conditions, 10% is the acceptable level of colony winter losses for beekeepers ³. However, the number of honey bee colonies in most European countries has been declining since 1965, and different causes for these losses have been proposed ⁴. The intricate co-exposure to various stressors, such as infectious and parasitic agents (IPAs), chemicals released in the environment (i.e. insecticides, herbicides, fungicides), the depletion and poor quality of food sources, inappropriate beekeeping practices, and meteorological or climatic causes, as well as the interaction between them, weakens bee health and generates a worrying colony over-mortality. Consequently, studying
73 the stress effects on bee health is becoming a crucial economic and ecological issue.
74 The advent of new technologies made possible for biologists the collection of large
75 datasets integrating information from multiple sources. Therefore, researchers need
76 methods that summarize the information contained in these large datasets in order to
77 perform statistical analyses and get meaningful insights from data.

78 Only a few methods are proposed in the literature to build indices for animal health.
79 In most animal breeding, many models developed to characterize population health
80 use methods that are either arbitrary ⁵ or not easily transferable to another animal
81 species. Indeed, these models are tailored to or heavily dependent on the species
82 characteristics ⁶⁻⁹. Especially in the bee world, in which most studies focused on
83 honey bee populations, no previous work has reported the development of an index
84 summarizing their exposure to multiple IPAs (the PRISMA bibliographic survey is
85 detailed in Appendix A, supplementary material). Only two publications described a
86 synthetic approach to explore the relationship between pesticide exposure and
87 pathogen prevalence in honey bees ¹⁰,¹¹. The study of large datasets is possible
88 using factor analyses, also known as component-based analyses ¹²,¹³. These
89 analyses are usually used as they provide synthetic graphical displays illustrating
90 relationships between variables and between observations. Principal Component
91 Analysis – PCA ¹²,¹⁴ or Multiple Correspondence Analysis – MCA ¹⁵,¹⁶ are some of
92 the most popular methods of this family. However, factor analyses, and more
93 precisely components sought from factor analyses, could also be used to seek
94 indices.

95 In this paper, we focused our statistical analyses on IPAs. Our goal was not to
96 include other factors (climate, landscape) in the indices. The proposed synthetic
97 indices are not designed to describe multiple stressors affecting bees coming from a
98 large dataset. We propose to apply a standard and transferable procedure based on
99 factor analysis and adapted to the configuration of biological data, to seek an index,
100 regardless of the nature and structure of the data. Depending on the needs of the
101 biologists and/or the goal of communication, this index can be either numerical,
102 ordinal (=ordered qualifier) or nominal (=unordered qualifier). These three types of
103 indices were calculated on data of honey bee exposure to IPAs produced by the
104 European project PoshBee on 128 sites (https://poshbee.eu/). We aimed at
105 summarizing the several variables (i.e. eight IPAs quantified at two sampling times)
106 into a single index that best described and discriminated the sites.
107

107

108

## 2. Materials and methods

109

### 2.1. Data

110 Data were produced in the frame of a Europe-wide project (PoshBee) which aimed at
111 assessing, monitoring and mitigating stressors on bee health. PoshBee involved 14
112 European Union member countries over a five-year period (2018-2023). By
113 "developing a site network for assessing exposure of bees to chemical, nutritional,
114 and pathogen stressors" (Work Package 1), and by "measuring chemical exposure,
115 pathogens and aspects of nutrition in honey bees, bumble bees and solitary bees"
116 (Work Package 2), PoshBee provided several datasets with unique, robust and
117 complete information ¹⁷. Data were collected in 128 sites, with two focal crops, across
118 eight European countries: Estonia (site numbers beginning with EST), Germany
119 (GER), Ireland (IRL), Italy (ITA), Spain (ESP), Sweden (SWE), Switzerland (CHE),

110

111

112

113

114

115

116

117

118

119
and United Kingdom (GBR). This means that eight sites per country for two crops were selected: 8 sites x 8 countries x 2 crops, leading to a total of 128 sites. At each site, three honey bee colonies were installed following the PoshBee protocols and over 50 measurements and samples related to the nutritional, toxicological, pathogenic and landscape components of the bees' environment were collected ¹⁸. In the following list, we provide examples of collected information: data on the abundance and diversity of field margin floral resources (floral survey); survey of wild pollinators as an indicator of general 'pollinator community health (pollinator survey); total and per hectare crop yields and inputs of agrochemicals (farmer survey) ¹⁸. Adult honey bees were sampled in each colony and nest two times in 2019 for the detection and quantification of various IPAs at two sampling times: before bees were deployed on the sites (T0) and after bloom of the focal crop (T1). Bees from the three colonies and nests were pooled respectively. Ahead of the real-time qPCR analysis, samples were ground for nucleic acid extraction. Nucleic acids (RNA and DNA) were purified on spin columns (Machery-Nagel) with the high-throughput automated extraction of nucleic acids (TECAN pipetting robot). For both viruses and bacteria, standard qPCR protocols validated by the European reference laboratory for diagnosis (www.eurl-bee.anses.fr) were optimised with the reagents and the parameters for high-throughput qPCR ¹⁹. For the three microsporidia species, a specific method was developed and validated using relevant molecular targets for an accurate quantification by real-time qPCR ²⁰. The loads of these IPAs were expressed in number of copies of the target IPA genome transformed into log₁₀ per honey bee (log₁₀ ge/bee).

The dataset contained 2,805 observations, corresponding to the loads of 11 IPAs (initially 11 and subsequently reduced to 8 IPAs, see below) quantified at two sampling times (T0 and T1) on 128 sites: six viruses (Acute bee paralysis virus - ABPV, Black queen cell virus - BQCV, Chronic bee paralysis virus - CBPV, Deformed wing virus type A - DWV-A, Deformed wing virus type B - DWV-B, Sacbrood bee virus - SBV), the bacterial agents of the American foulbrood (*Paenibacillus larvae*) and the European foulbrood (*Melisococcus plutonius*) and three *Nosema* microsporidia (*Nosema apis*, *N. ceranae* and *N. bombi*). Qualitative variables which described the country and the focal crop were also included in the dataset. The theoretical number of analyses was 2,816 (11 IPAs sampled in 128 sites at 2 times). After quality control, some samples were discarded (0.01%).

Amongst the 11 IPAs, *P. larvae*, *M. plutonius* and *N. bombi* were not detected in more than 90% of the samples at the two sampling times. Consequently, these three IPAs were not included in the development of the indices due to their very low presence (Table 1).

## 2.2. Statistical methods

It should be reminded that although the statistical methods presented here were not novel, this was the first time that these methods were used to characterize the exposure of bees to IPAs. This synthesis was a first step to interpret the exposure data. Subsequent steps may involve the use of structural equation modeling (SEM) for example.

The index was defined as a value assigned to a site, that summarises all the variables measured for that given site. Because components of factorial analyses are the best summary of several variables into a single one, these statistical methods
were applied. Indeed, components are built as linear combinations of the variables obtained from a Singular Value Decomposition of the overall dataset ²¹. Therefore, the first component explains the maximum amount of data variability.

The statistical procedure was organised in two steps. The first step consisted of reducing the dataset into a single (or into its main) component(s). The second step involved the calculation of the three indices from the first (or the main) component(s). Depending on the biologist's aim, three types of indices are proposed (numerical, ordinal or nominal). The first component was directly used to seek the numerical index, quantiles of the first component were used to seek the ordinal index, whereas several main components were used to seek the nominal index.

**Step 1: Seek the first or main component(s) that best summarise the data.** To summarise the large dataset into a single index, a factor analysis was performed and its first component was extracted.

For our data, the aim is to rank or classify into categories the 128 sites according to the bee exposure to the various IPAs. Since the variables are grouped into two sampling times (i.e. T0 and T1), a Multiple Factor Analysis (MFA) was applied ¹⁶.

**Step 2: Process the index from the component(s).** Three types of indices could be calculated and selected by the biologist depending on the aim pursued:

(i) a numerical index directly obtained from the first MFA component. It was the best summary of the whole dataset into a single dimension, the one that explained the greatest variability of the data;

(ii) an ordinal index from the quantiles of the first MFA component. This index was a transformation of the numerical index into ordinal levels using quantiles;

(iii) a nominal index produced from a clustering analysis (Hierarchical Clustering on Principal Components - HCPC ¹⁴) applied to the first main MFA components. This method gave the best multidimensional summary of the whole dataset. The choice of the final number of clusters to be retained was discussed among statistical, biological and bee health experts.

The MFA method was applied to the eight IPAs measured at two different times. Country and focal crop of each site were considered as supplementary variables in the analysis. They were not used to process the components but only for additional interpretation.

The factor analyses were implemented in the R software version 4.1.2 (R Core Team, 2017) from the FactoMineR package ²².

# 3. Results: application to honey bee exposure to IPAs

## 3.1. Calculation and interpretation of the numerical index

The coordinates of the 128 sites on the first component of the MFA (horizontal axis) corresponded to the numerical index values (Figure 1a and Appendix C). The index ranged from -3.58 for the 82th site to 2.32 for the 97th site, with a mean value of 0 (sd = 1.31). It explained 20.1% of the total data variability (Appendix B). The IPAs DWVB, SBV and *N. ceranae* were those which mostly contribute to build this
component/index (Figure 1b). In other words, the degree of exposure to IPAs that differentiates the 128 sampled sites was mostly characterised by these three IPAs. Moreover, the first component correlated congruently with the IPAs load on both sampling time, but not with the focal crop (Figure 1). Therefore, no interpretation concerning the focal crop would be made, and unless mentioned otherwise, all interpretations applied to either modality of the group variable (T0 and T1).

To extract the numerical index values, the site coordinates were projected onto the first component, as illustrated in red for 5 out of 128 sites. Using the numerical coordinates, sites were ranked according to their exposure to the most contributing IPAs. Insights into the characteristics of each site were obtained based on a juxtaposition of the representation of sites (Figure 1a) and variables (Figure 1b): for example, sampled sites with negative coordinates (so-called "low index values") were characterised by a higher load of DWV-B and lower load of SBV and *N. ceranae* compared to the sites that have an index value around 0; whereas sites with positive coordinates (so-called "high index values") had a lower load of DWV-B and higher loads of SBV and *N. ceranae*.

This index can also be interpreted with the supplementary variable 'country' (Figure 1a). The low-index-value sites were particularly common in Ireland, Spain and Sweden whereas the high-index-value sites were commonly found in Estonia, Italy, Switzerland and the United Kingdom. Most sites with a value around 0 were from Germany.

## 3.2. Calculation and interpretation of the ordinal index

In order to simplify the index interpretation, the latter numerical index could be transformed into an ordinal one by using quartiles for instance. The first quartile (Q1) corresponded to the limit below which 25% of the data is set. The second quartile (Q2) was the median (i.e., 50% of the data lies below this value). The third quartile (Q3) was the limit below which 75% of the data is set. By dividing the 128 sites into four equal quartiles (32 sites each), according to the relationships between each site coordinates and the cut-off value of each quartile, an ordinal index ranging from one to four was obtained. The first 'site' quartile clustered sites whose coordinates on the first component were strictly inferior to Q1=-1.10. The second quartile clustered sites whose coordinates were comprised between Q1=-1.10 and Q2=0.21. Sites with coordinates between Q2=0.21 and Q3=0.96 belonged to the third 'site' quartile and the last 'site' quartile grouped the remaining sites (Appendix C).

The meaning behind the score and each configuration of health status was the same as that of the numerical index. The higher the score, the lower the loads of DWV-B and the higher SBV and *N. ceranae* loads (Table 2 for the main features illustrated in Figure 2).

## 3.3. Calculation and interpretation of the nominal index

The first step consisted in selecting the number of components of the MFA to be included in the clustering. This step led to (i) take into account the block structure of the data (T0 block and T1 block), and (ii) use the relevant part of the information (and avoid noise) while selecting components. According to the decrease of inertia (Appendix B), 12 components (out of 16, the total number of IPAs for the two times
259
260

T0 and T1), for a cumulative 94% of the data inertia, were a relevant selection of the information contained in the data.

261
262
263
264

The second step consisted in selecting the best number of clusters. To do this, examination of criteria such as maximum loss of intra-class inertia and Ward's criterion suggested two-, three- or four-class clustering. According to the biologists' choice in relation with interpretation, a four-class nominal index was interpreted.

265
266
267
268
269

The nominal index was characterised mostly by the load of all the IPAs of interest at both sampling times. To emphasise the non-ordinal characteristics of the clusters, they were renamed with a letter instead of the default number (Figure 3). Table 3 (and Appendix C) summarised the characteristics of the clusters as well as the distribution per country of the sampling sites in each cluster.

270
271
272
273
274
275
276
277
278
279
280
281
282

The clusters A, B, C, and D included respectively 19, 34, 57 and 18 sites. Both clusters A and B had significantly lower loads of *N. ceranae* compared to the average load of all sites. Cluster A was mostly characterised by significantly lower loads of BQCV, CBPV and SBV than the average load. Cluster B had significantly higher loads of ABPV, BQCV, DWV-A and DWV-B at T1, as well as a significantly lower loads of *N. apis* at T1. Cluster C represented the majority of the sampled sites. As opposed to clusters A and B, it had significantly higher loads of ABPV, CBPV, SBV and *N. ceranae* at T0 and of CBPV at T1 than the average. Cluster C was also characterised by significantly lower loads of DWV-A and DWV-B at T0 and *N. apis* at both sampling times. The last cluster (D), was characterised by significantly higher loads of BQCV and *N. apis* at both sampling times and *N. ceranae* at T0, as well as significantly lower loads of DWV-B at T1 than the average. In summary Cluster C and D were characterised by high loads of *N. ceranae* as opposed to clusters A and B.

283
284
285
286
287
288

Interestingly, the country of origin was significantly related to the clustering as some countries were over-represented in some clusters (Table 3 for the main features illustrated in Figure 4). Cluster A grouped mostly Irish sites, while cluster B was dominated by Spanish and Swedish sites. The majority of sampled sites were found in cluster C with sites from Germany, Italy and Switzerland. The majority of British and Estonian sites were identified in cluster D.

289
290
291
292
293
294
295
296
297
298
299
300

Overall, we found the same antagonistic relationship between *N. ceranae* and DWV-B across these four clusters as observed for the numerical index (Figure 4). However, the information was completed with more nuances found for the six remaining IPAs as 94% of the information was contained in this index (instead of 20% for the numerical index). Among low-index-value sites, Ireland differed from Spain and Sweden in having an overall lower load of BQCV while the prevalence was rather high in those sites (cluster B). Among high-index-value sites, Estonia and United Kingdom differed from the other sites with a higher prevalence of *N. apis*. In addition, nuances related to the sampling times were better highlighted with the nominal index than with the numerical or the ordinal ones. A summary table (Appendix C) gave a full disclosure of the three types of indices for each PoshBee site.

301

302

# 4. Discussion and conclusion

303
304
305

In this paper, we described the calculation of three types of synthetic indices to characterise the honey bee exposure to IPAs. Depending on (i) the variable structure (i.e., variables organised into a single block of variables or into several ones) and (ii)
on the variable's nature (i.e., quantitative or qualitative values; or mixture of
quantitative and qualitative values), different factor analyses could be applied. More
precisely, PCA<sup>23</sup> was applied to numerical variables whereas MCA<sup>15,24</sup> was applied
to nominal ones. In the case of a mixture of numerical and nominal variables, the
FAMD method could be used<sup>13</sup>. When variables were organised into several blocks
(e.g., variables measured at different times or variables organised into meaningful
categories of characteristics such as IPAs, pesticides and nutritional resources), a
MFA could be used<sup>25</sup>.

These three indices (numerical, ordinal and nominal) summarised the rich information
of eight IPAs quantified at two time points into synthetic values. They provided
different yet complementary information. Similar to a toolbox, they could be chosen
according to the aim pursued. The **numerical index** reflecting the coordinates on the
first component was easy to calculate and implement in further statistical analyses
since each site had a distinct value of exposure on a continuum. It could be easily
included in all standard statistical analyses (i.e., regression). However, only a limited
amount of the original information (when rich and multidimensional) was contained in
this index. In our example, the 128 sampled sites were mostly characterised by the
loads of DWV-B, SBV and *N. ceranae*. If the aim was to include the complex data on
exposure to IPAs into further statistical analysis, the numerical index should be
chosen bearing in mind that the information contained in this index was highly
simplified when data are rich. In order to ease the interpretation of the numerical
index, an **ordinal index** using quartiles of the coordinates on the first component was
created. This ordinal index was much easier to interpret than the numerical one.
However, it had the same limits in the case of multidimensional data: it summarised
only a small amount of information contained in the data (20% of the information in
our study). DWV-B, SBV and *N. ceranae* were also the discriminative IPAs that
characterised low and high index values. In our example and for this index, each site
belonged to a 4-level ordered scale of exposure. The amount of information
summarised in the **nominal index** was the highest among the three types of indices
(94%) but the index was less practical to use in subsequent statistical analyses as
only a reduced number of statistical analyses are currently devoted to nominal data.
If the biologist needed to interpret data with deep insight, the nominal index would be
the most relevant tool. In our example, the dataset resulted in four distinct categories
of exposure configuration, including the same antagonistic relationship between *N. ceranae* and DWV-B as observed in the numerical index (Appendix D). However, the
information was completed with more nuances found for the remaining IPAs. These
results support the previous hypothesis that variation in tolerance might be much
more prevalent in systems with less pathogenic parasites (our study: BQCV, DWV,
SBV or *N. ceranae*) than in those with more pathogenic parasites (our study: ABPV,
CBPV or *N. apis*), and where the costs of tolerance are expressed not only in
resistance<sup>26</sup>. In other words, less pathogenic IPAs were found more often than IPAs
more lethal for their hosts. Therefore loads might be more relevant for assessing bee
health. Our results were consistent also with hypotheses of an antagonistic
relationship between *N. ceranae* and DWV in the absence of *Varroa destructor* mite
<sup>27</sup>. Doublet et al. (2015) experimentally explored this hypothesis with oral infections of
honey bees and found that the prior establishment of *N. ceranae* in midgut negatively
affected the development of DWV; but the competition between both pathogens was
asymmetric and the prior infection by DWV did not prevent the *N. ceranae* infection
<sup>28</sup>. Contrary to country, focal crops were not significantly linked to IPA exposure in
honey bees. Nevertheless, time of colony exposure on sites was linked to IPA
356 exposure. There was an evolution in pathogen loads between installation of colonies
357 on sites and end of the flowering bloom. One other study has also shown the impact
358 of IPA exposure on variation of morphological traits in honey bees at site level 29.
359 Interestingly, crop surrounding sites also impacted this variation. To conclude, we
360 proposed the calculation of three different types of indices to summarise a large
361 number of variables collected by biologists. This calculation could be applied to a wild
362 range of data; an example on bee data being given here. As factor analyses
363 summarise the information contained in the overall data into components (the first
364 one being the best summary into a single dimension, the first two ones into two
365 dimensions, and so on) they are relevant tools for biologists to summarise their data.
366 In addition, factor analyses are developed for various nature of variables (numerical,
367 nominal, mixture of numerical and nominal) and various structures of variables
368 (single or several blocks). Thus, these methods can be applied to all types of data to
369 summarise them into indices.

370 The factor analysis helps the biologist to make a biological interpretation of a set of
371 factors from the variable analyses. For honey bee heath, the field of applications of
372 these synthetic indices is huge and includes risk assessment data. Bees are exposed
373 to IPAs and pesticides, however these stressors do not present the same risk
374 depending on quantity and type of IPAs (virus, bacteria or else) or pesticides
375 (fungicides, insecticides or herbicides). Exposure data can be coupled with hazard
376 data (toxicological data) to express risk data. The calculation of indices described in
377 this paper can also be applied to datasets on the risks posed to honey bees by other
378 stressors (nutrition, access to habitats, others). Indices can be used for further
379 statistical analysis but can also be used by policy makers and public instances to
380 characterise a given situation.

381
382 **Acknowledgements**
383 This work would not have been possible without the tremendous work of the
384 numerous persons involved in setting up the sampling sites network and collecting,
385 sending and analyzing all the samples in the PoshBee project. The authors would
386 like to thank all of these collaborators.
387 This work has been performed in the frame work of the European Reference
388 Laboratory for Honey bee health.

389 **Funding**
390 This work was supported by the European Union's Horizon 2020 research and
391 innovation program under grant agreement No 773921.

392 **Author Contributions**
393 1 – Conceived and designed experiments: SB, MHTNN, ML, MPC
394 2 - Performed experiments: MHTNN, SB, ML, MPC
395 3 – Analysed and interpreted the data: MHTNN, AB, ED, MPR, ML, CD, SB, MPC
396 4 – Contributed reagents, materials, analysis tools or data: SB, MHTNN, AB, ED, CD,
397 ML, MPC
398 5 – Wrote the paper – MHTNN, SB, ML, AB, ED, MPR, MPC
399
400 **Competing interests**

401 The authors declare no competing interests.

402
<figure><img src="image_3.png" /><figcaption>**Appendix A.** Table A.1. Results of the bibliographic analysis following PRISMA method to look for existing indices describing honey bee exposure to stressors ³⁰.</figcaption></figure>

The scheme shows that only seven papers fulfilled the eligibility criteria.
410

**Appendix B.** Figure B.1. MFA scree plot of honey bee exposure to IPAs.

411

Percentages are proportion of the total variances explained by each dimension.

<figure><img src="image_4.png" /></figure>

412

413
414

415

416

**Appendix C. Table C.1.** The three types of indices calculated for each site to characterise honey bee exposure to IPAs.

<table><thead><tr><td>Site identification</td><td>Numerical index</td><td>Ordinal index</td><td>Nominal index</td></tr></thead><tbody><tr><td>GBR_OSR_01</td><td>0,99</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_02</td><td>1,29</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_03</td><td>2,17</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_04</td><td>2,12</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_05</td><td>1,88</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_06</td><td>1,8</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_07</td><td>1,28</td><td>4</td><td>D</td></tr><tr><td>GBR_OSR_08</td><td>1,94</td><td>4</td><td>D</td></tr><tr><td>GBR_APP_09</td><td>0</td><td>2</td><td>C</td></tr><tr><td>GBR_APP_10</td><td>0,68</td><td>3</td><td>C</td></tr><tr><td>GBR_APP_11</td><td>0,25</td><td>3</td><td>C</td></tr><tr><td>GBR_APP_12</td><td>0,73</td><td>3</td><td>C</td></tr><tr><td>GBR_APP_13</td><td>0,01</td><td>2</td><td>C</td></tr><tr><td>GBR_APP_14</td><td>-0,51</td><td>2</td><td>B</td></tr><tr><td>GBR_APP_15</td><td>0,66</td><td>3</td><td>C</td></tr><tr><td>GBR_APP_16</td><td>-0,3</td><td>2</td><td>B</td></tr><tr><td>ESP_OSR_01</td><td>-0,79</td><td>2</td><td>B</td></tr><tr><td>ESP_OSR_02</td><td>0,37</td><td>3</td><td>B</td></tr><tr><td>ESP_OSR_03</td><td>0,84</td><td>3</td><td>B</td></tr><tr><td>ESP_OSR_04</td><td>-1,18</td><td>1</td><td>B</td></tr><tr><td>ESP_OSR_05</td><td>-2,12</td><td>1</td><td>B</td></tr><tr><td>ESP_OSR_06</td><td>-1,67</td><td>1</td><td>B</td></tr><tr><td>ESP_OSR_07</td><td>-1,28</td><td>1</td><td>B</td></tr><tr><td>ESP_OSR_08</td><td>-0,32</td><td>2</td><td>B</td></tr><tr><td>ESP_APP_09</td><td>-0,64</td><td>2</td><td>B</td></tr><tr><td>ESP_APP_10</td><td>-1,12</td><td>1</td><td>B</td></tr><tr><td>ESP_APP_11</td><td>-1,79</td><td>1</td><td>B</td></tr><tr><td>ESP_APP_12</td><td>0,4</td><td>3</td><td>B</td></tr><tr><td>ESP_APP_13</td><td>-1,32</td><td>1</td><td>B</td></tr><tr><td>ESP_APP_14</td><td>-0,74</td><td>2</td><td>B</td></tr><tr><td>ESP_APP_15</td><td>-0,31</td><td>2</td><td>B</td></tr><tr><td>ESP_APP_16</td><td>-1,63</td><td>1</td><td>B</td></tr></tbody></table>
<table><tr><td>CHE_OSR_01</td><td>-0,05</td><td>2</td><td>C</td></tr><tr><td>CHE_OSR_02</td><td>0,53</td><td>3</td><td>C</td></tr><tr><td>CHE_OSR_03</td><td>0,57</td><td>3</td><td>C</td></tr><tr><td>CHE_OSR_04</td><td>0,05</td><td>2</td><td>C</td></tr><tr><td>CHE_OSR_05</td><td>0,38</td><td>3</td><td>C</td></tr><tr><td>CHE_OSR_06</td><td>1</td><td>4</td><td>C</td></tr><tr><td>CHE_OSR_07</td><td>1,41</td><td>4</td><td>C</td></tr><tr><td>CHE_OSR_08</td><td>0,43</td><td>3</td><td>C</td></tr><tr><td>CHE_APP_09</td><td>0,8</td><td>3</td><td>C</td></tr><tr><td>CHE_APP_10</td><td>1,32</td><td>4</td><td>C</td></tr><tr><td>CHE_APP_11</td><td>1,03</td><td>4</td><td>C</td></tr><tr><td>CHE_APP_12</td><td>0,81</td><td>3</td><td>C</td></tr><tr><td>CHE_APP_13</td><td>0,38</td><td>3</td><td>C</td></tr><tr><td>CHE_APP_14</td><td>0,13</td><td>2</td><td>C</td></tr><tr><td>CHE_APP_15</td><td>0,71</td><td>3</td><td>C</td></tr><tr><td>CHE_APP_16</td><td>0,51</td><td>3</td><td>C</td></tr><tr><td>ITA_OSR_01</td><td>1,66</td><td>4</td><td>C</td></tr><tr><td>ITA_OSR_02</td><td>0,73</td><td>3</td><td>C</td></tr><tr><td>ITA_OSR_03</td><td>0,26</td><td>3</td><td>C</td></tr><tr><td>ITA_OSR_04</td><td>0,15</td><td>2</td><td>C</td></tr><tr><td>ITA_OSR_05</td><td>1,3</td><td>4</td><td>C</td></tr><tr><td>ITA_OSR_06</td><td>0,71</td><td>3</td><td>C</td></tr><tr><td>ITA_OSR_07</td><td>2,07</td><td>4</td><td>C</td></tr><tr><td>ITA_OSR_08</td><td>0,77</td><td>3</td><td>C</td></tr><tr><td>ITA_APP_09</td><td>1,85</td><td>4</td><td>C</td></tr><tr><td>ITA_APP_10</td><td>1,54</td><td>4</td><td>C</td></tr><tr><td>ITA_APP_11</td><td>1,93</td><td>4</td><td>C</td></tr><tr><td>ITA_APP_12</td><td>0,27</td><td>3</td><td>C</td></tr><tr><td>ITA_APP_13</td><td>1,87</td><td>4</td><td>C</td></tr><tr><td>ITA_APP_14</td><td>1,21</td><td>4</td><td>C</td></tr><tr><td>ITA_APP_15</td><td>0,51</td><td>3</td><td>C</td></tr><tr><td>ITA_APP_16</td><td>0,91</td><td>3</td><td>C</td></tr><tr><td>GER_OSR_01</td><td>-0,79</td><td>2</td><td>C</td></tr><tr><td>GER_OSR_02</td><td>-0,07</td><td>2</td><td>C</td></tr><tr><td>GER_OSR_03</td><td>0,23</td><td>3</td><td>C</td></tr><tr><td>GER_OSR_04</td><td>0,34</td><td>3</td><td>C</td></tr></table>
<table><tr><td>GER_OSR_05</td><td>-0,91</td><td>2</td><td>A</td></tr><tr><td>GER_OSR_06</td><td>0,11</td><td>2</td><td>C</td></tr><tr><td>GER_OSR_07</td><td>-0,02</td><td>2</td><td>C</td></tr><tr><td>GER_OSR_08</td><td>-0,56</td><td>2</td><td>C</td></tr><tr><td>GER_APP_09</td><td>0,24</td><td>3</td><td>C</td></tr><tr><td>GER_APP_10</td><td>0,86</td><td>3</td><td>C</td></tr><tr><td>GER_APP_11</td><td>-0,75</td><td>2</td><td>C</td></tr><tr><td>GER_APP_12</td><td>0,49</td><td>3</td><td>C</td></tr><tr><td>GER_APP_13</td><td>0,18</td><td>2</td><td>C</td></tr><tr><td>GER_APP_14</td><td>1</td><td>4</td><td>C</td></tr><tr><td>GER_APP_15</td><td>0,4</td><td>3</td><td>C</td></tr><tr><td>GER_APP_16</td><td>-0,68</td><td>2</td><td>C</td></tr><tr><td>IRL_OSR_01</td><td>-1,34</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_02</td><td>-3,58</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_03</td><td>-2,35</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_04</td><td>-2,63</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_05</td><td>-0,91</td><td>2</td><td>A</td></tr><tr><td>IRL_OSR_06</td><td>-1,64</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_07</td><td>-2,81</td><td>1</td><td>A</td></tr><tr><td>IRL_OSR_08</td><td>-2,34</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_09</td><td>-2,41</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_10</td><td>-1,96</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_11</td><td>-0,91</td><td>2</td><td>A</td></tr><tr><td>IRL_APP_12</td><td>-1,33</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_13</td><td>-1,92</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_14</td><td>-1,69</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_15</td><td>-2,26</td><td>1</td><td>A</td></tr><tr><td>IRL_APP_16</td><td>-1,54</td><td>1</td><td>A</td></tr><tr><td>EST_OSR_01</td><td>2,32</td><td>4</td><td>D</td></tr><tr><td>EST_OSR_02</td><td>1,8</td><td>4</td><td>D</td></tr><tr><td>EST_OSR_03</td><td>1,8</td><td>4</td><td>C</td></tr><tr><td>EST_OSR_04</td><td>1,69</td><td>4</td><td>D</td></tr><tr><td>EST_OSR_05</td><td>-0,61</td><td>2</td><td>A</td></tr><tr><td>EST_OSR_06</td><td>1,34</td><td>4</td><td>C</td></tr><tr><td>EST_OSR_07</td><td>0,65</td><td>3</td><td>C</td></tr><tr><td>EST_OSR_08</td><td>1,83</td><td>4</td><td>D</td></tr></table>
<table><tr><td>EST_APP_09</td><td>-1,63</td><td>1</td><td>A</td></tr><tr><td>EST_APP_10</td><td>0,66</td><td>3</td><td>D</td></tr><tr><td>EST_APP_11</td><td>1,43</td><td>4</td><td>D</td></tr><tr><td>EST_APP_12</td><td>1,47</td><td>4</td><td>D</td></tr><tr><td>EST_APP_13</td><td>1,89</td><td>4</td><td>D</td></tr><tr><td>EST_APP_14</td><td>0,95</td><td>3</td><td>C</td></tr><tr><td>EST_APP_15</td><td>1,46</td><td>4</td><td>D</td></tr><tr><td>EST_APP_16</td><td>1,86</td><td>4</td><td>D</td></tr><tr><td>SWE_OSR_01</td><td>-1,48</td><td>1</td><td>B</td></tr><tr><td>SWE_OSR_02</td><td>-1,54</td><td>1</td><td>B</td></tr><tr><td>SWE_OSR_03</td><td>-1,09</td><td>2</td><td>B</td></tr><tr><td>SWE_OSR_04</td><td>-1,07</td><td>2</td><td>B</td></tr><tr><td>SWE_OSR_05</td><td>0,06</td><td>2</td><td>B</td></tr><tr><td>SWE_OSR_06</td><td>-2,05</td><td>1</td><td>B</td></tr><tr><td>SWE_OSR_07</td><td>-0,37</td><td>2</td><td>B</td></tr><tr><td>SWE_OSR_08</td><td>-0,46</td><td>2</td><td>B</td></tr><tr><td>SWE_APP_09</td><td>-0,36</td><td>2</td><td>B</td></tr><tr><td>SWE_APP_10</td><td>-1,26</td><td>1</td><td>B</td></tr><tr><td>SWE_APP_11</td><td>-0,72</td><td>2</td><td>B</td></tr><tr><td>SWE_APP_12</td><td>-1,25</td><td>1</td><td>B</td></tr><tr><td>SWE_APP_13</td><td>-1,36</td><td>1</td><td>B</td></tr><tr><td>SWE_APP_14</td><td>-1,24</td><td>1</td><td>B</td></tr><tr><td>SWE_APP_15</td><td>-1,19</td><td>1</td><td>B</td></tr><tr><td>SWE_APP_16</td><td>-1,44</td><td>1</td><td>B</td></tr></table>

417 The same dataset was summarised in three different ways. The numerical index provides a numerical value for each site on a continuum. The ordinal and the nominal indices provide four site clusters (continuous clusters 1, 2, 3 and 4 for the ordinal index and distinctive clusters A, B, C and D for the nominal index).

422
**Appendix D.** Table D.1. Contingency table with the number of sites reported for ordinal and nominal indices. The ordinal index classes are horizontal (from 1 to 4), the nominal index classes being vertical (A to D).

<table><thead><tr><th>ordinal ind.</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>nominal ind.</th><td></td><td></td><td></td><td></td></tr><tr><th>A</th><td>15</td><td>4</td><td>0</td><td>0</td></tr><tr><th>B</th><td>17</td><td>14</td><td>3</td><td>0</td></tr><tr><th>C</th><td>0</td><td>14</td><td>28</td><td>15</td></tr><tr><th>D</th><td>0</td><td>0</td><td>1</td><td>18</td></tr></tbody></table>

From this table, it follows that the nominal-cluster A was mainly associated with the ordinal-cluster 1 (low values of SBV and *N.ceranae*). The nominal-cluster B was associated with the ordinal-clusters 1 and 2 while the nominal-cluster C was associated with the ordinal-clusters 3 (high values of *N.ceranae*) and to a lesser extent with 2 and 4. The nominal-cluster D was associated with the ordinal-cluster 4 (high values of DWV-B and *N.ceranae*).

Author contribution statement:

Stéphanie Bougeard, Minh Huyen, Marion Laurent, Marie-Pierre Chauzat: Conceived and designed experiments, Performed experiment..</p>

Minh Huyen, Aurélie Babin, Eric Dubois, Marie-Pierre Rivière, Marion Laurent, Christine Druesne, Stéphanie Bougeard, Marie-Pierre Chauzat: Analysed and interpreted the data, Wrote the paper.</p>

Stéphanie Bougeard, Minh Huyen, Aurélie Babin, Eric Dubois, Christine Druesne, Marion Laurent, Marie-Pierre Chauzat: Contributed reagents, materials, analysis tools or data.</p>

Data availability statement:

The data that has been used is confidential.

Declaration of interest's statement:
456 The authors declare no competing interests.

457

458 Additional information:

459

460 No additional information is available for this paper.

461

462

463 References

464 1 Klein, A.-M. et al. Importance of pollinators in changing landscapes for world crops. *Proc.Roy.Soc.London Ser.B* **274**, 1-11, doi:10.1098/rspb.2006.3721 (2006).

465 2 Potts, S. G. et al. Safeguarding pollinators and their values to human well-being. *Nature* **540**, 220-229, doi:10.1038/nature20588 (2016).

468 3 Jacques, A. et al. A pan-European epidemiological study reveals honey bee colony survival depends on beekeeper education and disease control. *PLoS One* **12**, e0172591, doi:10.1371/journal.pone.0172591 (2017).

471 4 Potts, S. G. et al. Global pollinator declines: trends, impacts and drivers. *Trends Ecol Evol* **25**, 345-353, doi:10.1016/j.tree.2010.01.007 (2010).

473 5 Zavala-Cortés, A., Hernández, G. & Calderón-Salinas, J. V. An index for multidimensional assessment of swine health. *Trop. Anim. Health Prod.* **53**, doi:10.1007/s11250-020-02552-8 (2021).

476 6 Andraud, M. & Rose, N. Modelling infectious viral diseases in swine populations: A state of the art. *Porcine Health Management* **6**, doi:10.1186/s40813-020-00160-4 (2020).

479 7 Courtejoie, N., Zanella, G. & B., D. Bluetongue transmission and control in Europe: A systematic review of compartmental mathematical models. *Prev. Vet. Med.* **156**, 113-125, doi:10.1016/j.prevetmed.2018.05.012 (2018).

482 8 Rüegg, S. R., Heinzmann, D., Barbour, A. D. & Torgerson, P. R. Estimation of the transmission dynamics of Theileria equi and Babesia caballi in horses. *Parasitology* **135**, 555-565, doi:10.1017/S0031182008004204 (2008).

485 9 Sanchez-Vazquez, M. J., Nielen, M., Edwards, S. A., Gunn, G. J. & Lewis, F. I. Identifying associations between pig pathologies using a multi-dimensional machine learning methodology. *BMC Veterinary Research* **8**, doi:10.1186/1746-6148-8-151 (2012).

489 10 Alonso-Prados, E. et al. The toxic unit approach as a risk indicator in honey bees surveillance programmes: A case of study in *Apis mellifera iberiensis*. *Science of the Total Environnement* **698**, 134208, doi:10.1016/j.scitotenv.2019.134208 (2019).

492 11 Traynor, K. S. et al. Pesticides in honey bee colonies: Establishing a baseline for real world exposure over seven years in the USA. *Environ Pollut* **279**, 116566, doi:10.1016/j.envpol.2021.116566 (2021).

495 12 Abdi, H., Williams, L. J. & Valentin, D. Multiple factor analysis: principal component analysis for multitable and multiblock data sets. *WIREs Computational Statistics* **5**, 149-179, doi:10.1002/wics.1246 (2013).

498 13 Pagès, J. Analyse factorielle multiple appliquée aux variables qualitatives et aux données mixtes. *Revue de Statistique Appliquée* **50**, 5-37 (2002).

500 14 Husson, F., Lê, S. & Pagès, J. *Analyse de données avec R*., (Presses Universitaires de Rennes, 2016).
502 15 Benzécri, J. P. Histoire et préhistoire de l'analyse des données. Partie V L'analyse des correspondances. *Les Cahiers de l'Analyse des Données* **2**, 29 (1977).

504 16 Pagès, J. P., Cailliez, F. & Escoufier, Y. Analyse factorielle : un peu d'histoire et de géométrie. *Revue de Statistique Appliquée* **27**, 5-28 (1979).

506 17 Brown, M. et al. PoshBee: Pan-European Assessment, Monitoring, and Mitigation of Stressors on the Health of Bees. *ARPHA Preprints* **1**, ARPHA Preprints (2021).

508 18 Hodge, S. et al. Design and Planning of a Transdisciplinary Investigation into Farmland Pollinators: Rationale, Co-Design, and Lessons Learned. *Sustainability* **14**, 10549, doi:doi:10.3390/su141710549 (2022).

511 19 Babin, A. et al. Pan-european pathogen and parasite screening in the three sentinel bee species *A. mellifera*, *B. terrestris*, and *O. bicornis* upon exposure to crop fields *Paper submitted* (2023).

514 20 Babin, A., Schurr, F., Rivière, M.-P., Chauzat, M.-P. & Dubois, E. Specific detection and quantification of three microsporidia infecting bees, Nosema apis, Nosema ceranae, and Nosema bombi, using probe-based real-time PCR. *Eur. J. Protistol.*, 125935, doi:https://doi.org/10.1016/j.ejop.2022.125935 (2022).

518 21 Abdi, H. in *Encyclopedia of measurement and statistics* 907-912 (Thousand Oaks, Sage Publications, 2007).

520 22 Lê, S., Josse, J. & Husson, F. FactoMineR: An R Package for Multivariate Analysis. *Journal of Statistical Software* **25**, 1-18, doi:10.18637/jss.v025.i01 (2008).

522 23 Husson, F., Lê, S. & Pagès, J. Analyse de données avec R. (Presses Universitaires de Rennes, 2016).

524 24 Greenacre, M. J. *Theory and applications of correspondence analysis.*, (Academic Press, 1984).

526 25 de Tayrac, M., Lê, S., Aubry, M., Mosser, J. & Husson, F. Simultaneous analysis of distinct omics data sets with integration of biological knowledge: multiple factor analysis approach. *BMC Genomics* **10**, 32, doi:10.1186/1471-2164-10-32 (2009).

529 26 Best, A., White, A. & Boots, M. Maintenance of host variation in tolerance to pathogens and parasites. *Proceedings of the National Academy of Sciences* **105**, 20786-20791, doi:doi:10.1073/pnas.0809558105 (2008).

532 27 Costa, C., Tanner, G., Lodesani, M., Maistrello, L. & Neumann, P. Negative correlation between Nosema ceranae spore loads and deformed wing virus infection levels in adult honey bee workers. *J.Invertebr.Pathol.* **108**, 224-225, doi:10.1016/j.jip.2011.08.012 (2011).

536 28 Doublet, V., Natsopoulou, M. E., Zschiesche, L. & Paxton, R. J. Within-host competition among the honey bees pathogens Nosema ceranae and Deformed wing virus is asymmetric and to the disadvantage of the virus. *J. Invertebr. Pathol.* **124**, 31-34, doi:http://dx.doi.org/10.1016/j.jip.2014.10.007 (2015).

540 29 Gerard, M. et al. Stressful conditions reveal decrease in size, modification of shape but relatively stable asymmetry in bumblebee wings. *Sci. Rep.* **8**, 15169, doi:10.1038/s41598-018-33429-4 (2018).

543 30 Moher, D., Liberati, A., Tetzlaff, J. & Altman, D. G. Preferred reporting items for systematic reviews and meta-analyses: the PRISMA statement. *PLoS Med.* **6**, e1000097, doi:10.1371/journal.pmed.1000097 (2009).

546