Searching for Molecular Markers to Differentiate Bombus terrestris (Linnaeus) Subspecies in the Iberian Peninsula

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


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 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 nonnative 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 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. (1993Estoup et al. ( , 1996) ) adequate to analyze the genetic structure of bumblebee populations (Kraus et al., 2010, Dreier et al., 2014, Moreira et al., 2015).

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)

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).
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).

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 (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 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.Kraus et al., 2010, Lecocq et al., 2015), meaning that B. t. terrestris is naturalized, for now, in the south of Spain.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.

Fig 1 .
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.

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.

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).
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;