Holistic screening of collapsing honey bee colonies in Spain: a case study.

BACKGROUND: Here we present a holistic screening of collapsing colonies from three professional apiaries in Spain. Colonies with typical honey bee depopulation symptoms were selected for multiple possible factors to reveal the causes of collapse.
RESULTS: Omnipresent were Nosema ceranae and Lake Sinai Virus. Moderate prevalences were found for Black Queen Cell Virus and trypanosomatids, whereas Deformed Wing Virus, Aphid Lethal Paralysis Virus strain Brookings and neogregarines were rarely detected. Other viruses, Nosema apis, Acarapis woodi and Varroa destructor were not detected. Palinologic study of pollen demonstrated that all colonies were foraging on wild vegetation. Consequently, the pesticide residue analysis was negative for neonicotinoids. The genetic analysis of trypanosomatids GAPDH gene, showed that there is a large genetic distance between Crithidia mellificae ATCC30254, an authenticated cell strain since 1974, and the rest of the presumed C. mellificae sequences obtained in our study or published. This means that the latter group corresponds to a highly differentiated taxon that should be renamed accordingly.
CONCLUSION: The results of this study demonstrate that the drivers of colony collapse may differ between geographic regions with different environmental conditions, or with different beekeeping and agricultural practices. The role of other pathogens in colony collapse has to bee studied in future, especially trypanosomatids and neogregarines. Beside their pathological effect on honey bees, classification and taxonomy of these protozoan parasites should also be clarified.

Background

The beekeeping sector is suffering unexpected losses in many countries to such extent that pollination services in some cases are jeopardized. Several explanatory ‘drivers’ are known: an increasing number of pathogens, invasive species, exposure to pesticides, reduced genetic diversity and some apicultural practices. The driver ‘pathogens’ has received much attention, so scientists have sought for years to find the dangerous mixture of transmittable bee diseases. They often had to face conflicting data and it became increasingly clear that the involved pathogens may vary significantly in different regions.

The early detection of Acute Bee Paralysis Virus (ABPV) [1], Chronic Bee Paralysis Virus (CBPV) [2] and Deformed Wing Virus (DWV) [3] in collapsing colonies has determined the experimental design of many subsequent health studies. Later on, the set of target viruses has increased: so, Kashmir Bee Virus (KBV), Black Queen Cell Virus (BQCV) and Sacbrood Virus (SBV) became commonly examined [4]. An unbiased microbiome study aimed at finding the cause of the Colony Collapse Disorder (CCD) in the USA has extended this ‘short list’ with the Israeli Acute Paralysis Virus (IAPV) and the microsporidian parasite Nosema ceranae [5]. Further investigations either indicate [6, 7] or confirmed [8-11] the important role of N. ceranae in temperate areas of the world. Very recently, the trypanosomatid Crithidia mellificae was found to be a contributory factor to the colony losses in Belgium [12].

The driver ‘invasive species’ refers mainly to the ectoparasitic mite Varroa destructor. With respect to bee mortality, the Varroa-load seems to be one of the few decisive factors who stand across the borders. In the USA, the small hive beetle (Aethina tumida) seems to be another leading cause of mortality in beekeeping operations [13].

Among proposed causes of bee mortality, the exposure of bees to pesticides received much attention lately and to such extent that the European Commission adopted a proposal [14] to restrict the use of 3 pesticides belonging to the neonicotinoids family (clothianidin, imidacloprid and thiametoxam) for a two years period. However, nationwide monitoring programs of honeybees’ exposure level to these crop protection products are rather scarce. Moreover, the real involvement of these pesticides is controversial [15].

Although there is a general agreement that the bee mortality problem is multifactorial [13, 16], monitoring programs or case studies that go beyond screening of bee pathogens are rather limited. Besides, the few nationwide studies or clinical studies focusing on pathogens in combination with pesticide residue analysis [17-19] are restricted to the current ‘short list’ of pathogens (DWV, ABPV, V. destructor).

Here we present a holistic screening (a case study) of collapsing honey bee colonies from three professional Spanish apiaries with high colony losses. Colonies with typical depopulation were analysed for the presence of multiple putative drivers of collapse: honey bee viruses, Nosema spp. Varroa destructor, Acarapis woodi, trypanosomatids, neogregarines, neonicotionoid insecticides and foraging flora. This study could also reveal what causal factors should be included in future Spanish monitoring programs [20-28].

Results and discussion

The most widespread pathogens detected in the analyzed samples were N. ceranae and Lake Sinai Virus (LSV), found in all samples (100%, 10/10). Viruses of the ABPV complex, Acarapis woodi, Nosema apis, CBPV, SBV, SBPV and Varroa destructor were not detected (0%, 0/10) (Table 1).

Table 1

Results of the honey bee pathogen (viruses, parasites) screening and pollen analyses (pesticide residue; palynology) on samples from collapsing colonies in Spain

Sample information			Viruses				Parasites						Pesticide							Palinology
Sample ID	Apiary	Location	BQCV	DWV	LSV	ALPV	V. destructor	A. woodi	Trypanosomatids	Neogregarines	N. apis	N. ceranae	Dinotefuran	Nitenpyram	Thiametoxan	Clothianidin	Imidacloprid	Acetamiprid	Thiacloprid	Plant genus	Source
55	1	Guadalajara	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Raphanus, Rosmarinus, Salix, Brassica, Dorycnium	Wild spp.
56	1	Guadalajara	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Raphanus, Rosmarinus, Salix, Brassica, Dorycnium	Wild spp.
57	1	Guadalajara	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Raphanus, Rosmarinus, Salix, Brassica, Dorycnium	Wild spp.
58	1	Guadalajara	1	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Raphanus, Rosmarinus, Salix, Brassica, Dorycnium	Wild spp.
59	1	Guadalajara	1	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Raphanus, Rosmarinus, Salix, Brassica, Dorycnium	Wild spp.
1980	2	Vizcaya	1	1	1	1	0	0	1	0	0	1	0	0	0	0	0	0	0	Helianthemum, Raphanus, Rosmarinus	Wild spp.
324	3	Murcia	1	0	1	0	0	0	1	1	0	1	0	0	0	0	0	0	0	Thymus, Cistus, Salix, Lavandula	Wild spp.
325	3	Murcia	1	0	1	1	0	0	1	0	0	1	0	0	0	0	0	0	0	Thymus, Cistus, Salix, Lavandula	Wild spp.
328	3	Murcia	0	0	1	0	0	0	1	0	0	1	0	0	0	0	0	0	0	Thymus, Cistus, Salix, Lavandula	Wild spp.
329	3	Murcia	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	Thymus, Cistus, Salix, Lavandula	Wild spp.

0 = not detected; 1 = detected.

The high prevalence (100%) of N. ceranae in the present study confirmed earlier reports [25, 29]. Its ability to evoke the collapse of honey bee colonies and to create the symptoms described by veterinarians (see methods), have been previously reported [10, 18, 20, 30–34]. N. ceranae plays a controversial role in the worldwide colony losses phenomenon [16]. While in Mediterranean areas a direct link between this pathogen and the honeybee losses has been reported [6, 8, 9, 11, 30, 35, 36], it can be excluded as main cause of losses in colder areas or continental climates [17, 37–39]. Nevertheless, in Belgium it was found that the adverse effects of this microsporidian can be enhanced in combination with the trypanosomatid C. mellificae [12]. In the present study, they were found in 4 samples and 2 apiaries. The three isolates studied (324, 325 and 1980) displayed several haplotypes, with variable levels of intra-isolate diversity (Table 2). Overall, 18S rDNA was less variable (pooled total diversity π = 0.16 ± 0.06%; average ± SE) than GAPDH (π = 0.40 ± 0.14%) although the difference between both loci was not statistically significant. C. mellificae ATCC30254 also presented multiple haplotypes at the GAPDH locus, which reflects that this reference strain is not isogenic.

Table 2

Estimates of Trypanosomatid diversity at all sites ( ) expressed as percentage

Isolate	18S rDNA		GAPDH
	Average π	SE	Average π	SE
ATCC30254	0.00	0.00	0.68	0.23
324	0.06	0.06	0.86	0.35
325	0.30	0.11	0.29	0.11
1980	0.08	0.08	0.39	0.16
Pooled	0.16	0.06	0.40	0.14

SE standard error, calculated by a bootstrap procedure (3000 replicates) with MEGA5.

The pairwise comparison of the 18S rDNA sequences revealed that C. mellificae ATCC30254 (KJ704242 – KJ704251) exhibited a single mutation that was not present either in the presumed C. mellificae sequences deposited in GenBank (KF607064.1 and AB745488.1) or in ours (KJ704218 – KJ707241). This single difference between C. mellificae ATCC30254 and the rest of the sequences is extremely important since any nucleotide variant in a highly conserved marker like the 18S rDNA gene, which displays strong identity (about 99%) among different genus (Cepero et al., submitted), suggests that C. mellificae ATCC30254 and our sequences (which are mostly identical to KF607064.1 and AB745488.1; Additional file 1: Table S1) might represent genetically isolated organisms.

The analysis of the genetic distances between GAPDH sequences showed that C. mellificae ATCC30254 (KJ704273 – KJ704282) displayed a 6.89 ± 1.32% divergence with respect to the rest of the sequences, that included both the presumed C. mellificae sequences deposited in GenBank (AB716357.1, AB745489.1 and JF423199.1) and ours (KJ704252 – KJ704272), which, again, were in a large fraction identical to AB716357.1 and AB745489.1 (Additional file 2: Table S2). However, this finding is even more surprising if we bear in mind that total divergence estimates, as those calculated here, underestimate the neutral genetic distance between sequences (as replacement sites account for nearly 75% of the coding sequence and their evolutionary rate is severely limited by purifying selection). Consequently, to obtain more reliable estimates of divergence we performed pairwise comparisons at synonymous sites, which are considered neutral or nearly neutral. The outcome of this analysis is that there is a large genetic distance (23.47 ± 5.48%, Table 3) between C. mellificae ATCC30254, an authenticated cell strain since 1974, and the rest of the presumed C. mellificae sequences. This means that the latter group corresponds to a highly differentiated taxon that should be renamed accordingly. It is also worth noting the dispersal of Crithidia sp. sequences all over the tree (Figure 1), which questions their current taxonomic classification. In line with this, highly divergent C. mellificae lineages were previously mentioned [40].

Table 3

Pairwise estimates of synonymous substitutions per synonymous site between sequences

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
1
2	0
3	0	0
4	1.58	1.58	1.58
5	0	0	0	1.58
6	1.60	1.60	1.60	0	1.60
7	0	0	0	1.58	0	1.60
8	1.58	1.58	1.58	0	1.58	0	1.58
9	0	0	0	1.58	0	1.60	0	1.58
10	1.60	1.60	1.60	0	1.60	0	1.60	0	1.60
11	24.66	24.66	24.66	26.77	24.66	26.30	24.66	26.77	24.66	26.30
12	21.80	21.80	21.80	23.83	21.80	23.36	21.80	23.83	21.80	23.36	3.23
13	23.21	23.21	23.21	25.29	23.21	24.81	23.21	25.29	23.21	24.81	4.34	1.06
14	23.21	23.21	23.21	25.29	23.21	24.81	23.21	25.29	23.21	24.81	1.06	2.14	3.23

1: AB745489.1, 2: AB716357.1, 3: JF423199.1, 4: KJ704252 – 3, 5: KJ704254, 6: KJ704255, 7: KJ704256 – 63, 8: KJ704264, 9: KJ704265 – 69, KJ704271 – 72,10: KJ704270, 11: KJ704273, KJ704275, 12: KJ704274, KJ704281 – 82, 13: KJ704276 – 77, KJ704279, 14: KJ704278, KJ704280. Analyses were conducted in MEGA5 using the Nei-Gojobori model. All positions containing gaps and missing data were eliminated. Comparisons involving ATCC30254 and presumed C. mellificae sequences are highlighted in bold.

Figure 1

Map of Spain. Provinces of origin of the samples. In the north, Vizcaya. In the central area, Guadalajara. In the southeast, Murcia.

An unbiased metagenomic study in a Spanish non-professional apiary with collapsing colonies revealed the presence of several viruses, among which LSV complex and ALPV strain Brookings [41]. In this study, the former was present in all samples, and the latter in only a few (20%, 2/10). Sequencing of the LSV amplicons revealed that the strains from apiary 1 and 2 were almost identical (Genbank: KJ561228, KJ561229), but had a low resemblance to the strain from apiary 3 (Genbank: KJ561227). This strain has a high amino acid similarity (89%) with the Orf1 of LSV strain Navarra (Genbank: AGF84788). So, different LSV strains are present in Spain similar to the situation in Belgium [12] and the USA [40, 42]. The pathogenic implications of LSV in honey bee health status are under discussion. Cornman et al. [40] suggested a potential association between LSV and CCD colonies, although this observation was not found in a Belgian bee health screening [12]. In our work, LSV represents the most abundant virus in the analyzed samples collected during spring and summer. For these reason, the importance of LSV and the pathogenicity of different LSV strains should be further investigated.

Black Queen Cell Virus (BQCV) was found in 50% of samples (5/10) and in all apiaries. This result corresponds with those obtained before [24]. Aphid Lethal Paralysis Virus (ALPV) strain Brookings was found in 20% of samples (2/10) in two apiaries. The amino acid sequences were identical (>98%) to those detected previously in Belgium [12], Spain [41] and the USA [42] (Genbank: AGU62863, AGF84786, AEH26191). Curiously, Deformed Wing Virus (DWV) was only detected in one sample in one apiary, and this might be related to the low prevalence of V. destructor in the analyzed samples and therefore apiaries.

Neogregarines were detected in one sample. Direct sequencing indicated an infection by Apicystis bombi. However, the presence of overlapping peaks at particular points of the electropherograms suggested a potential mixture of templates, which was further investigated by cloning and sequencing. This process yielded sequences from another Apicomplexan parasite (99% identity with Eimeriidae or Cryptosporidiidae in Blastn), which was co-infecting the colony with A. bombi (whose presence was only confirmed by direct sequencing). Although A. bombi was thought to be mainly a bumblebee parasite [43, 44], there are increasing findings of that parasite in honey bees [45-47] since a molecular detection method became available [44]. Nevertheless, the amplification of other parasites with the same primers should be taken into account to avoid misdiagnosis in future studies.

Palinological analyses confirmed that the honeybees foraged on wild plants. As a consequence, no neonicotinoid residues in stored pollen were detected. Although these pesticides exhibited severe acute and sublethal effects on bees [48-50], their role as sole cause of colony loss is still not clear [13, 51]. Indeed, most published data have shown that acaricides, herbicides or other insecticide molecules, were more prevalent than neonicotinoids or phenilpirazoles in hives [22, 23, 52–55]. Our results confirm that neonicotinoids are not the only cause capable of causing the colony collapse of honey bee colonies.

Conclusion

Many predictive markers and drivers have been suggested for honey bee colony collapses [5, 13, 56, 57]. The collapses of honey bee colonies in our study were not related with the presence of neonicotinoids or V. destructor. Instead, N. ceranae seems to be the main culprit of the colony losses in this study as already suggested in previous investigations [20, 30]. The role of other pathogens in colony collapse has to be studied in future, especially trypanosomatids and neogregarines. Beside their pathological effect on honey bees, classification and taxonomy of these protozoans also should be clarified. The results of this study clearly demonstrate that the drivers of colony collapse may differ between different geographic regions (see 16, 40).

Methods

Samples were collected from three professional apiaries each located in a different region of Spain, all under different climatic and environmental conditions (Figure 2). The veterinarians responsible for these apiaries contacted the Centro Apícola Regional (CAR) pathology laboratory, because of alarming symptoms like depopulation and unusually high colony losses.

Figure 2

Phylogenetic analysis of GAPDH sequences of trypanosomatids. The tree was constructed using the Maximum Likelihood method (based on the General Time Reversible model, assuming a Gamma distribution and allowing some sites to be evolutionarily invariable). Bootstrap values (500 replicates) are shown next to the branches. Codon positions included were 1st + 2nd + 3rd + Noncoding. All positions containing gaps and missing data were eliminated. The clade of presumed Crithidia mellificae sequences, both from GenBank (AB745489.1, AB716357.1, JF423199.1) and from this work (KJ704252- KJ704272), was compressed for ease of visualisation and is displayed in bold. The same was done for C. mellificae ATCC30254 sequences (KJ704273 – KJ704282).

Apiaries

Apiary 1: It is located in the Center of Spain in the province of Guadalajara and consisted of 400 bee colonies wintered in 2012, distributed in five apiaries all sited in the same province. In early spring 2013, 70% of the colonies had collapsed, with clear symptoms of depopulation such as disappearance of adult bees and unattended brood. Chalkbrood and foulbrood were not reported by the responsible veterinarian. The 30% of surviving colonies in early spring had a small adult bee population and low vitality. To perform the analysis, 25 weak surviving bee colonies were sampled (worker bees and pollen) in August 2013.

Apiary 2 2: It is located in the province of Vizcaya, Northern Spain. In winter 2012, it had 300 bee colonies distributed in 3 apiaries very closed each other. In early spring 2013, all colonies showed clear signs of depopulation and imminent risk of collapse. The responsible veterinarian has not reported other disease symptoms, although the bee colonies showed a very marked lack of vitality. For analysis, ten collapsing colonies were sampled (worker bees and pollen) in August 2013.

Apiary 3: It is located in the Southeast of Spain, in the province of Murcia. This apiary wintered 175 bee colonies distributed in 4 apiaries in 2011 and the 90% of them had collapsed in early spring 2012. For the analysis, 16 weak surviving colonies were sampled (worker bees and pollen) in May 2012.

Honey bee sampling

Samples of worker honey bees (n = 100 worker bees) and stored pollen, were collected from the brood chamber by the veterinarians in charge and sent to the CAR bee Pathology Laboratory in spring (samples from Vizcaya and Murcia) o early summer (samples from Guadalajara). Samples from the apiary 1 (n = 25; five per apiary) were pooled together to obtain five pooled samples of five colonies each, one pool sample per apiary. Samples from apiary 2 (n = 10) were pooled into a single sample. Samples for apiary 3 (n = 16; four per apiary) were processed as the first apiary so there were four pooled samples from four colonies each (one pooled sample per apiary).

Varroa destructordetection

The presence of V. destructor in worker honey bee samples were analyzed as described previously [19, 20].

DNA and RNA extraction

All bee samples were macerated in AL buffer 50% (Qiagen) as previously described [24] using sterile bags with filter. Resulting pellets were used for DNA extraction and supernatants for RNA extraction. They were frozen at −80°C until extraction of the nucleic acids.

For DNA extraction, the pellet was resuspended in 3 ml sterile water and 400 μl was transferred into a 96-well plate (Qiagen) with glass beads (2 mm diameter, Sigma). Samples were then processed as previously described [25]. Briefly, the plates were shaken for 6 min at 30 Hz. Afterwards, 150 μL of each sample was transferred to a Deepwell plate (Eppendorf) with 30 μl of ATL buffer (Qiagen) and 20 μl of Proteinase K (Qiagen). After overnight incubation at 56°C, DNA was extracted using the BS96 DNA Tissue extraction protocol in a BioSprint machine (Qiagen). Plates were stored at −20°C until use.

RNA was extracted from the supernatant using the DNeasy Blood & Tissue kit (Qiagen), according to the manufacturer’s instructions for RNA isolation. Briefly, 200 μl PBS and 1 μl carrier RNA were added to 220 μl of the resulting supernatants. After 10 min incubation at 56°C, 230 μl ethanol was added. Subsequent to binding and washing, RNA elution was accomplished with 50 μl nuclease-free water.

The C. mellificae reference strain ATCC30254 was included in the analysis for sequence comparison with field isolates. This strain was first cultivated as recommended in ATCC medium 355. Further sub-cultivation was performed in Brain Heart Infusion (BHI) medium as described by Popp and Lattorf [58]. Visible and isolated colonies were taken and resuspended in milliQ water (PCR-quality) and the processed for DNA extraction as above described.

PCR and MLPA analysis

A broad pathogen screening of the honeybee samples was performed using published PCR assays for Acarapis spp. mites [27], Nosema apis and N. ceranae [25], trypanosomatids and neogregarines [44]. Given that the 18S rDNA gene alone is unsuitable to classify trypanosomatid/s infecting honeybee colonies (Cepero et al., submitted), these were also amplified at the GAPDH locus using CrGD-1 F (5′ GGTCGCCGTGGTGGAC 3′) and CrGD-1R (5′ CGTCGCCGTGTAGGAGTGA 3′). Since these oligos did not produce an amplicon from ATCC30254, this strain was PCR-amplified with Tryp-1 F (5′ CCGAGTACTTCKCSTACCAG 3′) and Tryp-1R (5′ AGCCGAGGATGCCCTTCAT 3′), a set of primers that should amplify most Crithidia and Leptomonas species published in GenBank to date. As template, we used 5 μl DNA in each reaction.

For virus analysis, the BeeDoctor test, a ‘multiplex-ligation probe dependent amplification’ (MLPA) based method capable of detecting CBPV, DWV-complex, ABPV-complex, BQCV, SBPV and SBV, was performed as described before [59], starting from 3 μl RNA. The amplified MLPA products were analyzed using 4% high resolution agarose gel electrophoresis.

In addition, we performed additional RT-PCR analyses for few viruses. Using random hexamer primers, 500 ng RNA was retro-transcribed with the RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific). For ALPV strain Brookings and LSV complex detection, we used 1 μl cDNA template in the PCR test described by Runckel et al., [42] and Ravoet et al. [12] respectively.

Sequencing and cloning

Positive samples of ALPV strain Brookings, LSV complex, trypanosomatids and neogregarines were re-analyzed using the HotStar HiFidelity Polymerase Kit (Qiagen) (viruses) or the Expand High FidelityPLUS PCR System (Roche Diagnostics) (protozoans).

Protozoan amplicons were cloned using the Topo TA Cloning Kit (Invitrogen) after purification with the QIAquick PCR Purification Kit (Qiagen). Plasmids were purified using the QIAprep Spin Miniprep Kit (QIAgen). Both strands were sequenced with M13 primers on an automated ABI3730XL sequencer using Big Dye (Applied Biosystems). ALPV strain Brookings and LSV complex amplicons were bidirectional sequenced on an ABI3130XL with gene-specific primers.

Sequences were checked for accurate base calling using CodonCode Aligner (CodonCode Corporation) and their identity verified by means of a nucleotide BLAST. Alignments were manually edited with BioEdit [60] and submitted to GenBank (KJ704218 – KJ704282).

Estimates of overall and pairwise nucleotide diversity were obtained with MEGA v.5.05 [61]. These were measured as π, applying the Jukes and Cantor (JC) correction [62], and their standard error (SE) calculated by a bootstrap procedure (3000 replicates). This software was also used to construct a maximum likelihood phylogeny under the GTR + G + I model (General Time Reversible + Gamma distributed + Invariable sites) that was selected as best model by applying the Akaike information criterion (AIC), as implemented in MEGA. The reliability of the tree topology was tested by bootstrap support (500 replicates).

Pollen sampling

In order to obtain a representative pollen samples from each colony, stored pollen (3 squares of approximately 10 × 10 cm each) collected from the brood chamber combs were extracted aseptically from five colonies of each apiary [20, 22]. Each pollen sample was divided in two aliquots. One aliquot of 100 g was used for neonicotinoid screening, while the other aliquot of 5–10 g was used for palinological analysis. Samples were stored at −20°C until further use.

Neonicotinoids analysis

Seven neonicotinoid insecticides (acetamiprid, clothianidin, dinotefuran, imidacloprid, nitenpyram, thiacloprid, and thiamethoxam) were analyzed pollen with a previous methodology developed by Yáñez et al. [63]. It should be mentioned that the equipment, methods and reagents were the same than the described in this previous research. The sample treatment consisted of a solid–liquid extraction of the neonicotinoids from pollen with dichloromethane, followed by evaporation and reconstitution. Briefly, 2 g pollen sample and 10 mL of dichloromethane were transferred to a centrifuge tube. The mixture was mechanically shaken for 10 min at 800 oscillations per minute in a Vibromatic and then centrifuged for 10 min at 25°C and 10,400 × g. Following this, the supernatant was collected, filtered through paper filter, transferred to a 25-mL conical flask, and then evaporated until dry in a rotary evaporator at 40°C. The dry extract was reconstituted with 1 mL of a water and acetonitrile mixture (50:50, v/v) and the resulting solution was passed through a syringe filter, after which a 15-μL aliquot was injected into a liquid chromatograph (LC) coupled to an electrospray ionization mass spectrometry detector (ESI-MS). Once the neonicotinoids were extracted, they were determined using an optimized LC-ESI-MS method, which was validated in terms of selectivity, linearity, precision and recovery. The limits of detection (LOD) and quantification (LOQ) were 0.4–2.8 μg/kg and 1.2–9.1 μg/kg, respectively, and the extraction recoveries were between 86 and 106% in all cases.

Palinological analysis

The palinological analysis was performed as described previously [20, 22]. Briefly, pollen grains were isolated from each sample and cleaned up by the Erdtman method [64]. Species identification was performed using a photographic atlas [65, 66] and the pollen slides reference collection at the CAR. Briefly, pollen were extracted by diluting 0.5 g in 10 ml of acidulated water (0.5% sulphuric acid) and then centrifuged at 2,500 rpm for 15 min. The pellet was washed with double-distilled (dd) water, centrifuged and resuspended in ddH2O. 200 μl of this suspension was placed onto a glycerin jelly slide and examined microscopically in order to identify the pollen.

Ethics statement

In Europe, the EU Directive 2010/63/EU on the protection of animals used for scientific purposes laid the down the ethical framework for the use of animals in scientific experiments. The scope of this directive also includes specific invertebrate species, i.e. cephalopods, but no insects. Thus, according to European legislation no specific permits were required for the described studies.

Additional file 1: Table S1: Pairwise estimates of evolutionary divergence between 18S rDNA sequences (expressed as %). (DOC 176 KB)

Additional file 2: Table S2: Pairwise estimates of evolutionary divergence between GAPDH sequences (expressed as %). (DOC 190 KB)