PCR-RFLP of mitochondrial DNA reveals two origins of Apis mellifera in Taiwan.

Beekeeping has been a highly valued industry in Taiwan. As a result, many subspecies of Apis mellifera have been introduced to Taiwan since 1911, leading to the hybridization of different subspecies. In order to know the matrilineal origins of Taiwan A. mellifera, a total of 280 samples collected from 33 apiaries throughout the island were examined. Using PCR-RFLP of four mitochondrial gene fragments, i.e., the non-coding region between tRNAleu and cytochrome c oxidase subunit II (intergenic tRNAleu-COII), cytochrome b (Cyt b), large subunit rRNA (Ls rRNA) and cytochrome c oxidase subunit I (COI), we only found two haplotypes exist in 280 samples. Haplotypes ababa and bbbaa account for 87% of these Western bees belonged to the Eastern European (C) lineage and 13% belonged to the Middle East (Z) lineage, respectively, with the latter being totally absent in northern Taiwan. African (A) and Mellifera (M) lineages, officially imported once in 1990s and 1930s respectively, were not detected. The identification of subspecies of A. mellifera and survey of their distribution on the island are expected to facilitate efficient breeding programs and establish a more booming beekeeping industry.

Introduction

Taiwan is located among the western pacific rim. To the north lie Japan; to the south are the Philippine Islands; to the west is Mainland China. The beekeeping industry in Taiwan owes its popularity and economic importance to the tropical/subtropical climate which provides an abundant supply of rich and diverse nectar plants (An et al., 2004). In addition to help pollinating crops, the 120,000 managed colonies of Apis mellifera on this island produced in 2015 ca. 12,000 tons of varietal honey as well as 430 tons of royal jelly (ranking the second in the world) (Bogdanov, 2015).

Official records have shown that A. mellifera subspecies, including armeniaca, capensis, carnica, caucasica, cypria, ligustica and mellifera were imported to improve the yield of various bee products during the period of 1911–1997 (Supplementary Information Table 1) (An et al., 2004, Sung et al., 2006). Since then, none has been imported officially. In Taiwan, these Western bees of some apiaries are known to give high yield of royal jelly while bees from other yards are better honey producer. It has been generally assumed that A. mellifera in Taiwan is composed of hybrids of different subspecies. Nevertheless, no study has been carried out up to the present to elucidate the makeup of honey bees managed by the beekeepers. We believe that identification of honey bee subspecies via molecular characterization is essential for more advanced management in apiculture in Taiwan.

Twenty-eight subspecies of A. mellifera have been identified morphologically. They are classified, according to their geographic origins, into four branches/lineages, i.e., Africa (A), Western Europe (M), South-Eastern Europe (C) and Middle East (O) (Ruttner, 1988, Engel, 1999, Gupta, 2014). However, phylogeographical studies have shown that variations of morphological characters of these subspecies resulted from adaption to local environments often complicate their proper identification. Therefore, molecular methods using DNA markers, such as mitochondrial DNA (mtDNA), microsatellites and single nucleotide polymorphisms (SNPs), have been developed. Among them, analysis of restriction site polymorphisms of mtDNA fragments, characterized by their maternal inheritance, has been widely used to identify subspecies of honey bees. The length polymorphism of tRNAleu-COII region of mtDNA, due to the absence/presence of P element and different copy numbers of Q element, which correlates with their evolutionary lineages, has been frequently used to classify A. mellifera lineages (Garnery et al., 1992, Garnery et al., 1993).

Among the five mitochondrial lineages identified in A. mellifera (Meixner et al., 2013), A, M and C roughly match with the morphological branches A, M and C (Cornuet and Garnery, 1991, Garnery et al., 1992). The original 4th mitochondrial lineage O from the Near East (Franck et al., 2000, Palmer et al., 2000) was split into mitochondrial C and Z lineages (Meixner et al., 2013), with a sub-lineage of A embedded in the latter (Alburaki et al., 2011, Alburaki et al., 2013). The 5th mitochondrial lineage Y has been identified in Ethiopia (Franck et al., 2001).

The existence of honey bees has been threatened by multiple factors, such as pesticide usage on crops, fragmentation and loss of habitat as well as pathogens and parasites (Antunez et al., 2009, Rosenkranz et al., 2010, Sanchez-Bayo et al., 2016). Colony Collapse Disorder (CCD), a combination of deadly effects, has widely occurred across Europe and America (Stokstad, 2007, van der Zee et al., 2012, Lee et al., 2015). Besides, the natural biodiversity of honey bees may be declined by the selective management of those subspecies that bring forth profits to the beekeepers. In view of the fact that CCD has not been observed in Taiwan and high yield of honey/royal jelly has been achieved, we intended to establish the matrilineal origins of these precious honey bees by analyzing four fragments of mtDNA, i.e., the non-coding region between tRNAleu and cytochrome c oxidase subunit II (intergenic tRNAleu-COII), cytochrome b (Cyt b), large subunit rRNA (Ls rRNA) and cytochrome c oxidase subunit I (COI) by PCR-RFLP and sequencing the tRNAleu-COII fragments of samples collected throughout the island.

Materials and methods

Sampling

A. mellifera samples representing 280 colonies were collected from northern (No. 1–7 apiary), central (No. 8–13 apiary), southern (No. 14–24 apiary) and eastern (No. 25–33 apiary) regions of Taiwan during September 2015–February 2016 (Fig. 1 and Supplementary Information Table 2). Five workers from each colony were preserved in 99% ethanol and stored at room temperature. One individual worker per colony was subjected for total DNA extraction.

Figure 1

The locations of thirty-three apiaries throughout Taiwan with the haplotype frequency. Total 280 A. mellifera colonies were sampled from northern (No. 1–7), central (No. 8–19), southern (No. 20–24) and eastern (No. 25–33) regions of Taiwan. Numbers in brackets represent number of colonies sampled at each apiary. The pie charts in four regions represent the haplotype frequency.

Total DNA extraction

Total DNA extraction was performed using Easy DNA High-Speed Extraction Tissue kit (Saturn Biotech, Australia) and following the manufacturer’s protocol. One foreleg from one individual bee was removed and dried for 5 min at 50 °C. Then it was homogenized in 100 μL mixture solution of Solution 1A, Solution 1B and Solution 2, and incubated for 30 min at 95 °C. After that, samples were spun down to remove cell debris and the supernatants with DNA were subjected to mtDNA analysis.

PCR amplification of mtDNA fragment

The fragments of the intergenic tRNAleu-COII region, Cyt b, Ls rRNA and COI of mitochondria were selected for PCR amplification with primer pairs of tRNAleu-COII region (Garnery et al., 1993), Cyt b (Crozier et al., 1991), Ls rRNA (with slight modification from Hall and Smith (1991)) and COI (Hall and Smith, 1991) (Table 1). PCR were performed using the thermal profile: 98 °C for 30 s, 35 cycles of 98 °C for 10 s, 40/45/50 °C for 15 s, 72 °C for 50 s and finished with 72 °C for 5 min. Each PCR reaction contains 5× PCR buffer, 0.1 μM of each primer pair, 10 mM dNTP, 2 μL of DNA and 0.25 μL of Phusion™ DNA Polymerase (NEB, USA) to make a final volume of 25 μL.

Table 1

Primers used in this study.

Mitochondrial genes	Primers	Tm (°C)	References
tRNAleu-COII	Forward (Position 3363): 5′ GGCAGAATAAGTGCATTG 3′	50	Garnery et al. (1993)
	Reverse (Position 3934): 5′ CAATATCATTGATGACC 3′
Cyt b	Forward (Position 11400): 5′ TATGTACTACCATGAGGACAAATATC 3′	45	Crozier et al. (1991)
	Reverse (Position 11884): 5′ ATTACACCTCCTAATTTATTAGGAAT 3′
Ls rRNA	Forward (Position 13708): 5′ CTATAGGGTCTTATCGTCCCATA 3′	40	Hall and Smith (1991)
	Reverse (Position 14443): 5′ GTACCTTTTGTATCAGGGTTGA 3′
COI	Forward (Position 1908): 5′ TTAAGATCCCCAGGATCATG 3′	45	Hall and Smith (1991)
	Reverse (Position 2951): 5′ TGCAAATACTGCACCTATTG 3′

(The position numbers of all primers are based on the mitochondrial sequence of A. mellifera ligustica L06178).

Digestion with restriction enzymes

Each PCR amplification product was digested with the appropriate restriction enzyme (Table 2) in a final volume of 20 μL, containing 10× reaction buffer, 17 μL of DNA and 1 μL of restriction enzyme (NEB, USA), with overnight incubation at 37 °C.

Table 2

Summary of restriction enzyme digested patterns of four mitochondrial lineages.

Mitochondrial genes	Restriction enzymes	Lineages of A. mellifera				References
		M	A	C	O
tRNAleu-COII	DraI	+	+	+	+	Garnery et al. (1993)
Cyt b	BglII	+	−	+	+	Crozier et al. (1991)
	HinfI	+	+	−	+	Suppasat et al. (2007)
Ls rRNA	EcoRI	−	−	+	−	Zaitoun et al. (2008)
						Suppasat et al. (2007)
COI	HinCII	+	−	−	−	Hall and Smith (1991)

Electrophoresis

The fragment sizes of PCR amplification and restriction enzyme digestion products were analyzed by the QIAxcel Advanced system (Qiagen, USA). DNA fragments were separated in QIAxcel gel cartridge using QIAxcel Method AM420 and the collected signal data were converted to gel image by the QIAxcel ScreenGel Software. The size marker and alignment marker used in this system are 100 bp/2.5 kp and 15 bp/5 kb, respectively.

Cloning and sequencing

The amplified gene fragments of tRNAleu-COII were cloned into pGEM-T vector (Promega, USA) and sequenced from both directions. The resulting sequences were compared to published sequences available in Genbank. DNA sequence alignment was performed, using the Clustal W algorithm of Vector NTI® software V10 (Invitrogen, USA). In the phylogenetic analysis, the nucleotide sequences were aligned to construct a phylogenetic tree by using the neighbor-joining method in the MEGA6 program with a Tamura 3-Parameter model, pairwise deletion of gaps, and 1000 bootstrap replicates (Tamura et al., 2013).

Results and discussion

PCR amplification of mtDNA fragments of A. mellifera

All four mtDNA fragments i.e., the intergenic tRNAleu-COII region, Cyt b, Ls rRNA, and COI were successfully amplified for all samples (Fig. 2). Two amplicon sizes, i.e., 580 and 930 bp for the tRNAleu-COII region were observed. The amplicon of ca. 580 bp appeared in 245 samples while 930 bp amplicon was observed in the remaining 35 samples. This result is similar to that of a study of Thailand Western bees which produced two amplicon sizes 573 and 837 bp of tRNAleu-COII fragments (Suppasat et al., 2007). The details of sequencing analysis of this region are provided in the following section. The amplicon sizes of Cyt b, Ls rRNA and COI are ca. 500, 800, and 1100 bp respectively, matching with those from previous studies (Crozier et al., 1991, Hall and Smith, 1991).

Figure 2

PCR amplification of four mtDNA fragments of A. mellifera. Four gene fragments, i.e., tRNAleu-COII, Cyt b, Ls rRNA and COI fragments were amplified from mtDNA. Lane 1: DNA size marker (100 bp/2.5 kp); Lane 2 & 3: tRNAleu-COII fragments; Lane 4: Cyt b; Lane 5: Ls rRNA; Lane 6: COI.

RFLP of four mtDNA regions in 280 samples of A. mellifera

Table 2 summarizes the patterns of restriction enzyme digests of four mtDNA regions used in the characterization of mitochondrial lineages of A. mellifera. The presence of a BglII restriction site in Cyt b fragment in all 280 samples excludes A lineage (pattern ‘b’, Fig.3B); and the absence of HinCII site in COI fragment of all samples indicates that none of them belongs to M lineage (pattern ‘a’, Fig.3E) (Crozier et al., 1991, Hall and Smith, 1991). This leaves only the possibility that these 280 samples belonged to either lineage C or Z; and the presence/absence of two restriction sites, HinfI site in Cyt b fragment and EcoRI site in Ls rRNA fragment, is used to differentiate these two lineages (Suppasat et al., 2007, Zaitoun et al., 2008). Thirty-five samples, with HinfI site in Cyt b fragment (pattern ‘b’) and a lack of EcoRI site in Ls rRNA fragment (pattern ‘a’), were designated as mitochondrial lineage Z (Lane 2 and Lane 1 in Fig.3C and D respectively); and the remaining 245 samples, having EcoRI site in Ls rRNA fragment (pattern ‘b’) but not HinfI site in Cyt b fragment (pattern ‘a’), belonged to C lineage (Lane 1 and Lane 2 in Fig.3C and D respectively). Additionally, the digestion patterns of 580 and 930 bp fragments of tRNAleu-COII with DraI (pattern ‘a’ and ‘b’, Fig.3A) match with the result of Thailand Western bees ThaiA and ThaiB, which have been demonstrated as East Mediterranean (i.e., lineage C) and Middle East (i.e., lineage Z) lineage respectively (Suppasat et al., 2007).

Figure 3

Restriction enzyme digested patterns of four mtDNA fragments of A. mellifera. (A) tRNAleu-COII fragments digested by DraI; (B) Cyt b digested by BglII; (C) Cyt b digested by HinfI; (D) Ls rRNA digested by EcoRI; (E) COI digested by HinCII. DNA size marker (100 bp/2.5 kp) was used in this work.

The patterns generated from all five restriction digests were combined into composite haplotypes. The composite haplotypes indicate the restriction patterns of five PCR product/restriction enzyme combinations: [tRNAleu-COII/DraI] [Cyt b/BglII] [Cyt b/HinfI] [Ls RNA/EcoRI] [COI/HinCII]. Only two composite haplotypes, i.e., ababa and bbbaa, were found in the 280 colonies of Taiwan A. mellifera, which further refer to lineage C and Z, respectively.

Sequence analysis of tRNAleu-COII region

Since PCR amplification of tRNAleu-COII fragment from mtDNA produced one amplicon (580 bp) in all 245 samples of lineage C and a amplicon size of 930 bp in all 35 samples of Z lineage, we decided to sequence the tRNAleu-COII region in some selected samples from both lineages, i.e., 14 samples of C lineages from 14 cities (one sample per city) and 16 samples of Z lineages from Center, South and East of Taiwan. A bootstrap consensus tree was constructed using aforementioned 30 samples with two reference genes (JQ977703 and FJ743633), and can be clustered into two groups, C and Z lineages (Fig.4A). The blast result shows that the 580 bp amplicon was identical to the corresponding sequence of lineage C (JQ977703) (Fig.4B) and the sequence identity among 14 samples of C lineage is 98%. The sequence of 930 bp amplicon was identical to lineage O (FJ743633) which subsequently merged into lineage Z (Franck et al., 2000, Franck et al., 2001, Alburaki et al., 2011) and the sequence identity among 16 samples of Z lineages is 98%. The presence in this amplicon of a 69 bp P element and the RFLP data presented above fully confirm that these 35 samples belong to lineage Z (Fig.4B).

Figure 4

Phylogenetic analysis of tRNAleu-COII gene fragments. (A) The tRNAleu-COII gene fragments of 32 samples, including 14 samples of C lineages, 16 samples of Z lineages, C lineage reference gene (JQ977703) and Z lineage reference gene (FJ743633) were aligned to construct a phylogenetic tree using the neighbor-joining method in the MEGA6 program with a Tamura 3-Parameter model, pairwise deletion of gaps, and 1000 bootstrap replicates. (B) Alignment analysis of 580 and 930 bp fragments of tRNAleu-COII. Highly conserved residues are marked in yellow and the P element is framed in red color.

Distribution and abundance of A. mellifera lineages C and Z

According to Ruttner’s scenario, the haplotype of the 17 imported A. mellifera subspecies in Taiwan consisted of 70% of C lineage, 18% of O lineage, 6% of M lineage and 6% of A lineage (Supplementary Information Fig. 1A) (Ruttner, 1988). Converting into the mitochondrial classification system and taking into consideration of the reclassification of subspecies O lineage into C and Z lineages, the haplotype frequencies of Taiwan Western bee would be 76% of C lineage, 6% of Z lineage, 6% of M lineage and 6% of A lineage (Supplementary Information Fig. 1B) (Meixner et al., 2013). In 280 samples examined in this study, 245 and 35 samples were identified, respectively, as mitochondrial C and Z lineages (Fig. 5). The bees of C lineage distribute throughout the entire island, while the bees of Z lineage, absent in the northern region, are mainly present in the southern and eastern regions (Fig. 1).

Figure 5

Haplotype frequency of 280 A. mellifera samples. The two composite haplotypes, i.e., ababa and bbbaa, refer to the mitochondrial C and Z lineage, respectively.

Conclusions

This is the first investigation of the origin of Western bees in Taiwan. The population of Taiwan Western bee mainly consists of Eastern Europe lineage which widely distributes throughout Taiwan. Additionally Taiwan Western bee population also contains a small proportion of Middle East lineage. According to haplotype frequency analysis, the Middle East lineage has increased from the initial 6% to 12.5%, indicating that this lineage population has expanded in Taiwan. The reason for the absence of Western Europe and Africa lineages is uncertain, though one might presume that it was related to the importation of low number of colonies or other beekeeping issues. Further studies will be performed to address the issues of expansion and distribution of Middle East lineage in Taiwan.