Risk factors associated with honey bee colony loss in apiaries in Galicia, NW Spain

A cross-sectional study was carried out in Galicia, NW Spain, in order to estimate the magnitude of honey bee colony losses and to identify potential risk factors involved. A total of 99 samples from 99 apiaries were collected in spring using simple random sampling. According to international guidelines, the apiaries were classified as affected by colony loss or asymptomatic. Each sample consisted of worker bees, brood and comb-stored pollen. All worker bees and brood samples were analysed individually in order to detect the main honey bee pathogens. Moreover, the presence of residues of the most prevalent agrotoxic insecticides and acaricides was assessed in comb-stored pollen. The general characteristics of the apiaries and sanitary information regarding previous years was evaluated through questionnaires, while the vegetation surrounding the apiaries sampled was assessed by palynological analysis of comb-stored pollen. The colony loss prevalence was 53.5% (CI 95% = 43.2-63.9) and Nosema ceranae was found to be the only risk factor strongly associated with colony loss. The decision tree also pointed out the impact of the Varroa mite presence while variables such as apiary size, the incorrect application of Varroa mite treatments, and the presence of Acarapis woodi and Kashmir bee virus (KBV) were identified as possible co-factors.


Introduction
Recent years have seen a global decline in pollinator populations (FAO, 2008), including surprisingly large scale losses of honey bee colonies worldwide. However, to date no consensus has been reached regarding the cause or origin of this phenomenon Johnson et al., 2010;vanEngelsdorp & Meixer, 2010;Cepero et al., 2014). Even further, a recent analysis based on FAO data set (Moritz & Erler, 2016) indicates that global scale do not show a general colony decline, while stablishing some relationship with changes on political or socioeconomic systems, as well as finding that honey trade is a dominating factor for the number of managed colonies in a country (strong import/export ratio is related with those suffering higher colony declines). Anyway all through past decade numerous reports on big colony losses are frequent and have been attributed to pesticides, pathogens, as well as climate change, landscape alteration or agricultural intensification (Gonzalez-Varo et al., 2013).
One of the earliest and most studied hypotheses is related to the action of neonicotinoids and pheylpyrazole pesticides on bees (Smirle et al., 1984;Halm et al., 2006;Johnson et al., 2010;Staveley et al., 2014). As bees encounter these pesticides while foraging for pollen or nectar, comb-stored pollen could represent a chronic source of toxicity to bees (Thompson & Maus, 2007). Nevertheless, it has not yet been possible to validate this hypothesis in field conditions (Nguyen et al., 2009;Higes et al., 2010b;Cepero et al., 2014;Dively et al., 2015).
Another hypothesis suggests that different pathogens (e.g., viruses, microsporidia, acari, co-infections etc.) might produce a rapid decline in the adult bee population of the colony (Higes et al., 2008;Cox-Foster et al., 2007;Blanchard et al., 2008;Ravoet et al., 2013). Several studies have directly linked parasitisation by Nosema ceranae with colony loss (Martín-Hernández et al., 2007;Bacandritsos et al., 2010;Villa et al., 2013;Wolf et al., 2014) and indeed, infection with N. ceranae has been shown to induce sudden collapse of bee colonies under field conditions (Higes et al., 2008Wolf et al., 2014, Bekele et al., 2015. The role of Varroa destructor in the decline and death of colonies is also well established in conditions in which this mite is poorly controlled (Rosenkranz et al., 2010), or it may act synergistically with other pathogens (Cox-Foster et al., 2007;Ravoet et al., 2013). Recently, trypanosomatids were found to be a contributory factor to the colony losses in some areas (Runckel et al., 2011;Ravoet et al., 2013) but their pathogenic effect is not yet well described and there is no consensus on the different species parasiting honeybees (Cepero et al., 2014). Even nutritional stress due to habitat loss was also proposed as an alternative basis for honey bee colony collapse in the USA (Naug, 2009).
As there is no consensus, the most outlined and supported hypothesis to explain these colony losses is the interaction of multiple risk factors, including those of pathogens and pesticides (vanEngelsdorp & Meixer, 2010) environmental factors (Nguyen et al., 2009) and the effect of the climate change (Le Conte & Navajas, 2008).
The present study focused on the region of Galicia (in northwest Spain) in an attempt to identify risk factors correlated to colony losses, which have become a serious local beekeeping problem in recent years (Arnaiz, 2008).

Study design
The target population consisted in all the 2007's registered apiaries by the Spanish Ministry of Agriculture at the Galicia region and the number of samples to obtain was calculated in relation to the number of beekeepers registered, with an expected prevalence of colony loss (CL) of 67% , an absolute error of <10% and a confidence level of 95%. A total of 99 samples were collected in spring using simple random sampling.
The survey was carried out between May and July 2008. The 99 colonies sampled belonged to 99 different apiaries (one colony per apiary) and were classified on the basis of: apiary affected by CL or asymptomatic (No CL) during previous winter (Fig. 1).
The case definition criteria (CL) was defined as the mortality of honey bee colonies in the selected apiaries during previous autumn/winter with disappearance of adult bees from the colonies with little or no build-up of dead bees (inside the colony or in front of the colony entrance) with no prior symptoms of any disease (AFSA, 2008;Higes et al., 2010b). Information on CL status was gathered through a questionnaire to beekeepers (Higes et al., 2006) which also contained questions on general and sanitary aspects for each apiary.

Sample collection
Samples were collected directly by the veterinarian, who was the only person responsible for sampling to minimize errors. From each selected colony (n=99), five sub-samples were taken: worker bees (>100 from brood combs), forager bees (>30 collected closing the hive entrance for 1 min), sealed brood cells n> 400 (no mating stages of development), comb-stored pollen (>100 g) and honey from brood chamber combs (>100 g).

Analysis of main honey bee pathogens
Varroa loads determination. For each colony sample, adult worker bees, forager bees and brood samples were analysed separately for the presence of V. destructor using the recommended OIE methods (http://www.oie.int/es/normas-internacionales/ manual-terrestre).

Sample preparation for molecular analyses.
Biological material was weighed (Sartorius, ± 0.1 mg) and 15 mL of 50% solution of AL buffer (Qiagen 1014604) + 1 µg RNA Carrier/mL (Qiagen; equivalent to 2 µg carrier/ mL AL) was added. Samples were macerated for 2 min at high speed in a Stomacher® 80 Biomaster (Seward), using bags with a filter (BA6040/STR, Seward), and the homogenate was transferred to a separate tube. The macerates were centrifuged at 514 g for 10 min, recovering the supernatant for viral RNA extraction and the sediment for DNA extraction.
DNA extraction. Sediment (150 µL) was transferred to a 96-well plate (Qiagen) containing glass beads (2 mm diameter, Sigma). One blank (negative extraction control) was added for every 20 samples. Plates were shaken (Tissuelyser, Qiagen) at 30 Hz for 6 min, 30 µL of ATL buffer (Qiagen 19076) and 20 µL of Proteinase K (Qiagen 19131) was added to each well. Plates were incubated overnight at 56ºC. DNA was extracted using a BS96 DNA tissue extraction protocol in a BioSprint apparatus (Qiagen) and stored at -20ºC until required.

RNA extraction.
For viral RNA extraction, 400 µL of supernatant was incubated with protease (1 hour at 56ºC) and the RNA was extracted as described above (BioSprint, Qiagen). Immediately after extraction, cDNA was synthesised using reverse transcriptase (Transcrip HiFi cDNA Synth. Kit, Roche 05081963001) and stored at -20ºC until required.

Molecular analyses for pathogen detection.
The list of pathogens is summarized in. All analyses were performed using previously published methods.
To detect viruses, qRT-PCR assays as described previously (Mumford et al., 2000;Chantawannakul et al., 2006)  In all cases, negative and positive controls were included and analyzed in parallel to detect possible contamination and the good performance of the techniques.

Pesticide and palynological analysis
Pesticide analysis. To each comb-stored pollen sample a multi-residue analysis was performed as described previously (Higes et al., 2008; García-Chao et al., Table 1. Summary of the items included in the questionnaire used to identify the risk factors for colony loss, and the list of pathogens and pesticides tested.

General items related to the installations (n=7)
Location, geographic coordinates (UTM), number of apiaries, number of honey bee colonies, type of hives, beekeeping production and migratory activity.

Sanitary information related to the previous year (n=4)
Varroosis control: Product applied, number of treatments/year, duration of treatment (weeks) and doses employed.

Palynological analysis.
From each comb-stored pollen sample, a portion of approximately 0.5 g was taken to perform palynological determinations of 89 pollen types (Table 1) allowing to identify the vegetation present in the vicinity of the apiary and to assess the affinity of honey bee foragers for the surrounding vegetation.

Data analysis
The data from the questionnaire and the laboratory results were entered into a data-entry database developed in EpiInfo TM for Windows (CDC, USA), and they were visually checked for implausible values that could be traced back to the original questionnaires or laboratory reports for clarification. Information related to Varroa control was used to create the new variable "correct use", whereby each register was classified (yes/no) based on the dose (nº strips) and duration of treatment by product according to the manufacturer´s instructions (Table 2).
A descriptive analysis was performed to estimate prevalence and confidence interval (CI 95% ) of CL and the variables included in the study using SAS 9.4. This step allowed identifying variables with a large number of missing observations or a low variability that may be of little interest for further studies. Fisher's exact test and the chi-square test were used to evaluate the association between CL and the other variables studied to identify risk factors. Odds ratio (OR) were calculated for accurate interpretation of the risk factors involved.
All variables were also included in a multivariate analysis using a decision-tree. Statistical significance was set at p<0.05 and restriction to new nodes at n=20. All calculations were performed using SPSS version 22.0.

Results
On the basis of case definition, 53 out of the 99 samples belonged to apiaries classified by the veterinarian as affected by CL during previous year (53.5%, CI 95% =43.2-63.9; Fig. 1). The observed prevalence of CL difference was not statistically significant with the expected one (67%; χ 2 = 2.90; p=0.087).
The beekeeping practices in the apiaries were very homogeneous: 95.9% of hives were Langstroth shape,  Table 2). The most frequently applied commercial product was Apistan®, although 6.1% used homemade formulas (e.g., fluvalinate). The commercial products were correctly applied by 37.1% of beekeepers, and Apistan® was the product with the most correct posology (41.4%). The majority of beekeepers (98%) did not use fumagilin, an antibiotic prescribed against N. ceranae.
Palynological results of pollen samples evidenced a total of 76 types of flora detected in the brood chamber, being 94.3% (CI 95% =88.6-98.3) classified as wild vegetation ( Table 2). The detection of pesticides in stored pollen was rare (Table 3), with fluvalinate as the most commonly identified compound (13%; CI 95% =5.9-20.3). No fipronil nor its metabolites, imidacloprid or thiametoxan were detected in any of the analyzed samples. Statistically, none of the detected compounds was correlated to CL.
The presence of V. destructor in exterior and interior bees was correlated to the presence of DWV in brood samples (Pearson=0.273, p<0.01 and Pearson=0.310, p<0.01, respectively) and the presence of V. destructor  Only the presence of N. ceranae (p<0.0001) and KBV (p<0.0325) were found to be significantly related to the CL (p<0.0001). All samples with presence of N. ceranae (=38) belonged to CL positive apiaries. Positive KBV samples (=5) were few and this might be taken into account for further interpretations.
The decision tree confirmed N. ceranae as the main risk factor (χ 2 =53.527; p<0.0001) to CL (Fig. 2). However for the samples with no presence of N. ceranae, the Varroa loads in brood (χ 2 =11.363; p=0.001) was the main explanatory variable yet no relation was found (p<0.3402) from the univariate analysis. Finally the Varroa loads in worker bees (χ 2 =4.375; p=0.036) might explain the CL for those samples with no N. ceranae and no Varroa in brood.

Discussion
The conducted study aimed to clarify the relations between CL and the presence of pesticides and the pathogen prevalence in colonies from Galicia. The study area represents approximately 4% of the total number of bee hives officially registered in Spain in 2007, but it is one of the regions in which the CL phenomenon has been discussed most actively in a wide range of forums, including beekeeper associations, as well as scientific and official institutions. This intense debate prompted the design of the present study to estimate the extent of the phenomenon and to identify risk factors involved.
Although CL is a worldwide observed phenomenon, its impact is not reflected in the official censuses in any country (vanEngelsdorp & Meixer, 2010), mainly due to the rapid replacement of colonies by means of beekeeping practices. Thus, prior to this study it was difficult to estimate the real impact of colony loss in the area under study, although an earlier estimate was used as a reference . The observed CL prevalence was 53.5% (CI 95% =43.2-63.9), lower than initially expected . This difference might be explained with the different sample collection strategies employed in each study.  studied samples sent by beekeepers whose honey bee colonies were suffering from CL, potentially leading to an overestimation of the extent of CL. In the present survey, samples were collected randomly thus representing a more complete and balanced analysis of the general situation.
When analysing the possible risk factors involved in the CL phenomenon, we found that variables regarding beekeeping management such as hive type, orientation or migration activity were largely homogeneous among apiaries, and thus no link to CL was established. Only a small apiary size was associated with a lower prevalence of CL, although this association could not be established statistically probably due to the small number of apiaries of this type found.
Several pathogens were widely distributed among the apiaries studied. The ubiquity of DWV and BQCV, both with a prevalence of 94.9%, was consistent with studies performed in France (Tentcheva et al., 2004), USA (Chen et al., 2006) and Austria (Berényi et al., 2006). As found in previous studies, we also observed correlations between each of these two viruses and V. destructor, which has been proposed to act as a mechanical viral vector (Ball, 1989;Tentcheva et al., 2004) or as a trigger for viral replication through its immunosuppressive activity (Yang & Cox-Foster, 2005). Honey bee virus infections are usually latent, and their presence is associated with the presence of V. destructor and/or N. ceranae (Tentcheva et al., 2004;Chen & Siede, 2007). No other viruses, such as ABPV, SBV, AIV, IAPV or CPV were detected in the present study. IAPV was initially linked to colony losses in the USA (Cox-Foster et al., 2007) and considered a significant marker for Colony Collapse Disorder (CCD). However, in agreement with the present findings, no such association was established between this virus and the severe winter losses in other European countries (Blanchard et al., 2008;Higes et al., 2008;Higes et al., 2009;Garrido-Bailón et al., 2010b), suggesting that other pathogens play an important role in CL.
In this study, KBV was significantly linked to the presence of CL when univariate analysis was performed. Such association could not be confirmed by the multivariate analysis, and so suggesting that the role of KBV in this study might be analysed carefully. The small number of positive samples may cause a bias in the interpretation of this result.
After the viruses, V. destructor was the most common pathogen identified in the survey with prevalences ranging from 35.4% in brood samples to 49.5% in house adult bees. This high prevalence of V. destructor, a well known pathogen for which a variety of commercial treatments are available (Rosenkranz et al., 2010), may be related to the strikingly poor observance of correct treatment application. Analysis of V. destructor control strategies, which determined whether correct doses and treatments were applied, revealed a higher proportion of CL in incorrectly treated apiaries. This association could not be statistically confirmed and so we can only suggest that poor management of V. destructor might be a risk factor in the CL phenomenon. This statement is supported with the came out of the multivariate analysis that pointed out that both Varroa loads in brood and worker bees can explain some of the CL according to this statistical model.
The use of homemade formulas for V. destructor treatment in 6% of the apiaries studied is also worrying, as these probably contain unregistered acaricides or compounds that are ineffective against this pathogen. Furthermore, the identification of fluvalinate in the 13.1% of the comb-stored pollen samples suggests the undeclared use of homemade products by some beekeepers. Similar findings have been reported in other countries (Nguyen et al., 2009) and these observations are consistent with those of Johnson et al. (2010), who wrote: "Beekeepers searching for the primary source of pesticides contaminating bee hives need only to look in a mirror". Nosema ceranae was identified in 38.4% of the apiaries sampled, a lower prevalence than that previously reported in other regions of Spain in 2005 . However, such figures are not surprising nowadays and according to the Brooks (1979) criteria for insect microsporidia, they represent epizootic levels of nosemosis due to N. ceranae, the agent of the nosemosis type C, so named to differentiate it from the disease caused by another honey bee microsporidia, N. apis (Higes et al., 2006). Decision tree strongly linked N. ceranae with CL as 100% apiaries that presented this microsporidium did suffer CL phenomenon. No association was found between the microsporidium species and the viruses. The most notable is the case of BQCV (Chen & Siede, 2007), possibly due to the facilitation of transenteric viral infection through the cytopathogenic effects of Nosema, but it did not appear in this study. CL in N. ceranae-positive apiaries was significantly greater than in those free of this pathogen, with an N. ceranae prevalence of 100% in the unhealthy colonies and an associated 4-fold increase in the risk of CL. This result concords with and extends the earlier findings in this country . While there has been debate regarding the role N. ceranae plays in CL, recent research has reached consensus that under certain conditions N. ceranae causes collapse to honey bee colonies (Betti et al., 2014;Bekele et al., 2015;Cavigli et al., 2016).
The observed prevalence of A. woodi (16.2%) was higher than expected and it was significantly greater in apiaries in which Varroa treatments were incorrectly applied, further emphasizing the dangers of inadequate mite control. Few studies of the last 15 years have investigated the role of A. woodi on bee colonies mortality in Europe, although McMullan & Brown (2009) qualitatively modelled the effect of tracheal mite infestation in colony death.
Some authors have linked the presence of trypanosomatids with CL (Runckel et al., 2011;Ravoet et al., 2013). However they were found in half of adult bee samples, then more research is needed to clarify this implication (Cepero et al., 2014).
Palynological analysis revealed that the flora surrounding the sampled hives was diverse but mainly of wild origin rather than managed crops. No association between surrounding flora and CL was established. This observation is crucial, given that CL has been linked with nectar and pollen contamination by agropesticides used in crop management close to apiaries (Brodschneider & Crailsheim, 2010;Higes et al., 2010b). Neither imidacloprid nor fipronil were detected in comb-stored pollen samples, despite the low detection limits of the techniques used in our study (Nguyen et al., 2009), further implying that crop pesticide contamination does not influence the CL phenomenon in the study region.
With all the information gathered, N. ceranae is the main explanatory variable associated with CL for this study. This pathogen causes chronic affection of the honey bee colonies and has been identified worldwide as primarly cause of deaths. Pathological effects of N. ceranae in Spain had been associated with the increased virulence of the Spanish strain of N. ceranae (Huang et al., 2012), or even that Apis mellifera iberiensis might be more susceptible to infection. However, the development of N. ceranae infection, bee survival and spore loads were similar in French and Spanish bee strains (Dussaubat et al., 2013), and only bees from the Netherlands appear to survive better than Spanish bees (Van der Zee et al., 2014). However, recent results showed that other beekeeping factors may have an important role in Dutch bee survival, such as the use of oxalic acid for varroa treatment (Nannetti et al., 2015) or the beneficial beekeeping management of changing queens every season (Botias et al., 2012) to avoid N. ceranae acceleration of age polyethism of young bees, causing them to exhibit behaviours typical of older bees and early death (Lecoq et al., 2016). Other issues to take into account is that the prevalence and virulence of this pathogen is higher in temperate areas like Southern Europe , while its prevalence is lower in colder areas (Fries, 2010). This can be explained by the ability of its spores to survive high temperatures and desiccation better than low temperatures, and the ability of N. ceranae to complete its life cycle more efficiently at high temperatures (Fenoy et al., 2009;Sanchez-Collado et al., 2014).
The results of this study identify primary the presence of N. ceranae as the main explanatory variable associated with CL. However, the prevalence of V. destructor in brood cells and house bees also influences the presence of this phenomenon. The scenario found in this representative region of Spain, with a high prevalence of CL, infection with viruses, V. destructor and A. woodi, and the absence of crop pesticides, suggests that this phenomenon is related to the poor control of infectious agents rather than to toxic environmental effects. This work is in agreement with the previous observations in a temperate climate country of the lethal effect of pathogens such as Nosema ceranae and Varroa destructor and can explain the local beekeeping big problems detected in Northern Spain region at the beginning of the XXI century.