GENETIC RELATIONSHIPS BETWEEN COMMERCIALLY PRODUCED AND NATURAL POPULATIONS OF BOMBUS TERRESTRIS DALMATINUS IN TERMS OF MITOCHONDRIAL COI AND CYTB

Bombus terrestris dalmatinus is naturally common in many countries, including Turkey, and is also used commercially for the pollination of greenhouse plants. Intensive commercial production and international trade in many countries are considered as reasons for the disappearance of some natural populations. Hybridization of native bumble bees with those produced commercially, but having escaped from greenhouses and colonization of these commercial bees in natural habitats are cause for concern. In order to assess this concern, B. t. dalmatinus workers were collected from twelve different populations: five commercial producers, three surrounding greenhouse centers, three natural areas at least 30 km away from greenhouses, and one more recent greenhouse zone in Antalya, Turkey. The genetic variations and relationships among the twelve populations were estimated using SNP haplotypes determined in mitochondrial COI and CytB. Twenty and sixteen haplotypes were obtained for COI and CytB, respectively. A single haplotype, H1, was widespread with a high frequency in all individuals for both genes. Individuals collected from around greenhouse centers and commercial companies had more common haplotypes. The genetic variations of intra-populations were higher than the inter-populations in both COI (65.41%>34.59%) and CytB (72.47%>27.53%). The natural and commercial populations were genetically more distant from each other considering F st values. However, samples from near the greenhouses had a higher similarity with the commercially produced samples, while the natural populations far away from greenhouses still retained their genetic distinctiveness.


INTRODUCTION
Bombus species are adapted to different climatic and environmental conditions, having a wide range of habitats from sea level up to 5800 m. Bombus terrestris (Linnaeus, 1758) shows abilities for high migration, early colony development, to benefit from a large number of flowers and to withstand low temperatures (Goka et al., 2001;Goulson et al., 2008). This species is one of the most important pollinators of both natural and cultivated plants (Michener, 2000) and constitutes about 90% of commercially produced colonies (Velthuis & van Doorn, 2006). It is now used for pollination in over fifty-seven countries including Turkey and has not been found naturally in sixteen of them. Bombus terrestris dalmatinus Dalla Torre is one of the most suitable subspecies for commercial breeding because of its high adaptability and reproductive performance under captive conditions (Velthuis & van Doorn, 2006). The many benefits of commercially produced Bombus colonies including increased yield and quality of fruit in undergrowth cultivation and reduced use of hormones and chemical drugs have increased demand globally since the 1980s (Ono, 1998;Gürel et al., 1999;Goulson & Hanley, 2004;Hingston, 2006;Velthuis & van Doorn, 2006;Inoue et al., 2008;Williams & Osborne, Genetic relationships among Bumblebee populations 2009). However, bees emerging from commercial colonies are reported to cause problems (Goulson & Hughes, 2015) such as competition with local Bombus species for nesting and ecotypes (Goulson, 2003;Inoue & Yokoyama, 2010;Aizen, 2018) because it has escaped from greenhouses (Seabra et al., 2019) and colonized the environment (Kraus et al., 2011). Trillo et al. (2019) reported the presence of managed bumblebees escaping from greenhouses and feeding on the same flowering plants with the native ones, especially in winter. Thus, there has likely been genetic contamination due to hybridization (Rhymer & Simberloff, 1996;Ings et al., 2006;Kanbe et al., 2008;Goulson, 2010). In addition, this species, which is highly invasive (Hingston et al., 2002), is thought to cause the spread of exotic diseases (Whitehorn et al., 2013). Due to such reasons, natural bumblebee populations are thought to be negatively affected and even at risk of extinction in some places (Goulson, 2003;Cejas et al., 2018;Tsuchida et al., 2019). However, there is a lack of knowledge about the current state of most populations (Chandler et al., 2019). Genetic studies on the detection of introgression between commercial and native Bombus populations are still scarce (Kraus et al., 2011;Seabra et al., 2019;Cejas et al., 2020). Each year, approximately two million commercially produced Bombus colonies are estimated to be shipped globally (Lecocq et al., 2016), and almost 300,000 colonies are used per year, especially in tomato pollination, in Turkey alone (Cilavdaroglu & Gurel, 2020). In greenhouses in Turkey there are no measures in place to prevent the escape of bumblebees into the natural habitat. Furthermore, as is also usually practiced in Portugal, hive boxes at the end of their useful life are left outside of greenhouses with bees still in residence (Seabra et al., 2019). In many species, mitochondrial DNA (mtDNA) is used for identifying the inter and intra specific differences due to high evolution ratio compared to nuclear DNA and to non-recombining and maternal inheritance (Rubinoff & Holland, 2005). In this context, animal mtDNA is a reliable and valid marker for population genetics and phylogenetic studies. Although the B. terrestris subspecies differ morphologically, reliable molecular tools are needed because the features are insufficient in some subtype and ecotype level definitions. mtDNA is widely used in phylogenetic studies, in particular, COI and CytB sequencing data accumulated for bumblebee species are used to reveal genetic relationships between bumblebee populations (Pedersen, 1996;Widmer & Schmid-Hempel, 1999;Yoon et al., 2003;Murray et al., 2008;Kim, 2009;Tokoro, 2010;Williams et al., 2012;Han et al., 2019). The genetic structure of commercial colonies taken from companies in Antalya and their possible genetic effects on local genotypes are not fully known. Therefore, the genetic diversity of B. t. dalmatinus populations in habitats may be at risk of being lost without being identified. The information obtained about the genetic basis and the population structure of the geographical variation will serve as an essential resource for the use and safeguarding of bumblebees, but for now molecular genetic studies on bumblebees are limited in Turkey. In order to assess this concern, this study was carried with the use of COI and CytB haplotypes to evaluate whether there has been a hybridization between commercially reared and indigenous populations of B. t. dalmatinus in Antalya.

Sample Collection
Fifteen B. t. dalmatinus workers from each of the twelve different locations (five commercial producers, seven collected from natural fields including three nearby greenhouse centers and four natural districts) were sampled from the Antalya province. Three foreign and two domestic commercial companies engaged in B. t. dalmatinus breeding. Three of seven different fields are greenhouse centers (Aksu, Demre, and Kumluca) where intensive greenhouse activities are carried out, and a large number of commercially produced B. t. dalmatinus colonies are used. Bees easily escape from these greenhouses because there are no preventive measures in general. Bumblebees were collected at a distance of about 200 meters from the greenhouses in these centers. The other three of the seven field collections are natural districts (Bayatbademler, Phaselis, and Termessos) which are at least 30 km away from the nearest greenhouse centers and isolated with natural barriers. Bayatbademler is a plain surrounded by mountains, Phaselis an ancient city at sea level and Termessos a national park top of the mountain. The last field site, Geyikbayiri, is a plateau where some greenhouse activities have started in recent years. Also, an outgroup (Lecocq et al., 2013) was used as a reference to determine the phylogenetic relationships of the populations.
The bumblebees were caught with use of an entomological net. The coordinates and altitude information of the sampling locations are given in Fig. 1. The samples were stored in collection tubes with pure ethanol at +4°C until DNA extraction.

DNA
extraction, amplification, and sequencing The thorax of each bumblebee was placed in an individual tube and crushed through the application of cold nitrogen. The CTAB method was used to extract DNA (Doyle, 1990). The quantities and qualities of DNA were determined through BioDrop spectrophotometer. In addition, DNA

Processing of data and statistical analyses
The COI and CytB gene region sequences of mtDNA were examined for each individual with the use of ChromasPro version 1.7.4 (Technelysium Pty. Ltd. Australia), and incorrect and uncertain reads resulting from sequencer were edited based on visual examination of the chromatograms. The ends of the reads presenting lower quality bases were trimmed and thus the length of the sequence used in further analysis has been shortened. Analyzes were continued with DNA sequences from 157 individuals for COI (204 bp:base pairs) and 170 individuals for CytB (373 bp). At several base positions there were double peaks in the chromatograms, and examples are shown in the Sup. Fig. 1. Since all suspicious artifacts were removed from the sequences, the sequences with double peaks were considered heteroplasmic. Thus, sequences belonging to all individuals were used twice, whether heteroplasmic or not, in order not to make a propor-  tional error. So, including these sites, haplotype analysis was performed on total of 314 (2*157) DNA sequences for the COI and 340 (2*170) DNA sequences for the CytB. Sequenced DNA regions were searched against BLAST (GenBank: JQ820651. 1 for COI and GenBank: JQ820853. 1 for CytB). Sequence alignment was conducted with MEGA 6 (Molecular Evolutionary Genetics Analysis, version 6.06.) software (Tamura et al., 2013). SNPs were identified, and the haplotypes for each region were generated with DnaSP (version 5. 10.01) software (Rozas, 2010). Also, such diversity parameters as haplotype (gene) diversity (Hd) and nucleotide diversity (π) were calculated for the studied gene regions. Arlequin (version 33. 11) software was used to calculate genetic differentiation within and among the populations (Excoffier et al., 2007) with the unbiased fixation index (F ST ) based on information of haplotypes (p<0.01, p<0.05). Molecular variation (AMOVA) was analyzed to calculate genetic variation values for both studied gene regions through the use haplotype sequences detected in COI and CytB. Phylogenetic analyses were made with TASSEL (version 5.2.43) software (Bradbury et al., 2007) through the use of the genetic distance matrices among populations. Mutational relationships among haplotypes were represented by median-joining (MJ) Network (version 5.0.0.3) software (Bandelt et al., 1999).

Sequence analysis and haplotype distributions
In all analyzed individuals, the COI fragment sequence had 204 base pairs, and the CytB fragment sequence had 373 base pairs. According to the results of the BLAST analyses, the DNA nucleotide sequences obtained from mtDNA COI and CytB gene fragments were determined to belong to B. t. dalmatinus and   Other haplotypes were restricted to one, two or five populations and had a low frequency (Tab. 2 and Tab. 3). For example, H3, H4, H6, H9, and H10 in terms of COI and H2, H7, H8, H9, and H19 in terms of CytB were only in samples taken from companies (Tab. 2 and Tab. 3). After haplotype (gene) and nucleotide diversity parameters were calculated, the highest values were found in natural areas for both genes (Tab. 4 and Tab. 5). The distribution of haplotypes was visualized with the Median-Joining Network ( Fig. 2a-b), but both network patterns showed that there was some differentiation between groups in general ( Fig. 2a-b). These figures (2a-b) show that greenhouse areas and commercial companies had more common haplotypes, and conversely natural areas were visibly separated from the others. Bees collected from greenhouse areas were closer to commercially produced bees genetically according to network pattern.

Genetic relationships
As a result of AMOVA (Tab. 4-5), the percentage of variation within populations was higher than among populations for both COI and CytB gene regions in all populations and each group. The mean F ST values were calculated as 0.346 for COI and 0.275 for CytB. Also, the highest and the lowest pairwise F ST values for COI were detected as 0.78 between Aksu and Phaselis and as 0.00 between Company_2 and Termessos, respectively (Tab. 6). The highest and the lowest F ST values for CytB were detected as 0.50 between Company 4 and Termessos and as 0.02 between Geyikbayiri and Phaselis, respectively (Tab. 6). Genetic distances between populations are lower among individuals who are commercially produced and caught around the greenhouse regions.
In order to determine the genetic relationships among the twelve different populations of B. t. dalmatinus, a phylogenetic-tree were constructed based on the Neighbor-Joining (Saitou & Nei, 1987) method with the use of TASSEL software v. 5.2.59. Bees from twelve different sources were clearly distinguished according to haplotypes (Sup. Fig. 2a-b). The closer relationship of commercially produced bumblebees with the ones caught in the greenhouse region suggests that they are in a closer genetic relationship than the bees caught from natural areas.    Seabra et al. (2019) supports this result reporting a relatively common haplotype in commercial individuals and in wild individuals collected nearby the greenhouses where the commercial hives are used. The presence of jointly shared haplotypes indicates an affinity between commercially produced individuals and individuals collected from the greenhouse environment. This result suggests that individuals can go leave greenhouses and survive around greenhouses in the natural environment. Candidate hybrids were detected in the wild, as well as putatively escaped commercial bumblebees, some being potentially fertile males (Seabra et al., 2019). In another study, commercial individuals escaping from greenhouses were shown to hybridize with individuals in the vicinity. In addition, if natural bees are far enough from greenhouse areas, they show the ability to preserve their original genetic structure (Kraus et al., 2011). Cejas et al. (2018) showed 16s haplotype as evidence of integration between B. terrestris subspecies. As seen in Tables 2 and 3, similar haplotype numbers (COI=16 and CytB=20) were determined in this present study. The similarity of the number of haplotypes supports that we can use it as evidence to demonstrate integration among sub-populations in the Antalya region. Pedersen (1996) identified the nucleotide sequence of the fragment of 532 bp in the mitochondrial COI gene of bumblebees and revealed the phylogenetic relationships of eleven Bombus species. Tokoro et al. (2010) studied the mitochondrial COI (1048 bp) gene and identified fifteen haplotypes in order to determine the risk of the genetic breakdown in local Japanese bumblebees and the outsourced gene flow of samples collected from Japan, China and Korea. When the literature is examined, the number of haplotypes calculated in this study is notably higher than those found in other studies. This situation is due to Turkey's very different geographical and climatic conditions, the intersection of Europe, Asia and Africa in Turkey and the high genetic diversity of Turkey's B. terrestris populations. AMOVA analysis revealed that genetic variation within populations is greater than among populations. One reason for this situation may be hybridization between populations. Individuals carrying a genomic structure from one population to another may have increased intrapopulation variation due to hybridization, while reducing genetic variation between populations. The B. t. dalmatinus species in Antalya may be genomically a large population and the differences between sub-populations are still low. Genetic similarities and differences were calculated through the pairwise F ST values for mtDNA COI and CytB gene fragments of twelve populations. According to these results, individuals gathered from around the greenhouses were genetically closer to commercially produced individuals than natural ones. When Tab. 4 and Tab. 5, which show the results of AMOVA, were examined, the remarkable results were that the intra-genetic variations among the commercial populations and among the populations from the near the greenhouses were lower than the natural populations far from greenhouse areas. Understandably commercial populations produced thorugh culturing and the populations from the near the greenhouses showed a more homogeneous genetic structure. Moreira et al. (2015) investigated the genetic structure in twenty-two natural and two commercial B. terrestris L. populations using eight microsatellites and two mitochondrial genes (COI and CytB). The barrier formed by the Irish Sea and the dominant south winds is thought to prevent gene flow between Ireland and Britain, and also heterozygosity in Ireland and Isle of Man was reported to be lower than in the European continent and commercial populations. The data suggest that B. terrestris individuals in the west of Ireland, where the use of managed bumblebees is rare, associated with populations from the European mainland in terms of COI. Estoup et al. (1996) stated that B. terrestris populations in the European continent have a relatively homogeneous genetic structure due to the low overall geographical barriers. Based on this idea, we think that individuals who do not encounter any barrier in greenhouses without precautions easily affect the genetics of natural bees. Although few studies have been conducted in Turkey, up to fifty different Bombus species have been identified. Based on world-species distribution, Turkey is understood to be an important genetic resource in terms of bumblebees (Özbek, 1997;Barkan & Aytekin, 2013;Meydan et al., 2016), but it is not clear whether there is a hybridization between natural populations and commercial ones. Thus, the genetic diversity of Turkey's local B. t. dalmatinus populations must be determined and utilized in protection strategies as the source of genes. The data revealed by this completed study will be useful for this purpose. The Antalya province is where B. t. dalmatinus is widely found in natural habitats, and many sub-species are also encountered (Barkan & Aytekin, 2013), which is where very intensive greenhouse activities are conducted. Intensive use of commercial B. t. dalmatinus colonies are thought to increase the numbers of escaped bumblebees from greenhouses, but the spread of a large number of commercially produced individuals leads to increased competition with natural populations for nesting and feeding, spread of diseases and uncontrolled crossbreeding with local genotypes (Dafni, 1998;Goka et al., 2001;Whitehorn et al., 2013;Aizen, 2018). Furthermore, the introduction of foreign Bombus terrestris (Linnaeus) has resulted in a decline in native bumblebee populations in such countries as Japan, Chile and Argentina (Cejas et al., 2018).
Although it is difficult to prove that commercially produced bees and wild bees are hybridized, bees fleeing from greenhouses can compete with natural bees in terms of nest and nutrient sharing. The invasive potential of commercially produced bumblebees has been clearly identified in several studies. Ings et al. (2006) used a paired design to compare the nectar-foraging performance and reproductive outputs of commercial and native colonies growing under identical field conditions, and they found that the commercial colonies have high reproductive success, superior foraging ability and large colony size. In a study conducted in Turkey, Gösterit (2017) determined that the colonies founded by commercial queens produced more than twice the number of gynes (82. 11±9.32) than colonies founded by native queens (32.85±3.99). Pirounakis et al. (1998) reported that the CytB gene could be successfully used to identify geographical subspecies of bumblebees. In another study using the CytB gene, the genetic structure of B. terrestris bees on two islands in the African Gulf was found to be heterogeneous (Widmer et al., 1998). Koulianos and Hempel (2000) investigated the genetic differences and kinship relationships using mitochondrial CytB and COI genes in nineteen bombus species collected from sixteen regions in Europe and three areas in North and South America. Morath (2007) used the mtDNA CytB gene to determine genetic variation in Bombus impatiens. In addition, the genetic structure of seven mainland and island Asian populations of Bombus ignitus was investigated with the use of nine microsatellite markers, and the sequences of part of the mitochondrial CytB gene were determined. The results of this study, carried out with the samples collected from most of the companies in Antalya and the areas where the bumblebees are concentrated, overlap with many of the previous studies and summarize the genetic relationships among the populations in the Antalya region. Although the investigation of heteroplasmy was Genetic relationships among Bumblebee populations not the subject of this study we surprisingly encountered individuals with double peaks, which were too obvious to be considered as artifacts in sequencing of both directions. (Sup. Fig. 1). As mentioned in previous similar studies (Wernick et al., 2016;Williams et al., 2019;Pizzirani et al., 2020;Ricardo et al., 2020;Tikochinski et al., 2020) double peaks on sequence traces with both alleles from one individual were identified as mitochondrial heteroplasmy and a single individual with two haplotypes. At least one heteroplasmy was detected in 93 of 157 individuals for COI and at least 105 of 170 individuals for CytB in the current study. Ricardo et al. (2020) reported that heteroplasmy was detected in individuals from all the ten sampled locations with an average of six heteroplasmic haplotypes per individual. Moreover, they found that some of these heteroplasmic haplotypes are shared between individuals from different locations. This situation regarding heteroplasmy in this study suggested that new studies are needed to prove it. Similarly, after the initial definition of heteroplasmy for Bombus morio (Francoso et al., 2016) further investigation was necessary (Ricardo et al., 2020). The bee samples from near the greenhouses have been determined to have high similarity with the commercially produced individuals. These results suggest that commercial colonies escaped from greenhouses and can cause genetic hybridization with local populations around greenhouses. However, the natural populations collected from different native habitats far away from greenhouses constitute different branches of phylogenetic tree from the commercial. This shows that they can still retain their genetic constructs, but uncontrolled and imprudent use of commercial colonies may soon be a serious threat to the original genetic stock of B. t. dalmatinus. In such areas as Antalya where hundreds of thousands of commercial colonies are produced every year and used in greenhouses, very serious measures must be taken.

Disclosure statement
The authors declare that they have no conflict of interest.