DOI: 10.13102/sociobiology.v65i4.3442Sociobiology 65(4): 558-565 (October, 2018) Special Issue Open access journal: http://periodicos.uefs.br/ojs/index.php/sociobiology ISSN: 0361-6525 Searching for Molecular Markers to Differentiate Bombus terrestris (Linnaeus) Subspecies in the Iberian Peninsula Introduction The average pollination dependency of crops was estimated to around 300 000 million of euros in 2009 (Lautenbach et al., 2012). Among the insect pollinators, the bumblebee Bombus terrestris (Linnaeus) is used in greenhouses all over the world because they are better pollinators than honeybees for crops as strawberry, tomato or melon, due to their large body and buzzing capacity. The economic value of this service, just for tomato, has been measured in 12 000 million of euros per year (Velthuis & Doorn, 2006). Nevertheless, in the current scenario of bee decline, the introduction of non-native Bombus species and subspecies without a proper legislation may affect local bumblebee biodiversity (Goulson et al., Abstract Bumblebees (genus Bombus Latreille) are pollinator insects of great ecological and economic importance, which commercial use for pollination has increased since the 80s. However, the introduction of foreign Bombus terrestris (Linnaeus) has resulted in a decline of native bumblebee populations in Japan, Chile and Argentina among others. To study the potential introgression of commercial B. terrestris into the Iberian endemic subspecies Bombus terrestris lusitanicus Krüger it is necessary to find a precise molecular marker that differentiates both subspecies. For this purpose, comparative analyses were carried out between B. t. lusitanicus and Bombus terrestris terrestris (Linnaeus) from Spain and from Belgium by sequencing the nuclear genes elongation factor 1-α and arginine kinase and the mitochondrial gene 16S ribosomal RNA, and genotyping with eleven microsatellite loci. No differentiation was observed at the nuclear level, but haplotypes found within the 16S sequence correlated with the morphological characterization of the subspecies. In a case study including individuals sampled before the establishment of bumblebee rearing companies and others from recent samplings, we detected hybrid individuals (those with non-matching morphological subspecies and 16S haplotype) more frequently in the south supporting the naturalization of commercial B. t. terrestris and introgression events between both subspecies. This marker should be used in Iberian populations with the aim to support management and conservation actions in endemic populations of B. t. lusitanicus. Sociobiology An international journal on social insects D Cejas1, C Ornosa2, I Muñoz1, P De la Rúa1 Article History Edited by Cândida Aguiar, UEFS, Brazil Received 10 May 2018 Initial acceptance 13 June 2018 Final acceptance 03 August 2018 Publication date 11 October 2018 Keywords Bumblebees, genetic diversity, ArgK, EF1, 16S, mitochondrial DNA, Iberian Peninsula, microsatellites. Corresponding author Pilar De la Rúa Área de Biología Animal Departamento de Zoología y Antropología Física Facultad de Veterinaria Universidad de Murcia 30100 Murcia, Spain. E-Mail: pdelarua@um.es 2015). Such introductions suppose a risk for the conservation of endemic species and subspecies in many countries (Lecocq et al., 2015) to the extent that invasions of commercial non- native bumblebees have been reported as one of the 15 emerging issues for global conservation in 2017 (Sutherland et al., 2016). In this sense, Bombus terrestris (Linnaeus) can leave greenhouses and colonize the environment (Kraus et al., 2010), presenting some negative effects on native bees, for instance, displacing wild species while competing for resources such as pollen or nesting places (Matsumura et al., 2004; Aizen et al., 2018). This colonization might affect native species by spreading exotic diseases and parasites (Whitehorn et al., 2013) and changing plant-pollinator interactions in non-native environments, impacting crops, native plants 1 - Departamento de Zoología y Antropología Física, Facultad de Veterinaria, Universidad de Murcia, Murcia, Spain 2 - Departamento de Biodiversidad, Ecología y Evolución, Facultad de Ciencias Biológicas, Universidad Complutense, Madrid, Spain RESEARCH ARTIClE - BEES mailto:pdelarua@um.es Sociobiology 65(4): 558-565 (October, 2018) Special Issue 559 and pollinators (Morales et al., 2013; Aizen et al., 2014). B. terrestris have also better competitive capacities such as large foraging ranges and a broad diet, and they emerge early in the season making them adaptable to different environments (Matsumura et al., 2004). Currently, the invasive distribution of B. terrestris is increasing due to direct human intervention. Bumblebee introduction has impacted local bee populations in countries as Chile, China, Israel, Japan, Mexico, South Africa, South Korea, Taiwan, and New Zealand, where they are displacing native bee species to the edge of local extinction (Inoue et al., 2008; Schmid-Hempel et al., 2014; Acosta et al., 2016). The effect of reared B. terrestris in areas with consubspecific populations has been less investigated (Lecocq et al., 2016). B. terrestris presents nine subspecies classified by their body hair color pattern and distribution range (Rasmont et al., 2008). Among them, Bombus terrestris terrestris (Linnaeus) and Bombus terrestris dalmatinus Dalla Torre are the two most widely used in artificial rearing (Velthuis & Doorn, 2006). The endemic subspecies Bombus terrestris lusitanicus Krüger inhabits the Iberian Peninsula and reaches southern France where there is a natural contact zone with B. t. terrestris (Ornosa & Ortiz-Sánchez, 2004). During the last years, B. t. terrestris has been also detected at the south of the Peninsula (Ornosa, 1996; Vargas et al., 2013), mainly in Almeria where more than 30 000 hectares of greenhouses are located. We hypothesize that the contact between the endemic and the introduced subspecies in this region leads to introgression events that may yield to a loss of genetic diversity, or even displacement of the local endemic subspecies by individuals of the commercial subspecies. The taxonomy of these two B. terrestris subspecies has been discussed before, but neither the barcoding approach (based on mitochondrial cox1 gene fragment variation) nor the cephalic labial gland secretions have been able to differentiate them (Coppée et al., 2008; Williams et al., 2012; Lecocq et al., 2016). In this work, we aim to evaluate the sequence variation of two nuclear genes and one mitochondrial gene as potential subspecific markers to differentiate B. t. terrestris and B. t. lusitanicus. Nuclear arginine kinase (ArgK) and elongation factor-1α (EF-1α) genes include introns that have been proven useful to differentiate Bombus species (Kawakita et al., 2003) and to characterize geographic variations (Hines et al., 2006). The mitochondrial 16S ribosomal RNA (16S) gene has a higher substitution rate than barcoding fragment cox1, what results useful to lower-level (genera and species) analysis on the Hymenoptera (Rokas et al., 2002; Hines et al., 2006; Cameron et al., 2007). Furthermore, we have also used microsatellite loci developed by Estoup et al. (1993, 1996) adequate to analyze the genetic structure of bumblebee populations (Kraus et al., 2010, Dreier et al., 2014, Moreira et al., 2015). Material and Methods Sampling Preliminary sampling For the first experiment searching for molecular markers to differentiate B. terrestris subspecies in the Iberian Peninsula, we included 10 female individuals of B. t. lusitanicus and 20 B. t. terrestris sampled in 2011-2015. B. t. lusitanicus individuals were sampled along the bank of Duero River (Soria, Spain), National Park (N. P.) of Guadarrama (Madrid, Spain) and the Sierra Nevada N. P. (Granada, Spain). B. t. terrestris individuals were sampled also in Sierra Nevada N. P. and in Belgium (Table S1). The species of every individual was confirmed through barcoding (Murray et al., 2008), whereas the subspecies was determined through morphological characters. Individuals were classified in three groups according to their subspecies and geographic origin: TLS for B. t. lusitanicus from Spain, TTS for B. t. terrestris from Spain and TTB for B. t. terrestris from Belgium (Table S1). Samples were maintained in absolute ethanol and stored at -20 ºC in the laboratory until processed. Case study In order to confirm the discriminative subspecific power of the 16S haplotypes we designed a case study by adding 123 individuals to the prior sampling: 48 B. t. lusitanicus individuals collected in Spain between 1977 and 1985 before the settlement of bumblebee rearing companies plus 40 B. t. lusitanicus from different Spanish locations, 29 B. t. terrestris (21 from the north of France and eight from Sierra Nevada N. P., Spain) and six individuals considered morphologically as hybrids from the north of Spain (Huesca) and surroundings of Sierra Nevada N. P. (Spain). These 75 individuals were collected between 2013 and 2015 (Table 1 and Table S2 in Supplementary Material) and their species status was also confirmed through barcoding (Murray et al., 2008) Table 1. Subspecies, number of Bombus terrestris individuals (N), country and years of sampling in the case study. Distribution of 16S haplotypes detected is also shown. Subspecies N Country Years 16S haplotype A G B. t. lusitanicus 48 Spain 1977/1985 1 47 50 Spain 2013/2015 7 43 Total 98 8 90 B. t. terrestris 18 Spain 2013/2014 7 11 10 Belgium 2013 10 - 21 France 2015 21 - Total 49 38 11 Hybrids 6 Spain 2014 2 4 D Cejas; C Ornosa; I. Muñoz; P. De la Rúa – B. terrestris subspecies molecularly characterized in Spain560 DNA extraction, amplification and sequencing Total DNA was extracted from a hind leg of every individual following the protocol of Ivanova et al. (2006). As a first step, the tissue was digested with proteinase K during six hours at 56 ºC in continuous shaking. PureTaqTM Ready- To-Go PCR beads (GE Healthcare) were used for DNA amplification. PCR profile consisted on initial denaturation at 94 ºC for 3 min followed by 36 cycles of denaturation at 94 ºC for 1 min, annealing temperature at 48 ºC (ArgK), 53 ºC (EF-1α) and 49 ºC (16S) for 1 min, elongation at 68 ºC (72 ºC for EF-1α) for 1 min and a final extension at 72 ºC for 10 min (Hines et al., 2006). Efficacy of PCR reactions was checked in 1.5% agarose gel stained with Redsafe (Ecogen). Amplicons were sequenced using forward primers for EF-1α and ArgK fragments, and with both forward and reverse primers for the 16S fragment (Secugen, Madrid, Spain). Genotyping Microsatellites loci were amplified in two PCR multiplex reactions: RB1 (B10, B11, B100, B96, B124 and B126), and RB2 (B118, B119, B121, B131 and B132), with 0.2 μM concentration for each primer. Reactions also included MgCl2 (1.2 mM), dNTPs (0.3 mM), BSA (1.2 mg/ ml) and Kappa DNA polymerase (1.5 U, Kapa Biosystems®). PCR profiles for both multiplex were set following Cejas et al. (personal information) as follows: initial denaturation at 95 ºC for 5 min followed by 30 cycles of denaturation at 92 ºC for 30 s, annealing temperature at 54 ºC for 30 s, elongation at 72 ºC for 30 s and a final extension at 72 ºC for 30 min. PCR products were sent to the Servei Central de Suport a la Investigació Experimental (University of Valencia, Spain). Allele scoring was performed using GENEMAPPER 4.8 (Applied Biosystems Inc.) by comparing alleles with an internal size standard (GeneScan-500 Liz, Applied Biosystems Inc.). Binsets were built to correct genetic analyzer errors. Quality control was made manually. Data analysis DNA sequences were edited and aligned using Geneious 7.1.3 (http://www.geneious.com, Kearse et al., 2012). B. terrestris sequences of ArgK (AF492888.1), EF-1α (AF492955.1) (Kawakita et al., 2003) and 16S (AY737386.1) (Hines et al., 2006) genes were downloaded from Genbank (http://www.ncbi.nlm.nih.gov/genbank/) and included for comparison purposes. MEGA6 (Tamura et al., 2013) was used for sequence similarity analysis. Individuals with less than nine microsatellites amplified were discarded. Sibling workers from the same location inferred with Colony 2.0.6.4 (Wang, 2012) were excluded from further analyses. Colony parameters were selected as default for a haplodiploid species. Genetic parameters were calculated taken the three groups as three independent populations. GenAlex 6.5 (Peakall & Smouse, 2012) was used to obtain the number of alleles (Na), effective number of alleles (Ne), private number of alleles (Np) and observed (Ho) and expected (He) heterozygosity values. R program (R Core Team, 2013), with the adegenet package (Jombart, 2008; Jombart & Ahmed, 2011) was used to perform a discriminant analysis of principal components (DAPC), that is a multivariate method which allows probabilistic assignment of individuals to different clusters. Information about the subspecies (B. t. lusitanicus or B. t. terrestris), country, year of sampling and 16S haplotype of each individual in both samplings (preliminary and case study) was used to calculate (1) Kendall´s tau to estimate the correlation between subspecies and haplotype (hybrids were not taken into account, N = 147) and (2) Pearson´s chi-squared test with Yates continuity correction to analyze whether differences of 16S haplotype frequency in past and present samplings of B. t. lusitanicus were significant (N = 98). Statistical analyses were carried in R with package stats version 3.4.3 (R Core Team, 2013). Results Sequence analyses In total, 23 ArgK, 27 EF-1α and 30 16S sequences were obtained in the preliminary study. After trimming the sequences, the size of the ArgK fragment was 755 bp long. Seven point-mutations were detected (Table 2): five transitions (4 A↔G and 1 C→T) and two transversions (1 C→G and 1 C→A), but none of them were subspecies or geographic specific. The fragment amplified of the EF-1α gene was 522 bp long. Sequence alignment did not show any point mutations. The 16S trimmed fragment was 504 bp long. Two point-mutations were detected (Table 2): a transversion (A→T) and a transition (A→G). The transversion was only found in one B. t. terrestris individual (TTB.02), whereas the transition yielded two haplotypes that discriminate between the subspecies: individuals from the TLS group (Iberian B. t. lusitanicus) presented a guanine (haplotype-G) and those from TTB (Belgian B. t. terrestris) presented an adenine (haplotype-A). In the TTS group (B. t. terrestris from Spain), both haplotypes were found in five individuals each. Case study The final sequence set included 153 individuals (Table 1, Table S2 in Supplementary Material). In relation to B. t. lusitanicus 47 out of the 48 (97.92%) sampled in the 80s and 43 out of the 50 (86%) from the recent samplings presented the G-haplotype. The individual from the old sampling with G-haplotype was located at the north of Spain near the Pyrenees, whereas the individuals from the most recent surveys were distributed at the north (2), center (2) and south (3) of Spain (Fig S1 in Supplementary Material). All the French and Belgian B. t. terrestris individuals (100%) presented the A-haplotype; however, seven of the 18 B. t. terrestris individuals Sociobiology 65(4): 558-565 (October, 2018) Special Issue 561 (38.89%) recently sampled in Spain showed the A-haplotype, all of them from southern Spain (Sierra Nevada N. P.). Individuals considered morphologically as hybrids presented both haplotypes: two bear A-haplotype and four G-haplotype. Overall, 91.84% B. t. lusitanicus individuals showed G-haplotype, while 77.56% B. t. terrestris individuals showed A-haplotype. Kendall´s tau coefficient showed a significant correlation between subspecies and haplotype (tau = 0.705, p < 2.2e-16). However, when comparing French and Belgian B. t. terrestris populations with Spanish B. t. lusitanicus from the 80s, the correlation between both variables is almost one (tau = 0.973, p < 2.2e-16). Pearson´s chi-squared test for Iberian B. t. lusitanicus populations showed no significant differences (p = 0.074) between old and current individuals in haplotype segregation. Genotyping and clustering results One individual (TTB.05) was removed due to low amplification quality. Two possible siblings (TTS.01 and TTS.03, p = 1) were detected, therefore one sample (TTS.01) was randomly chosen and removed from posterior analyses (N = 28). Values of the number of alleles and heterozygosity (Table 3) were similar between TLS and TTB and slightly higher in TTS. The number of private alleles was different in the three groups, again higher in TTS. Bayesian Information Criterion (BIC) calculated prior to the DAPC yielded one as the most probable number of clusters. Even though, DAPC models with K = 2 and K = 3 were studied, but no correlation was found between the formed clusters with TLS, TTS and TTB groups (data not shown). Discussion Sequences of the two nuclear genes EF-1α and ArgK did not discriminate between the subspecies B. t. terrestris and B. t. lusitanicus in congruence with previous studies based on DNA barcoding (Williams et al., 2012) and the cephalic labial gland secretions (Coppée et al., 2008; Lecocq et al., 2016). Though some diversity was found along the ArgK sequence, it did not allow grouping the individuals into the two morphologically described subspecies. These genes were selected because they were previously used by Kawakita et al. (2003) and Hines et al. (2006) for analyzing Bombus intraspecificity. In fact, Hines et al. (2006) were able to distinguish between Bombus lucorum (Linnaeus) and Bombus Table 2. Point mutations found in the alignments of the two genes (ArgK and 16S). GenBank sequences from Bombus terrestris were concatenated and taken as reference (ArgK: AF492888.1, 16S: AY737386.1) (bp = position of the mutation in base pairs; . = same nucleotide as in the reference sequence; ― = missing data). ArgK 16S bp 69 269 393 558 705 706 749 161 263 reference G C G C A G C A A TLS.01 . . . . . . . G . TLS.02 A . A . . . . G . TLS.03 . . A . . . . G . TLS.04 ― ― ― ― ― ― ― G . TLS.05 A . . . . . . G . TLS.06 . . . . . . . G . TLS.07 . . . . . . . G . TLS.08 . . . . . . . G . TLS.09 A . A . . . . G . TLS.10 A . A . . . . G . TTS.01 . . . . . . . G . TTS.02 A . A . . . . . . TTS.03 A . A . . . . . . TTS.04 A . . . . . . G . TTS.05 ― ― ― ― ― ― ― . . TTS.06 A T A . . . . G . TTS.07 A . . . . . . . . TTS.08 . . . . . . . . . TTS.09 A . A . . . . G . TTS.10 A . A . . . . G . TTB.01 A . . . . . . . . TTB.02 ― ― ― ― ― ― ― . . TTB.03 A . . . . . G . . TTB.04 ― ― ― ― ― ― ― . . TTB.05 ― ― ― ― ― ― ― . . TTB.06 ― ― ― ― ― ― ― . . TTB.07 . . . A G A . . G TTB.09 A . . . . . . ― . TTB.10 A . A . . . . . . Group N Na Ne Np Ho He TLS 10 5.182±0.749 3.606±0.556 5 0.543±0.075 0.619±0.078 TTS 9 6.273±1.000 4.668±0.783 13 0.646±0.108 0.669±0.083 TTB 9 5.727±0.821 3.842±0.645 8 0.645±0.093 0.621±0.081 Table 3. Population genetic parameters across the three Bombus terrestris groups (TLS = B. t. lusitanicus from the Spain; TTS = B. t. terrestris from Spain; and TTB = B. t. terrestris from Belgium) based on eleven microsatellite loci (N = number of individuals, Na = number of alleles, Ne = number of effective alleles, Np = number of privative alleles, Ho = observed heterozygosity and He = expected heterozygosity). D Cejas; C Ornosa; I. Muñoz; P. De la Rúa – B. terrestris subspecies molecularly characterized in Spain562 hypnorum (Linnaeus) populations from Europe and China. In the same sense, the DAPC analysis based on microsatellite did not show a clustering agreement with the previously morphological defined groups. This low genetic differentiation between Bombus subspecies based on microsatellite variation has been reported before in mainland Europe (Estoup et al., 1996; Moreira et al., 2015), as a consequence of the high gene flow in B. terrestris populations, probably intensified due to contact and hybridization between commercial and wild B. terrestris populations. Interestingly, the mitochondrial gene sequenced here showed a variation which allows discriminating two haplotypes. The unique haplotype of the mitochondrial gene 16S retrieved in the Iberian individuals was not found in the Belgian and French B. t. terrestris individuals used as a reference. Furthermore, a significant correlation was found between haplotypes and subspecies in the case study, so it might be considered that G-haplotype identified B. t. lusitanicus and A-haplotype to B. t. terrestris. This is supported by the fact that all except one of the individuals collected between 1977 and 1985 before the establishment of the rearing companies and the start of the commercial use of B. terrestris (Velthuis & Doorn, 2006), showed the G-haplotype. The presence of individuals of B. t. terrestris in Spain may be due to two situations: on one hand, a natural contact area at the north of the Iberian Peninsula may have been established given the overlapping distribution range of both subspecies in the Pyrenees (Ornosa & Ortiz-Sánchez, 2004; Rasmont et al., 2008; Ornosa et al., 2017); on the other hand the presence of B. t. terrestris in the south of the Iberian Peninsula may be due to individuals that have escaped from greenhouses where they are widely used for pollination of crops such as tomatoes (Fig 1). The detection of hybrids and individuals with non-matching morphological subspecies and 16S haplotype is an evidence of subspecific introgression. The existence of individuals with B. t. terrestris phenotype and B. t. lusitanicus haplotype suggests that hybridization between commercially reared bumblebees and wild B. t. lusitanicus populations has occurred. The opposite case has also been observed: individuals with B. t. lusitanicus phenotype and B. t. terrestris haplotypes. Given the maternal inheritance of the mitochondrial molecule we can deduce that both commercial queens and males have colonized the environment, which is still a discussed topic (Velthuis & Doorn, 2006; Kraus et al., 2010, Lecocq et al., 2015), meaning that B. t. terrestris is naturalized, for now, in the south of Spain. Fig 1. Frequency of Bombus terrestris lusitanicus, Bombus terrestris terrestris and hybrids (those with morphological subspecies identification and 16S haplotype that did not match i. e. B. t. lusitanicus with A-haplotype and B. t. terrestris with G-haplotype) in the National Park of Sierra Nevada (southern Spain). Size of the circles is proportional to the total number of individuals sampled at each location. The high density of greenhouses can be seen in the lower right corner of the figure. Ortophotograph obtained from http://www.juntadeandalucia.es. Sociobiology 65(4): 558-565 (October, 2018) Special Issue 563 B. terrestris subspecies usually do not interbreed due to the specificity of the hormones secreted by the male cephalic labial glands (Coppée et al., 2008). Yet, in the case of B. t. terrestris and B. t. lusitanicus, the secretions are chemically very similar (Coppée et al., 2008) allowing more usual copulation between them than with other B. terrestris subspecies. In other attempts of colonizing new lands by B. terrestris (Ings et al., 2010), the subspecies B. t. sassaricus stopped to be seen in the south of France two years after the end of its importation, and no hybrids were found although no molecular analyses were performed. In the case study, naturalized individuals of B. t. terrestris on the Iberian Peninsula may have colonized the environment in the south where the breeding companies are based and, given their ability to breed with B. t. lusitanicus, the number of hybrid individuals in this area may increase in the future. The effects that this hybridization is having on wild populations are unknown (Facon et al., 2011), as are the effects that these introduced populations may have on native flora (Aizen et al., 2018). However, according to the expansion models of Lecocq et al. (2015), B. t. terrestris should not be able to continue its colonization on the Iberian Peninsula, as it is not suitable for its climatic conditions. Future studies including more localities within the Iberian Peninsula and using preferably the 16S marker and other highly variable nuclear markers such as SNPs (Single Nucleotide Polymorphisms) will show the degree of introduction of commercial B. terrestris over the territory, what is of conservation concern due to a possible impact on the genetic diversity of Iberian populations or even to the displacement of the native wild populations. Acknowledgments We would like to thank Kevin Maebe and Guy Smagghe for providing samples needed for this study, Laura Jara, Vicente Martínez and Ana Isabel Asensio for helping with the analyses and Carlos Ruiz for productive discussion and two anonymous reviewers for comments on a previous version. Sampling permissions were obtained from the corresponding authorities (National Parks of Sierra Nevada and Guadarrama). This work has been supported by the Ministry of Economy and Competitiveness (BIOBOMBUS project CGL2012-34897), INIA-FEDER (RTA2013-00042- C10-05) and Fundación Séneca (Plan of Science and Technology of the Region of Murcia, grant of Regional Excellence 19908/GERM/2015). Irene Muñoz is supported by Fundación Séneca (Murcia, Spain) through the post-doctoral fellowship “Saavedra Fajardo” (20036/SF/16). Supplementary Material DOI: 10.13102/sociobiology.v65i4.3442.s2222 Link: http://periodicos.uefs.br/index.php/sociobiology/rt/supp Files/3442/0 Authors’ Contribution PDlR and CO designed the experiments. CO organized the sampling and identified the taxonomic status of individuals. DC made the experimental work. Data analysis and results interpretation were carried on by DC and IM. First draft of this paper was redacted by DC, and final draft was revised and approved by PDlR, IM and CO. References Acosta, A.L., Giannini, T.C., Imperatriz-Fonseca, V.L. & Saraiva, A.M. (2016). Worldwide Alien Invasion: A Methodological Approach to Forecast the Potential Spread of a Highly Invasive Pollinator. PLoS ONE, 11: e0148295. doi: 10.1371/journal.pone.0148295 Aizen, M.A., Morales, C.L., Vázquez, D.P., Garibaldi, L.A., Sáez, A. & Harder, L.D. (2014). When mutualism goes bad: density- dependent impacts of introduced bees on plant reproduction. New Phytologist, 204: 322-328. doi: 10.1111/nph.12924 Aizen, M.A., Smith-Ramírez, C., Morales, C.L., Vieli, L., Sáez, A., Barahona-Segovia, R.M., Arbetman, M.P., Montalva, J., Garibaldi, L.A., Inouye, D.W. & Harder, L. (2018). Coordinated species importation policies are needed to reduce serious invasions globally: The case of alien bumblebees in South America. Journal of Applied Ecology, doi: 10.1111/1365-2664.13121 Cameron, S., Hines, H.M. & Williams, P.H. (2007). A comprehensive phylogeny of the bumble bees (Bombus). Biological Journal of the Linnean Society, 91: 161-188. doi: 10.1111/j.1095-8312.2007.00784.x Coppée, A., Terzo, M., Valterova, I. & Rasmont, P. (2008). Intraspecific Variation of the Cephalic Labial Gland Secretions in Bombus terrestris (L.) (Hymenoptera: Apidae). Chemistry & Biodiversity, 5: 2654-2661. doi: 10.1002/cbdv.200890219. Dreier, S., Redhead, J., Warren, I., Bourke, A., Heard, M., Jordan, W., . . . Carvell, C. (2014). Fine-scale spatial genetic structure of common and declining bumble bees across an agricultural landscape. Molecular Ecology, 23: 3384-3395. doi: 10.1111/mec.12823 Estoup, A., Solignac, M., Harry, M. & Cornuet, J. (1993). Characterization of (GT)n and (CT)n microsatellites in two insect species: Apis mellifera and Bombus terrestris. Nucleic Acids Research, 21: 1427-1431. Estoup, A., Solignac, M., Cornuet, J., Goudet, J. & Scholl, A. (1996). Genetic differentiation of continental and island populations of Bombus terrestris (Hymenoptera: Apidae) in Europe. Molecular Ecology, 5: 19-31. Facon, B., Crespin, L., Loiseau, A., Lombaert, E., Magro, A. & Estoup, A. (2011). Can things get worse when an invasive species hybridizes? The harlequin ladybird Harmonia axyridis D Cejas; C Ornosa; I. Muñoz; P. De la Rúa – B. terrestris subspecies molecularly characterized in Spain564 in France as a case study. Evolutionary Applications, 4: 71- 88. doi: 10.1111/j.1752-4571.2010.00134.x. Goulson, D., Nicholls, E., Botías, C. & Rotheray, E. (2015). Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science, 347: 1255957. doi:10.1126/science.1255957 Hines, H., Cameron, S. & Williams, P. (2006). Molecular phylogeny of the bumble bee subgenus Pyrobombus (Hymenoptera: Apidae: Bombus) with insights into gen utility for lower-level analysis. Invertebrate Systematics, 20: 289- 303. doi: 10.1071/IS05028 Ings, T.C., Ings, N.L., Chittka, L. & Rasmont, P. (2010). A failed invasion? Commercially introduced pollinators in Southern France. Apidologie, 41: 1-13. doi: 10.1051/ apido/2009044 Inoue, M., Yokoyama, J. & Washitani, I. (2008). Displacement of Japanese native bumblebees by the recently introduced Bombus terrestris (L.) (Hymenoptera: Apidae). Journal of Insect Conservation, 12: 135-146. doi: 10.1007/s10841-007-9071-z Ivanova, N., Dewaard, J. & Herbert, D. (2006). An inexpensive, automation-friendly protocol for recovering high-quality DNA. Molecular Ecology Notes, 6: 998-1002. doi: 10.1111/j.1471-8286.2006.01428.x Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics, 24: 1403-1405. doi:10.1093/bioinformatics/btn129 Jombart, T. & Ahmed, I. (2011). adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics, 24: 1403–1405. doi: 10.1093/bioinformatics/btr521 Kawakita, A., Sota, T., Ascher, J., Ito, M., Tanaka, H. & Kato, M. (2003). Evolution and Phylogenetic Utility of Alignment Gaps Within Intron Sequences of Three Nuclear Genes in Bumble Bees (Bombus). Molecular Biology and Evolution, 20: 87-92. doi: 10.1093/molbev/msg007 Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., … Drummond, A. (2012). Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics, 28: 1647-1649. doi: 10.1093/bioinformatics/bts199 Kraus, F., Szentgyörgyi, H., Rozej, E., Rhode, M., Moron, D. & Woyciechowski, M. (2010). Greenhouse bumblebees (Bombus terrestris) spread their genes into the wild. Conservation Genetics, 12: 187-192. doi:10.1007/s10592-010-0131-7 Lautenbach, S., Seppelt, R., Liebscher, J. & Dormann, F.C. (2012). Spatial and Temporal Trends of Global Pollination Benefit. PloS One, 7: e35954. doi: 10.1371/journal.pone.0035954 Lecocq, T., Rasmont, P., Harpke, A. & Schweiger, O. (2015). Improving International Trade Regulation by Considering Intraspecific Variation for Invasion Risk Assessment of Commercially Traded Species: The Bombus terrestris Case. Conservation Letters, 9: 281-289. doi:10.1111/conl.12215 Lecocq, T., Coppée, A., Michez, D., Brasero, N., Rasplus, J.Y., Valterová, I. & Rasmont, P. (2016). The alien’s identity: consequences of taxonomic status for the international bumblebee trade regulations. Biological Conservation, 195: 169-176. doi: 10.1016/j.biocon.2016.01.004 Matsumura, C., Yokoyama, Y., & Whasitani, I. (2004). Invasion Status and Potential Ecological Impacts of an Invasive Alien Bumblebee, Bombus terrestris L. (Hymenoptera: Apidae) Naturalized in Southern Hokkaido, Japan. Global Environmental Research, 8: 51-66. Morales, C.L., Arbetman, M.P., Cameron S.A. & Aizen, M.A. (2013). Rapid ecological replacement of a native bumble bee by invasive species. Frontiers in Ecology and the Environment, 11: 529-534. doi: 10.1890/120321 Moreira, A., Horgan, F., Murray, T. & Kakouli-Duarte, K. (2015). Population genetic structure of Bombus terrestris in Europe: Isolation and genetic differentiation of Irish and British populations. Molecular Ecology, 24: 3257-3268. doi: 10.1111/mec.13235 Murray, T.E., Fitzpatrick, U., Brown, M.J.F. & Paxton, R.J. (2008). Cryptic species diversity in a widespread bumble bee complex revealed using mitochondrial DNA RFLPs. Conservation Genetics, 9: 653-666. doi: 10.1007/s10592- 007-9394-z Ornosa, C. (1996). Una nota de atención sobre la introducción artificial de subespecies foráneas de abejorros polinizadores en la Península Ibérica (Hymenoptera, Apidae, Bombinae). Boletín de la Asociación Española de Entomología, 20: 259-260. Ornosa, C. & Ortiz-Sánchez, F. (2004). Hymenoptera: Apoidea I, Fauna Ibérica (Vol. 23). Madrid: Museo de Ciencias Naturales, CSIC, 553 p Ornosa, C., Torres, F. & De la Rúa, P. (2017). Updated list of bumblebees (Hymenoptera: Apidae) from the Spanish Pyrenees with notes on their decline and conservation status. Zootaxa, 4237: 41-77. doi: 10.11646/zootaxa.4237.1.3. Peakall, R. & Smouse, P. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics, 28: 2537-2539. doi: 10.1093/bioinformatics/bts460 R Development Core Team (2013). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Rasmont, P., Coppée, A., Michez, D., & De Meulemeester, T. (2008). An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae). Annales de la Société Entomologique de France (N.S.), 44: 243-250. doi: 10.1080/ 00379271.2008.10697559 Sociobiology 65(4): 558-565 (October, 2018) Special Issue 565 Rokas, A., Nylander, J., Ronquist, F. & Stone, G. (2002). A maximum-likelihood analysis of eight phylogenetic markers in gallwasps (Hymenoptera: Cynipidae): implications for insect phylogenetic studies. Molecular Phylogenetics and Evolution, 22: 206-219. doi: 10.1006/mpev.2001.1032 Schmid-Hempel, R., Eckhardt, M., Goulson, D., Heinzmann, D., Lange, C., Plischuk, S., . . . Schmid-Hempel, P. (2014). The invasion of southern South America by imported bumblebees and associated parasites. Journal of Animal Ecology, 83: 823- 837. doi: 10.1111/1365-2656.12185 Sutherland, W. J., Barnard, P., Broad, S., Clout, M., Connor, B., Côté, I. M., … Ockendon, N. (2016). A 2017 Horizon scan of emerging issues for global conservation and biological diversity. Trends in Ecology and Evolution, 32: 31-42. doi: 10.1016/j.tree.2016.11.005 Tamura, K., Stecher, G., Peterson, D., Filipski, A. & Kumar, S. (2013). MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology and Evolution, 30: 2725-2729. doi:10.1093/molbev/mst197 Vargas, P., Ornosa, C., Blanco-Pastor, J.L., Romero, D., Fernández-Mazuecos, M. & Rodríguez-Gironés, M.A. (2013). Searching for areas of genetic diversity in Sierra Nevada: analyses of plants and bees. In L. Ramírez & B. Asensio (Eds.), Proyectos de investigación en parques nacionales: 2009-2012 (pp. 123-142). Naturaleza y Parques Naturales. Velthuis, H.H..W & Doorn, A. (2006). A century of advances in bumblebee domestication and the economic and environmental aspects of its commercialization for pollination. Apidologie, 37: 421-451. doi:10.1051/apido:2006019 Wang, J. (2012). Computationally efficent sibship and parentage assignment from multilocus marker data. Genetics, 191: 183-194. doi: 10.1534/genetics.111.138149 Whitehorn, P., Tinsley, M., Brown, M. & Goulson, D. (2013). Investigating the impact of deploying commercial Bombus terrestris for crop pollination on pathogen dynamics in wild bumble bees. Journal of Apicultural Research, 52: 149-157. doi: 10.3896/IBRA.1.52.3.06 Williams, P. H., Brown, M. J., Carolan, J. C., An, J., Goulson, D., Aytekin, A. M., ... & Huang, J. (2012). Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Systematics and Biodiversity, 10: 21-56. doi: 10.1080/14772000.2012.664574 _Hlk513543699