HRM Analysis of Spermathecal Contents to Determine the Origin of Drones that Inseminated Honey Bee Queens

Western honey bee (Apis mellifera L.) has a wide biodiversity and many of its subspecies have adapted to the habitats of Africa, Asia and Europe. The subspecies have differentiated in terms of morphological, molecular, behavioral and physiological characteristics, and to date approximately thirty western honey bee subspecies have been described (Meixner et al., 2011; Ruttner, 1988; Sheppard et al., 1997; Sheppard & Meixner, 2003; Chen et al., 2016). Analyses showed that the morphological characteristics of western honey bee subspecies were distributed into five morphogenetic lineages: Tropical Africa (A), North and Western Europe (M), South-Eastern Europe (C), Near East (O) and Yemenitica (Y) (Kauhausen-Keller, Ruttner, & Keller, 1997; Ruttner, Tassencourt, & Louveaux, 1978; Ruttner, 1988).

Four western honey bee subspecies, Carniolan bee (A. m. carnica), European Dark bee (A. m. mellifera), Italian bee (A. m. ligustica) and Caucasian bee (A. m. caucasica), are of current economic importance, and the most common subspecies breeding/produced all over the world. Since honey bees adapt to the particular geographic area and conditions in which they live, each subpopulation has potential value with regard to their unique genetic and phenotypic traits, and it is important to protect them (Rúa et al., 2013).

Such basic reproductive characteristics of honey bees as haplodiploid sex determination mechanism, mating of the queens with more than one drone (polyandry) at one or several nuptial flights, drones dying immediately after their single mating (monogamy) and storing live sperm for fertilization in spermatheca during their life (Winston, 1987; Baer, 2005) make it difficult to conserve genetically distinct subpopulations, but the most important are polyandry and in-flight mating at unpredictable drone congregation areas (Baudry et al., 1998; Koeniger et al., 2005). Honey bee queens mate with between seven and eighteen drones (average ten drones) during one or several nuptial flights (Woyke, 1962). Estimated mating frequencies in various studies range from 6.5 (Taber, 1954) to 41.3 (Kraus et al., 2004).

The conservation of honey bee genetic sources are possible through the controlled mating within closed populations or through instrumental insemination to prevent introgression between subpopulations. Besides the main challenges of honey bee reproductive biology, anthropogenic and beekeeping activities including migratory beekeeping, use of foreign queens and selective breeding also affect genetic resources by increasing genetic diversity or mixing honey bee populations (Harpur et al., 2012, 2013; Rúa et al., 2009 and 2013). The combination of these factors increases the degeneration (or speed of the genetic admixture) of endemic honey bee subspecies and ecotypes with foreign subspecies. The loss of valuable traits along with the degeneration of native honeybee populations poses a threat to sustainable beekeeping and agriculture. Therefore, successful conservation strategies are under development to prevent introgression between native and foreign honey bee populations (Rúa et al., 2009).

Practical molecular genetic tools should be developed to mating control and determine whether there is gene flow from foreign sources in honey bee populations (Zayed, 2009). In this study, a simple and fast method is suggested for conservation approaches of native honey bee subspecies. A SNP-based method, Real-Time PCR-HRM, was implemented to determine the genetic sources of the spermathecal contents of honey bee queens. For this purpose, an mtDNA SNP found to vary across two subspecies (Caucasian and Italian) was identified. An mtDNA SNP was chosen due to concerns about analytical sensitivity during HRM analysis, though the future integration of additional nuclear markers would constitute an advance by enabling a more complete characterization of lineage introgression. Since sperm cells includes one copy of nuclear DNA as opposed to multiple copies of mtDNA (Meusel & Moritz, 1993), this initial work focused on an mtDNA marker. After SNP determination, queen bees were instrumentally inseminated with the semen mixture of two subspecies. Then, isolated DNA from the spermathecal content was genotypically identified by Real-Time PCR-HRM analysis based on mtDNA SNP marker.

This research was carried out at the Honey Bee and Beekeeping Laboratory, Department of Animal Science, Faculty of Agriculture, Ankara University. Field studies were done in 2013 and laboratory tests in 2014. In the spring season, one Caucasian (A. m. caucasica) and one Italian (A. m. ligustica) colony with a strong bee population were chosen and given to empty drone comb within a plastic queen excluder cage for drone rearing. The queen of the colony was confined to the plastic cage on the drone comb and was provided for egg laying into the empty drone cells. Drones were observed to emerge from sealed cells and were marked with paint (Edding 751) on the thorax. They were kept in a large drone banking colony for the sexual maturation period (10–12 days of age). The marked drones were used for semen collection and DNA sequencing.

The queen bees were reared through grafting. Upon reaching sexual maturity (6 days of age), they were instrumentally inseminated with 8 μl semen mixture of Caucasian and Italian drones (1:1) in order to obtain heterospermic spermathecal content. A Gilmont micrometer connected to a Harbo syringe was used for measuring semen volume. Instrumentally inseminated queens were grabbed after egg laying in the mating boxes and kept in ethyl alcohol (96%, −20ºC) until the dissection process. DNA was isolated via a kit (Invitrogen, PureLink® Genomic DNA kit) from both thorax of drones and spermathecal contents of queens. Sperms in the spermathecae of instrumentally inseminated queens (25 queen bees) were dissected under the stereo microscope (Leica, Z16 APO). The concentration and purity of isolated DNA were measured by Nanodrop Spectrophotometer (ND 2000), and DNA integrity was tested by running with 1% agarose gel electrophoresis. DNA samples were diluted to 7 ng/μl concentration with ddH₂O before PCR.

The positions of SNPs on mtDNA were determined by DNA sequencing in order to distinguish the Caucasian drones from the Italian drones. Eight Caucasian and six Italian drones were used for DNA sequencing. The SNPs were searched on seven mtDNA region. Sequences of the forward and reverse primer were for amplifying mtDNA regions; 5-GGCAGAATAAGTGCATTG-3 and 5-CAATATCATTGATGACC-3 (Garnery et al., 1993) for tRNA^leu-COΙI, 5-TTTTGTACCTTTTGTATCAGGGTTG-3 and 5-CTATAGGGTCTTATCGTCCC-3 (Hall & Smith, 1991), for LrRNA, 5-TTAAGATCCCCAGGATCATG-3 and 5-TGCAAATACTGCACCTATTG-3 (Hall & Smith, 1991) for COΙ, 5-TATGTACTACCATGAGGACAAATATC-3 and 5-ATTACACCTCCTAATTTATTAGGAAT-3 (Crozier, Koulianos, & Crozier, 1991) for Cytb, 5-TGATAAAAGAAATATTTTGA-3 and 5-TGAAACTATTATATAAATTG-3 (Arias & Sheppard, 1996) for ND2, 5-CAACATCGAGGTCGCAAACATC-3 and 5-GTACCTTTTGTATCAGGGTTGA-3 (Nielsen, Page, & Crosland, 1994) for 16s rRNA, 5-TCGAAATGAATAGGATACAG-3 and 5-GGTTGAGATGGTTTAGGATT-3 (Bouga et al., 2005) for ND5. Each PCR was performed in 30 μl volume containing dNTP (0.167 mM, Thermo Scientific dNTP Set), forward and reverse primers (0.167 μM, IDT), Taq Buffer (1x Thermo Scientific Fermentas), Taq polymerase (1 U, Thermo Scientific Fermentas) and 30 ng/reaction DNA. The thermal cycling parameters of PCR amplification were thirty-five cycles at 94ºC for 1 min., at 50–55ºC for 1 min., at 72ºC for 3 min. and at 72ºC for 15 min.

After SNPs were found on mtDNA regions of Caucasian and Italian drones, a new primer (5-TGAATTTGAGGTGGATTTTCA-3 and 5-CCAAGAGGATTAGATGATCCAG-3) was designed for the region where the SNP (T/C 11606 bp, Cau/Lig) were found in the Cytb region of the drones for use in HRM analyses. HRM analysis was performed on the Qiagen Rotor-Gene Q Real-Time Thermal Cycler (36 wells). PCR prior to HRM analysis was set up in 20 μl total volume. Each reaction contained a mixture of Roche-LightCycler® 480 High Resolution Melting Master, (1X), forward and reverse primer (0.2 μM), MgCl₂ (2.0 mM) and DNA (14 ng/reaction) isolated from spermathecal contents. The conditions of the Real-Time thermal cycler were initial denaturation at 94°C for 10 min., 50 cycles of denaturation at 94°C for 10 sec., annealing at 59°C for 15 sec., extension at 72°C 20 sec. And the final extension at 72°C for 2 min. After PCR, the HRM analysis was performed at 470 nm (excitation) and 510 nm (emission) wavelength channel. Before the HRM analysis, PCR products were kept to ensure heteroduplex formation 94ºC for 1 min., at 40ºC for 1 min. Melting curves were then generated by increasing the temperatures in steps 0.1ºC for two sec between 65 and 75ºC. Besides the HRM analysis of DNA isolated from the spermathecal contents, we examined the DNA obtained from Caucasian and Italian drone thoraces artificially mixed at different proportions (100% Caucasian, 75% Caucasian + 25% Italian, 50% Caucasian + 50% Italian, 25% Caucasian + 75% Italian, 100% Italian) for sensitivity of HRM analysis and the melting curve behaviors obtained from HRM analysis (homoduplex/heteroduplex).

The sequence results of the Caucasian and Italian drones were visualized through FinchTV (1.4.0) and the base sequences were aligned with Bioedit (Hall, 1999). Raw data of Real-Time PCR-HRM analyses was obtained from Rotorgene Q (2.1.0) software. The raw HRM data was transferred to the Gnumeric (1.12.35) software, geometric means were calculated from three points and melting curve graphics were generated.

One SNP (T/C 11606 bp, Cau/Lig) was found to discriminate the Caucasian and Italian drones. Fig. 1 shows the aligned nucleotide sequence of the mtDNA cytb region in two drone genotypes. The image of PCR products (143 bp) which contain detected SNP on the mtDNA Cytb region is shown in Fig. 2.

Fig. 1

Fig. 2

Isolated DNA from Caucasian and Italian drones was used as pure (navy blue: 100% Caucasian, claret red: 100% Italian) and as mixed in different proportions (pink: 75% Caucasian + 25% Italian, yellow: 50% Caucasian + 50% Italian, light blue: 25% Caucasian + % 75 Italian) to test the analytical sensitivity of HRM analyses. HRM analysis was found to be applicable for detecting the presence of foreign drone DNA. The melting curves of pure Caucasian and Italian samples (homozygotes) showed melting curves with one peak (homozygous). Fig. 3 shows the different melting curves of the pure Caucasian (navy blue) and Italian (claret red) DNA fragments amplified by qPCR. A difference of about 1°C was observed at the peak of the melting curves with one peak between pure Caucasian (navy blue) and Italian (claret red) DNA fragments. Caucasian alleles consisted of a T nucleotide at position 11606 in mtDNA leading to fragments with lower melting temperature (69°C), while the Italian alleles had C nucleotide, leading to fragments with higher melt peak (70°C). Due to imperfect binding, all mixed samples of two genotypes (heterozygotes) showed another early peak in the melting curves (heteropolymers) regardless of the mixture proportion (Fig. 3). The results of HRM analysis clearly showed promise for determining genotypes in spermatheca depending on the melting curve shape.

Fig. 3

The HRM analysis of isolated DNA from spermathecal contents of queen was successful in detecting the described SNP diversity of mtDNA associated with Caucasian and Italian genotypes. The spermathecal contents of twenty-five queen bees were analyzed with the Real-Time HRM method. The derivative curves obtained from the raw data in the HRM analysis of the isolated DNA from the spermathecal contents are presented in Fig. 4. All spermathecal DNA isolates of queens inseminated with the Caucasian and Italian sperm mixtures showed an early peak in the curves (blue arrow shows early peaks due to heterozygosity).

Fig. 4

In this study, a Real-Time PCR-HRM analysis was developed for targeting a genetic marker for the rapid discrimination of two drone genotypes in spermathecal contents in order to determine mating with drones or foreign maternal ancestry. Thus, it has been tested whether HRM analysis could be used as a genomic tool to determine gene flow to native honey bee populations from foreign subspecies. HRM analysis is widely used in heteroduplex analysis, SNP analysis, mutation screening and DNA fingerprint analysis. This method in which the melting thermodynamics of the double helix structure is monitored with the help of intercalating fluorescent dyes which can bound between double-stranded DNA.

The HRM analysis of all DNA samples isolated from the spermathecae were heteroduplexed. The melting curves of spermathecal contents were presented without any normalization procedure. HRM analysis results can be evaluated successfully without applying normalization to melting curves data (Druml & Cichna-Markl, 2014). In this research, all the queens were instrumentally inseminated with 50% Caucasian + 50% Italian sperm mixture (8 μl). Also, artificial mixtures of isolated DNA from Caucasian and Italian drone thoraces were studied for tests of HRM analysis. However, the sensitivity of the HRM analysis was not tested under the 25:75% DNA mixture (Fig. 3). Although the results may change according to the amplicon length and SNP class, heteroduplex structure is known to be detected in DNA mixture ratios as low as 1–5% by HRM analysis (Carillo et al., 2011; Do & Dobrovic, 2009; Ganopoulos et al., 2012; Krypuy et al., 2007).

Nucleotide variation within the mitochondrial genome is frequently used in the identification of subspecies and linage of honey bees (Collet, Arias, & Lama, 2007; Cornuet & Garnery, 1991; Smith et al., 1997; Solorzano et al., 2009). AT-rich mtDNA (~ 80%) (Crozier & Crozier, 1993) limits the success of primer design for HRM analysis (tm = ~ 60°C and fragment length <150 bp). Additionally, the class of the SNP to be studied should be considered, because HRM analysis shows higher success in Class 1 SNP mismatches (T/C and A/G) that change more the stability of the DNA chain (Liew, 2004; Ramirez et al., 2010).

One of the critical points in this research is DNA isolation from the contents of spermatheca, which is a small organ (~1 mm³). The small spermathecal content is a risk for obtaining sufficient quantity and quality of DNA to carry out molecular genetic analyses, but it was possible to obtain DNA in quantity and quality that would allow the execution of the studies (10–30 ng/μl). Microsatellites were used in past works that focused on the identification of drone subspecies that inseminated honey bee queens, (Oleksa et al., 2013). An additional study was conducted to develop a SNP panel for the detection of Africanization in honey bees (Chapman et al., 2017). While this assay is commendable, our understanding of honey bee population genetics should be substantially improved in order to apply this analysis type to the lineages of interest in this work. An HRM-based analysis may be more feasible in some cases as well given potential limitations on available laboratory technology and resources. Controlling the colony-level genetics of honey bee populations is important not only in terms of biology but also in determining whether the practices of breeding and conservation genetics are being done effectively. The monitoring of gene flow from foreign genetic sources to native honeybee populations depends on the identification of markers at the molecular level. In this research one SNP on mtDNA was used as a marker. The results showed that a genetic mixture from foreign sources to native honey bee populations could be quickly determined based on SNP genotyping with Real-Time HRM analysis in spermathecae.

Future works would benefit from nuclear markers as well as a broader set of markers. In this context, published genomes and SNP or indel (insertion/deletion) variant libraries of honey bee subspecies will be important. With additional genomic data, this method can be extended for screening additional genotypes, although HRM analysis typically requires an individual PCR reaction per marker for such genotyping. The use of more DNA regions for each genotype comparison will allow the power of the HRM genotyping method to be increased and the method may be made available for screening the success of honey bee conservation breeding efforts.