DOI: 10.13102/sociobiology.v65i3.2876Sociobiology 65(3): 482-490 (September, 2018) Open access journal: http://periodicos.uefs.br/ojs/index.php/sociobiology ISSN: 0361-6525 Genetic Variation in Iranian Honey bees, Apis mellifera meda Skorikow, 1829, (Hymenoptera: Apidae) Inferred from PCR-RFLP Analysis of two mtDNA Gene Segments (COI and 16S rDNA) Introduction Intraspecific taxonomy of the honey bee Apis mellifera L. has been based mainly on morphology. At present, 29 subspecies of A. mellifera are recognized on the basis of morphometric characters (Ruttner, 1988, 1992; Sheppard et al., 1997; Sheppard & Meixner, 2003; Arias & Sheppard, 2005; Meixner et al., 2011 Rahimi et al., 2017). The western honey bee originated in Asia and entered Africa and Europe in three distinct evoluationary branches: branch (A), which included the subspecies from Africa (A. m. lamarckii, A. m. Abstract In this study, the genetic structure of Iranian honey bee (Apis mellifera meda) populations, mainly obtained from all of regions, were investigated at two different mitochondrial regions. A total of 300 worker bees were collected from 20 different populations in 20 different locations. Portions of the mitochondrial 16S ribosomal RNA (16S rDNA) and cytochrome C oxidase I (COI) genes were amplified by PCR and then subjected to RFLP pattern analysis using 8 restriction enzymes. Nucleotide polymorphisms were revealed using restriction enzyme Sau3A I, Ssp I and Taq I in COI and Bsp143I, Ssp I and Dra I in the 16S rDNA gene segment. In this study, 3 novel composite genotypes (haplotypes) were found in Iranian honey bee populations. The average haplotype diversity (h) within populations was 0.0405. Heterozygosity values, Shannon index and the number of alleles of Iranian honey bee populations were low that could be caused by low definite geographic structure of Iranian honey bee populations. Genetic distance (D) values were found to be low (0.0–0.0011) within Iranian honey bee populations. Cluster analysis based on UPGMA method revealed that all populations and samples groups be in one cluster. Also, the phylogenetic tree based on Neighbor-joining method divided 29 subspecies of honey bee to 5 distinct clusters. The Iranian subspecies honey bee composed of a shared clade with subspecies of Eastern Mediterranean, Near East and Eastern parts of Middle East (O branch). This result is very useful for the control of conservation of local honey bees, as the movement of colonies across the border line of these neighboring countries, may affect the genetic structure of honey bee populations. Sociobiology An international journal on social insects A Rahimi1,2, A Mirmoayedi1, D Kahrizi3, L Zarei3, S Jamali1 Article History Edited by Marco Costa, UESC, Brazil Received 04 February 2018 Initial acceptance 03 June 2018 Final acceptance 05 September 2018 Publication date 02 October 2018 Keywords Apis mellifera; Genetic structure, mtDNA, Iran. Corresponding author Ataollah Rahimi Department of Plant Protection Faculty of Agriculture, Campus of Agriculture and Natural Resources Razi University, Kermanshah, Iran. E-Mail: Rahimi.ata.1@gmail.com yemenitica, A. m. scutellata, A. m. litorea, A. m. adansonii, A. m. capensis), branch (M) which included the subspecies of North African and West European (A. m. mellifera, A. m. iberica, A. m. intermissa and branch (C) which included the subspecies from Eastern Europe, Northern Mediterranean and Middle East. Afterwards, subspecies in branch C were divided into two groups, branch C included A. m. carnica, A. m. ligustica, A. m. macedonica, A. m. cecropia and A. m. sicula, branch O included the Near and Middle Eastern subspecies (A. m. caucasica, A. m. armeniaca, A. m. meda, A. m. anatoliaca, A. m. syriaca, A. m.cypria, A. m. adami) (Ruttner, 1988). The 1- Department of Plant Protection, Faculty of Agriculture, Campus of Agriculture and Natural Resources, Razi University, Kermanshah, Iran 2- Department of Animal Science, College of Agriculture - Kifri, Garmian University, Kalar, As Sulaymaniyah, KRG of Iraq 3- Department of Plant Breeding and Biotechnology, Faculty of Agriculture, Campus of Agriculture and Natural Resources, Razi University, Kermanshah, Iran RESEARCH ARTICLE - BEES Sociobiology 65(3): 482-490 (September, 2018) 483 phylogenetic relationships based on molecular data (Cornuet & Garnery, 1991; Garnery et al., 1992; Arias & Sheppard, 1996) agree in general with those obtained with confirmed using mitochondrial and microsatellite variability (Franck et al., 2000; Palmer et al., 2000; Kandemir et al., 2006; Rahimi, 2014a,b). Southwest Asia included Iran is a region of high morphological diversification and evolution for honey bees. One clearly distinct race has been evolved in this region, which include a diversity of habitats. Honey bee race in this region include the subspecies A. m. meda (Ruttner, 1988). Honey bee subspecies from Iran was studied extensively using morphometric and izoenzymic analysis (Moradi & Kandemir, 2004; Rahimi & Asadi, 2010; Rahimi & Mirmoaydi, 2013; Rahimi et al., 2014a,b; Rahimi et al., 2015a; Rahimi et al., 2017); molecular markers, mitochondrial DNA (mtDNA) analysis (Rahimi et al., 2014a,b; Rahimi et al., 2015b; Rahimi et al., 2016); microsatellite and RAPD analysis (Royan et al., 2007; Kamrani et al., 2012; Rahimi et al., 2014a,b). More recently, genetic systems such as allozymes (Nunamaker & Wilson, 1982; Badino et al., 1988), nuclear DNA (Hall, 1990; Tarès et al., 1993), mitochondrial DNA (mtDNA) (Moritz et al., 1986; Smith et al., 1989, 1991; Hunt & Page, 1992; Garnery et al., 1993; Oldroyd et al., 1995; Arias and Sheppard, 1996; Pedersen, 1996; De la Rùa et al., 2000) and microsatellites (Estoup et al., 1993; Garnery et al., 1998) have been used to study honey bee diversification. Such analysis of population genetic differentiation at the molecular level contributes to better understanding of honey bee population structure by allowing the comparison of morphological, behavioural, geographical and molecular variation. The mtDNA is a favourite tool in systematic and population biology. It is generally maternally inherited without recombination. Only maternal inheritance of mtDNA has been demonstrated for honey bees (Smith, 1991; Meusel & Moritz, 1993; Arias & Sheppard, 1996; Francisco et al., 2001; Pinto et al., 2003). Iran is a vast country comprising a land area about 1648195 km2. From the geographical point of view, Iran has a very diverse climatic condition with diverse patterns of plant communities forming four climatic zones simultaneously; therefore it has a high potential for beekeeping and honey production. Therefore, genetic structure maintenance and preservation of Iranian honey bees is the first step to explore of this potential. The main aim of the present research was to determine the level of genetic differentiation among Iranian honey bee populations as discriminated using PCR-RFLP pattern analysis of the COI and 16S rDNA gene regions of mtDNA, and also the results of this study were compared with the results of other earlier mitochondrial studies of honey bees, such comparison thereby allowing much more complete estimation of the genetic structure of Iranian honey bee populations than until now previously possible using morphometrics alone. Fig 1. Geographical locations of the 20 sampled honey bee populations in Iran. N is the number of sampled bees per population. A Rahimi, A Mirmoayedi, D Kahrizi, L Zarei, S Jamali – Genetic variation in Iranian honey bee484 Materials and Methods Sampling A total of 300 young adult worker bees from 150 Apis mellifera meda colonies were sampled in 100 different localities from 20 Iranian provinces during June to October 2014 (Fig 1, Table 1). Samples were taken from honey bee colonies in most active apiaries from five cities in each province, one apiary in each city and one to three hive per apiary. Young adult worker bees directly were collected from the broad areas on the combs. The samples were immediately killed by immersion in absolute ethanol and kept at –20 ºC until DNA extraction. Molecular analysis DNA extraction Total genomic DNA was extracted from the head and thorax sections of each honey bee worker, using the salting out method described by Aljianabi and Martinez, (1997) with slight modifications and then stored at −20 ºC. RFLP Analysis The mtDNA variation was analyzed by RFLP’s, performed on PCR amplified products. Mitochondrial regions were amplified according to Bouga et al. (2005). Two sets of primers were used for amplifying 16S rDNA and COI gene regions. These were 5’- CAACATCGAGGTCGCAAACATC-3’ and 5’-GTACCTTTT GTATCAGGGTTGA-3’ for 16S rDNA and 5’-GATTA CTTCCTCCCTCATTA-3’ and 5’-AATCTGGATAGTCTG AATAA-3’for COI segment. PCR was run in a total volume 25 µl of the following reaction mixture: 2.5 µl of 10 X reaction buffer with KCl as provided by the manufacturer (Fermentas Life Sciences, Vilnius, Lithuania), 1 mM MgCl2, 0.5 mM of dNTP mix, 1 µM of each primer, 2 U of Taq polymerase and 20 ng of total purified honey bee DNA. For each primer pair, the following reaction profile was used: initial denaturation 94 °C for 4 min, 35 cycles of 94 °C for 1 min, annealing at 55 °C for 1 min, and extension at 72 °C for 2 min, followed by a final extension step at 72 °C for 10 min. The amplified products obtained were next electro- phoresed on 1% agarose gel to verify the size of the fragment. Amplified mtDNA regions from one individual of each bee were digested with 8 restriction enzymes to check for the presence of recognition sites. The informative restriction enzymes were then analyzed using one individual from each colony. The informative restriction enzymes used for the 16S rDNA gene fragment were: Sau3A I, Ssp I, Dra I, and EcoR I, and for the COI gene fragment Sau3A I, Ssp I, Taq I and NcoI. The digested fragments were separated electrophore- tically on 2% or 3% agarose gels in 1× TBE buffer, stained with Table 1. Sampling localities, geographical positions and the number of honey bees used for PCR-RFLP analysis in Iranian honey bee. Province/ locations Abbreviation of the province/ locations The number of apiary The number of colonies The number of bees Geographical position Latitude Longitude altitude Kermanshah KER 5 7 14 34° 29′ 15.62 N 047° 05′ 57.65 E 1664 Hamadan HAM 5 8 16 34° 43′ 24.64 N 048° 31′ 01.27 E 2464 Ilam ILA 5 5 10 33° 34′ 31.97 N 046° 27′ 26.44 E 1548 Lorestan LOR 5 6 12 33° 32′ 01.81 N 048° 14′ 12.03 E 1802 Esfahan ESF 5 10 20 32° 36′ 56.09 N 051° 34′ 04.20 E 1618 Chaharmahal and Bakhtiari CHA 5 7 14 32° 28′ 19.40 N 050° 06′ 30.53 E 2582 Fars FAR 5 10 20 28° 59′ 19.33 N 053° 39′ 56.47 E 1563 Sistan and Baluchestan SIS 5 5 10 29° 30′ 53.89 N 060° 50′ 17.17 E 1407 Kerman KRM 5 6 12 29° 18′ 19.25 N 057° 08′ 29.72 E 2838 Kurdistan KUR 5 7 14 34° 42′ 48.06 N 046° 51′ 34.16 E 1775 West Azerbaijan WEA 5 10 20 37° 34′ 04.14 N 044° 44′ 17.19 E 2138 East Azerbaijan EAA 5 10 20 38° 06′ 13.12 N 046° 11′ 10.99 E 1345 Ardabil ARD 5 9 18 38° 13′ 58.28 N 048° 10′ 05.22 E 1455 Zanjan ZAN 5 6 12 36° 42′ 37.00 N 048° 22′ 25.99 E 1546 Tehran TEH 5 6 12 35° 43′ 44.11 N 052° 07′ 00.70 E 2233 Razavi Khorasan REK 5 7 14 36° 25′ 56.06 N 059° 31′ 26.22 E 922 North Khorasan NOK 5 6 12 37° 29′ 16.88 N 057° 09′ 56.89 E 1454 Golestan GOL 5 5 10 36° 50′ 47.46 N 054° 40′ 12.88 E 181 Mazandaran MAZ 5 10 20 36° 25′ 22.19 N 052° 14′ 39.71 E 135 Gilan GIL 5 10 20 36° 50′ 23.28 N 049° 25′ 57.83 E 329 Total 150 300 Sociobiology 65(3): 482-490 (September, 2018) 485 ethidium bromide, visualized under UV light and photographed using a Vilber Lourmat gel imaging system. The sizes of DNA fragments were compared to the PCR marker (Promega G316A, Promega Corp.) run on the same gel and were calculated using DNA frag 3.03 (Nash, 1991) program. Data analysis Composite genotypes for each individual were then defined from all the restriction patterns of the two mtDNA segments. The restriction fragment data were converted to restriction data (gain or loss of restriction site). The POPGENE software version 1.31 was used for estimating population genetic structure. The genetic distance from restriction site data (Nei & Tajima, 1981; Nei & Miller, 1990) was estimated using the REAP computer package (McElroy et al., 1991). Phylogenetic tree was constructed by the Neighbor-joining method, based on results of genetic distance from restriction site data, using the PAST package version 2.02. Results from analogous study on honey bees from different areas for other subspecies of honey bee were included in the above mentioned statistical process. Results The sizes of the PCR-amplified mtDNA regions for all populations studied were found to average 964 bp for 16S rDNA and 1028 bp for COI mtDNA gene regions. Eight restriction enzymes were found to have at least one restriction site in the amplified 16S rDNA and COI regions, respectively. Observed fragment patterns generated by each restriction enzyme for the two mtDNA regions are summarized in Table 2. Sau3A I, Ssp I and Taq I restrictions in COI and Sau3A I, Ssp I and Dra I restriction in the 16S rDNA gene each generated two different restriction profiles in Iranian honey bee populations (Table 2). Diagnostic and novel patterns were revealed in the East Azerbaijan (EAA), West Azerbaijan (WEA) and Sistan and Baluchestan (SIS) populations after the digestion of 16S rDNA and COI segment with the restriction enzymes Bsp143I, Ssp I and Dra I and Sau3A I, Ssp I and Taq I, respectively (pattern type A). The three different and novel haplotypes (composite genotypes) which were detected in the twenty populations studied and the haplotype frequencies and haplotype diversity values are presented in Table 3. The three novel haplotypes were found in East Azerbaijan (EAA), West Azerbaijan (WEA) and Sistan and Baluchestan (SIS) samples. The average haplotye diversity (h) within populations was 0.0405 (Table 3). The number of alleles (Na) and the effective number of alleles (Ne), the observed and expected heterozygosities (Ho and He) and Shannon index (I) per COI and 16S rDNA COI 16S rDNA Restriction enzyme Patterns observed (bp) Restriction enzyme Patterns observed (bp) Sau3A I A: 338, 315, 168, 75, 72, 60 B: 398, 315, 168, 75, 72 Sau3A I A: 964 B: 545, 419 Ssp I A: 520, 210, 170, 80, 48 B: 520, 258, 170, 80 Ssp I A: 515, 335, 114 B: 630, 334 Taq I A: 350, 262, 244, 172 B: 432, 352, 244 Dra I A: 314, 212, 118, 115, 90, 70, 45 B: 359, 212, 118, 115, 90, 70 NcoI 1028 EcoR I 494, 470 Table 2. Restriction fragment patterns generated from analysis of COI and 16S rDNA gene segments of Iranian honey bee. Fig 2. Cluster analysis of 300 honey bee individuals from twenty provinces of Iran based on UPGMA method. A Rahimi, A Mirmoayedi, D Kahrizi, L Zarei, S Jamali – Genetic variation in Iranian honey bee486 mtDNA gene and per restriction enzyme in Iranian honey bee populations are shown in Table 4. The genetic distance (D) values were found to be low (0.0 – 0.0011) within Iranian honey bee populatioins (Table 5). Phylogenetic trees are generally plotted to find the genetic distances between the collected samples of honey bees. Cluster analysis based on UPGMA method between all samples Haplotype Composite genotype Haplotype diversity (h) NCOI 16S rDNA Sau3A I Ssp I Taq I NcoI Sau3A I Ssp I Dra I EcoR I Kermanshah A A A A A A A A 0.00 14 Hamadan A A A A A A A A 0.00 16 Ilam A A A A A A A A 0.00 10 Lorestan A A A A A A A A 0.00 12 Esfahan A A A A A A A A 0.00 20 Chaharmahal and Bakhtiari A A A A A A A A 0.00 14 Fars A A A A A A A A 0.00 20 Sistan and Baluchestan B B B A B B B A 0.18 10 Kerman A A A A A A A A 0.00 12 Kurdistan A A A A A A A A 0.35 14 West Azerbaijan B B B A B B B A 0.28 20 East Azerbaijan B B B A B B B A 0.00 20 Ardabil A A A A A A A A 0.00 18 Zanjan A A A A A A A A 0.00 12 Tehran A A A A A A A A 0.00 12 Razavi Khorasan A A A A A A A A 0.00 14 North Khorasan A A A A A A A A 0.00 12 Golestan A A A A A A A A 0.00 10 Mazandaran A A A A A A A A 0.00 20 Gilan A A A A A A A A 0.00 20 Haplotype diversity (h) 0.0405 Table 3. Composite genotypes (haplotypes), haplotype diversity and sample size of all the populations studied. mtDNA gene Restriction enzymes Na Ne HO He I COI NcoI 1 1 0.001 0.001 0.001 TaqI 3 1 0.001 0.001 0.001 SspI 5.1 1.001 0.002 0.002 0.0009 Sau3AI 6.1 1.002 0.004 0.004 0.0033 16S rDNA EcoRI 2.2 1.003 0.006 0.006 0.0016 DraI 6.1 1.003 0.003 0.003 0.0033 PstI 2.1 1.003 0.005 0.005 0.0033 TaqI 2.1 1.003 0.005 0.005 0.0033 Table 4. Genetic parameters: the number of alleles (Na) and the effective number of alleles (Ne), the observed and expected heterozygosities (Ho and He) and Shannon index (I) per COI and 16S rDNA mtDNA gene and per restriction enzymes in Iranian honey bee populations. in this study shown in Fig 2. The results of this cluster showed that all populations and samples groups be in one cluster. Phylogenetic tree based on RFLP analysis was drawn using Neighbor-joining method to compare phylogenetic relationships Iranian subspecies honey bee with other honey bee subspecies. Phylogenetic tree showed that 29 honey bee subspecies could be divided into 5 distinct groups (Fig 3). Discussion In a recent study of Iranian honey bee samples, we earlier showed that all honey bee samples then tested belonged to the East Mediterranean (C) lineage as found using several restriction enzymes and morphological characters (Rahimi, 2015b; Rahimi et al., 2017). In the present study, we have used different two mtDNA regions with different samples of the mitochondrial genome to verify the nature and distribution of genetic variation within and between Iranian honey bee populations. The results of Fig 3 showed that Iranian honey bee belong to evolutionary lineages (C) and confirmed the results of previous studies. It is remarkable that several technical improvements introduced in beekeeping management may have interfered with the natural distribution of populations. The importation of foreign Sociobiology 65(3): 482-490 (September, 2018) 487 W E A E A A A R D K U R Z A N G IL M A Z G O L N O K R E K K E R H A M IL A T E H L O R E SF C H A FA R K R M SI S W E A ** * E A A 0. 00 00 ** ** A R D 0. 00 03 0. 00 03 ** ** K U R 0. 00 00 0. 00 00 0. 00 03 ** ** Z A N 0. 00 04 0. 00 04 0. 00 08 0. 00 03 ** ** G IL 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 ** ** M A Z 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 ** ** G O L 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 ** ** N O K 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 ** ** R E K 0. 00 00 0. 00 00 0. 00 04 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** K E R 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** H A M 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** IL A 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** T E H 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** L O R 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** E SF 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** C H A 0. 00 00 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 0. 00 00 ** ** FA R 0. 00 03 0. 00 00 0. 00 10 0. 00 06 0. 00 11 0. 00 03 0. 00 00 0. 00 03 0. 00 00 0. 00 04 0. 00 03 0. 00 00 0. 00 04 0. 00 03 0. 00 00 0. 00 04 0. 00 01 ** ** K R M 0. 00 01 0. 00 01 0. 00 03 0. 00 01 0. 00 05 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 01 0. 00 07 ** ** SI S 0. 00 06 0. 00 06 0. 00 10 0. 00 06 0. 00 11 0. 00 06 0. 00 06 0. 00 06 0. 00 06 0. 00 07 0. 00 06 0. 00 06 0. 00 06 0. 00 06 0. 00 06 0. 00 06 0. 00 01 0. 00 06 0. 00 01 ** ** T ab le 5 . G en et ic d is ta nc e (D ) v al ue s be tw ee n pa ir s of th e Ir an ia n ho ne y be e po pu la tio ns . Fig 3. Phylogenetic tree between 29 honeybee subspecies based on Neighbor-joining method. queens and the practice of moving colonies several times per year are factors that can affect the genetic structure of a local honey bee population through genetic introgression (Garnery et al., 1998). In comparison our results with Bouga et al. (2005), Kekeçoglu et al. (2009) and Ozdil et al. (2012), we find similar restriction profiles with generally small base differences in the fragments, except for Sau3A I, Ssp I and Taq I and Bsp143I, Ssp I and Dra I digestions in COI and 16S rDNA gene regions, respectively (Table 2). In these regions, different restriction profiles were detected in our samples or populations. Sau3A I, Ssp I and Taq I restriction in COI revealed six or five, five or four and four or three restriction sites in this study, respectively, whereas in others subspecies such as A. m. anatoliaca, A. m. caucasica and A. m. meda (in Turkey), restriction analysis were previously reported to have revealed four or six, three or two and two restriction sites for these enzymes, respectively (Bouga et al., 2005, Kekeçoglu et al., 2009, Ozdil et al., 2012). Correspondingly, the restriction of the 16S rDNA gene region with Sau3A I, Ssp I and Dra I revealed two or one, three or two and seven or six sites (Table 2), compared with the three or two, two and two or seven as reported in Turkish honey bees (Bouga et al., A Rahimi, A Mirmoayedi, D Kahrizi, L Zarei, S Jamali – Genetic variation in Iranian honey bee488 2005, Kekeçoglu et al., 2009, Ozdil et al., 2012). Comparing our findings with these of Ozdil et al. (2012), honey bees in Iran and Turkey (specially A. m. meda) were found to be well discriminated; it is noted that for the first time are reported diagnostic patterns that distinguish the honey bee populations of these neighboring countries. Honey bees from Iran are also discriminating from other subspecies of honey bees in Turkey, Iraq and Syria. This study indicated that Iranian honey bee populations were characterized by a low variability in the number of alleles (Na) and the effective number of alleles (Ne), the observed and expected heterozygosities (Ho and He) and Shannon index (I). According to results of cluster analysis (Fig 2), it should be concluded that the honey bees of these twenty populations have formed one population, and there is no special divisions among them. Also, the cluster analysis showed that honey bees from different cities and provinces are scattered across the branches of the phylogenetic tree and apparently belong to the same population. This information showed that genetic diversity between Iranian honey bee populations is low, so we need a comprehensive breeding program for breeding of Iranian honey bee. However, the above mentioned results are very useful for the conservation of Iranian honey bee, but further follow up studies are needed to characterize mitochondrial DNA variation in the honey bees of these regions along with morphometric analyses in order to accurately characterize the particular regional honey bee populations. Acknowledgements The authors appreciate the kindness of those beekeepers who allowed us to collect honey bee samples from their apiaries and equally thank the authorities of Plant Protection and Biotechnology departments of Razi university who let us a free access to laboratory facilities which led to fulfill of the present study. References Aljianabi, S.M. & Martinez, I. (1997). Universal and rapid salt extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Research, 25: 4692-4693. Arias, M.C. & Sheppard, W.S. (1996). Molecular phylogenetics of honey bee subspecies (Apis mellifera L.) inferred from mitochondrial DNA sequence. Molecular Phylogenetics and Evolution, 5: 557-566. doi: 10.1006/mpev.1996.0050. Arias, M. & Sheppard, W. (2005). Phylogenetic relationships of honey bees (Hymenoptera: Apinae: Apini) inferred from nuclear and mitochondrial DNA sequence data. Molecular Phylogenetics and Evolution, 37: 25-35. doi:10.1016/j. ympev.2005.02.017. Badino, G., Celebrano, G., Manino, A. & Ifantidis, M. D. (1988). Allozyme variability in Greek honeybees (Apis mellifera L.). Apidologie, 19: 377-386. Bouga, M., Harizanis, P.C., Kilias, G. & Alahiotis, S. (2005). Genetic divergence and phylogenetic relationships of honeybee Apis mellifera (Hymenoptera: Apidae) populations from Greece and Cyprus using PCR-RFLP analysis of three mtDNA segments. Apidologie, 36: 335-344. doi: 10.1051/apido:2005021. Cornuet, J. M. & Garnery, L. (1991). Genetic Diversity in Apis mellifera. In: Smith, DR. Ed. Diversity in the genus Apis, Westview Press, Boulder, Co. De La Rúa, P., Simon, U.E., Tide, A.C., Moritz, R.F.A. & Fuchs, S. (2000). MtDNA variation in Apis cerana populations from the Philippines. Heredity, 84: 124-130. Estoup, A., Solignac, M., Harry, M. & Cornuet, J.M. (1993). Characterization of (GT)n and (CT)n microsatellites in two insect species: Apis mellifera and Bombus terrestris. Nucleic Acids Research, 21: 1427-1431. Garnery, L., Cornuet, J.M. & Solignac, M. (1992). Evolutionary history of the honeybee (Apis mellifera L.) inferred from mitochondrial DNA analysis. Molecular Ecology, 3: 145-154. Garnery, L., Solignac, M., Celebrano, G. & Cornuet, J.M. (1993). A simple test using restricted PCR-amplified mitochondrial DNA to study the genetic structure of Apis mellifera L. Experientia, 49: 1016-1021. Garnery, L., Franck, P., Baudry, E., Vautrin, D. & Cornuet, J.M. (1998). Genetic biodiversity of the West european honey bee (Apis mellifera and A.m. iberica) I: Mitochondrial DNA. Genet Selection Evolution, 30: S31-S47. Francýsco, F. O., Silvestre, D. & Arias, M. C. (2001). Mitochonrial DNA characterization of five species of Plebeia (Apidae: Meliponini): RFLP and resitriction maps. Apidologie, 32: 323-332. Franck, P., Garnery, L., Solignac, M. & Cornuet, J. M. ( 2000). Molecular confirmation of a fourth lineage in honeybees from the Near East. Apidologie, 31: 167-180. doi.org/10.1051/ apido:2000114. Hall, H.G. (1990): Parental analysis of introgressive hybridization between African and European honeybees using nuclear DNA RFLPs. Genetics, 125: 611-621. Hunt, J.G. & Page Jr, E.R. (1992). Patterns of inheritance with RAPD molecular markers reveal novel types of polymorphism in the honey bee. Theoretical and Applied Genetics, 85: 15-20. Kandemir, I., Kence, M., Sheppard, W. S. & Kence, A. (2006). Mitochondrial DNA variation in honeybee (Apis mellifera L.) population from Turkey. Journal of Apicultural Research, 45: 33-38. doi: 10.3896/IBRA.1.45.1.08. Kamrani, B., Pirany, N., Hashemi, A. & Kamrani, M. (2012). Genetic Characterization of Honey bees (Hymenoptera: Apidea) Populations from North West of Iran Using RAPD Sociobiology 65(3): 482-490 (September, 2018) 489 Markers. Technical Journal of Engineering and Applied Sciences, 2: 430-435. Kekecoglu, M., Bouga, M., Soysal, M.I. & Harizanis, P. (2009). Genetic divergence and phylogenetic relationships of honey bee populations from Turkey using PCR-RFLP’s analysis of two mtDNA segments. Bulgarian Journal of Agricultural Science, 15: 589-597. Meixner, M.D., Leta, M.A., Koeniger, N., Fuchs, S. (2011). The honey bees of Ethiopia represent a new subspecies of Apis mellifera – Apis mellifera simensis ssp. Apidologie, 42: 425-437. doi: 10.1007/s13592-011-0007-y. Meusel, M. S & Moritz, R. F. A. (1993). Transfer of paternal mitochondrial DNA in fertilization of honeybees (Apis mellifera L.) eggs. Current Genetic, 24: 539-543. Moradi, M. & Kandemir, A. (2004). Morphometric and Allozyme Variability in Persian Bee Population from the Alburz Mountains, Iran. Iranian Journal of Science and Technology, 5: 151-166. Moritz, R.F.A., Hawkins, C.F., Crozier, R.H. & McKinley, A.G. (1986). A mitochondrial DNA polymorphism in honey bees (Apis mellifera L.). Experientia, 42: 322-324. Mcelroy, D., Moran, P., Bermingham, E. & Kornfield, I. (1991). Reap: The Restriction Enzyme Analysis Package, Version 4.0. University Of Maine, Orono. Nash, J.H.E. (199). DNAfrag, program version 3.03, Institute for Biological Sciences National Research Council of Canada, Ottawa, Ontario, Canada. Nei, M. & Tajima, F. (1981). DNA polymorphism detectable by restriction endonucleases. Genetics, 97: 145-163. Nei, M. & Miller, J. C. (1990). A simple method for estimating average number of nucleotide substitutions within and between populations from restriction data. Genetics, 125: 873-879. Nunamaker, R. A., Wilson, W. T. & Haley, B. E. (1984). Electrophoretic detection of Africanized honey bee (Apis mellifera scutellata) in Guatemala and Mexico based on malate dehydrogenase allozyme patterns. Journal of the Kansas Entomological Society, 57: 622-631. Ozdil, F., Aytekin, I., Ilhan, F. & Boztepe, S. (2012). Genetic variation in Turkish honeybees Apis mellifera anatoliaca, A. m. caucasica, A. m. meda (Hymenoptera: Apidae) inferred from RFLP analysis of three mtDNA regions (16S rDNA- COI-ND5). European Journal of Entomology, 109: 161-167. Oldroyd, B.P., Cornuet, J.M., Rowe, D., Rinderer, E.T. & Crozier, R.H. (1995). Racial admixture of Apis mellifera in Tasmania, Australia: similarities and differences with natural hybrid zones in Europe. Heredity, 74: 315-325. Palmer, M. N., Smith, D. R. & Kaftanoglu, O. (2000). Turkish Honeybees: Genetic variation and evidence for a fourth lineage of Apis mellifera mtDNA. Journal of Heredity, 91: 42-46. Pedersen, B.V. (1996). On the phylogenetic position of the Danish strain of the black honeybee (the Laeso bee), Apis mellifera mellifera L. (Hymenoptera: Apidae) inferred from mitochondrial DNA sequences. Entomologica Scandinavica, 27: 241-250. Pinto, M. A., Johnston, J. S., Rubink, W. L., Coulson, R. N., Patton, J. C. & Sheppard, W. S. (2003). Identification of Africanized honey bee (Hymenoptera: Aphidae) mitochondrial DNA: Validation of a Rapid Polymerase Chain Reaction- Based Assay. Annals of the Entomological Society of America, 96: 679-684. Rahimi, A. & Asadi, M. (2010). Morphological characteristics of Apis mellifera meda (Hymenoptera: Apidae) in Saghez (west of Iran). Nature Montenegro, 10: 101-107. Rahimi, A. & Mirmoayedi, A. (2013). Evaluation of morphlogical characteristics of honey bee Apis mellifera meda (Hymenoptera: Apidae) in Mazandaran (North of Iran). Technical Journal of Engineering and Applied Sciences, 3: 1280-1284. Rahimi, A., Asadi, M., Abdolshahi, R. (2014a) A. Genetic diversity of honey bee (Apis mellifera meda) populations using microsatellite markers in Jiroft. Journal of Science and Tactic of Honey Bee, 6: 26-33. Rahimi, A., Miromayedi, A., Kahrizi, D., Abdolshahi, R., Kazemi, E., Yari KH. (2014b) B. Microsatellite genetic diversity of Apis mellifera meda skorikov. Molecular Biology Reports, 41: 7755-7761. doi: 10.1007/s11033-014-3667-7. Rahimi, A., Hasheminasab, H., Azati, N. (2015a). Predicting honey production based on morphological characteristics of honey bee (Apis mellifera L.) using multiple regression model. Ecology-Environment-Conservation, 21: 29-33. Rahimi, A (2015b). Study of the genetic diversity of Iranian honey bee (Apis mellifera meda Skorikow, 1829) populations using the mtDNA COI–COII intergenic region. Biologija, 61: 54-59. Rahimi, A., Mirmoayedi, A., Kahrizi, D., zaraei, L. & Jamali, S. (2016). Genetic diversity of Iranian honey bee (Apis mellifera meda Skorikow, 1829) populations based on ISSR markers. Cellular and Molecular Biology, 62 (4): 53-58. Rahimi, A., Mirmoayedi, A., Kahrizi, D., zaraei, L. & Jamali, S. (2017). Morphometric diversity and phylogenetic relationships among Iranian honey bee (Apis mellifera meda Skorikow, 1829) populations using morphological characters. Sociobiology, 64: 33-41. doi: 10.13102/sociobiology.v64i1.1179. Royan, M., Rahimi, G., Esmaeilkhanian, S. & Ansari, Z. (2007). A study on the genetic diversity of the Apis mellifera meda population in the south coast of the Caspian Sea using microsatellite markers. Journal of Apicultural Research and Bee World, 46: 236-241. doi: 10.3896/IBRA.1.46.4.05. A Rahimi, A Mirmoayedi, D Kahrizi, L Zarei, S Jamali – Genetic variation in Iranian honey bee490 Ruttner, F. (1988). Biogeography and Taxonomy of Honeybees, Springer-Verlag, Berlin, 284 p. Ruttner, F. (1992). Naturgeschichte der Honigbienen. Ehrenwirth Verlag. M¸nich. Germany. 357pp. Sheppard, W.S., Arias, M.C., Grech, A. & Meixner, M.D. (1997). Apis mellifera ruttneri, a new honey bee subspecies from Malta. Apidologie, 28: 287-293. Sheppard, W.S. & Meixner, M.D. (2003). Apis mellifera pomonella, a new honey bee subspecies from Central Asia. Apidologie, 34: 367-375. doi: 10.1051/apido:2003037 Smith, D. R. (1991). Mitochondrial DNA and honey bee biogeography. In: Smith, DR. (ed) Diversity in the genus Apis. Boulder, CO Westview, pp. 131-176. Smith, D.R., Taylor, O.R. & Brown, W.M. (1989). Neotropical Africanized honeybees have African mitochondrial DNA. Nature, 339: 213-215. Tarès, S., Cornuet, J.M. & Abad P. (1993). Characterization of an Unusually Conserved Alu I Highly Reiterated DNA Sequence Family From the Honeybee Apis mellifera. Genetics, 134: 1195-1204.