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 16 S rDNA )

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

Iran is a vast country comprising a land area about 1648195 km 2 .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.

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 electrophoresed 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 electrophoretically on 2% or 3% agarose gels in 1× TBE buffer, stained with 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   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  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., 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.

Fig 1 .
Fig 1. Geographical locations of the 20 sampled honey bee populations in Iran.N is the number of sampled bees per population.

Fig 2 .
Fig 2. Cluster analysis of 300 honey bee individuals from twenty provinces of Iran based on UPGMA method.
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 Genetic distance (D) values between pairs of the Iranian honey bee populations.

Table 1 .
Sampling localities, geographical positions and the number of honey bees used for PCR-RFLP analysis in Iranian honey bee.

Table 2 .
Restriction fragment patterns generated from analysis of COI and 16S rDNA gene segments of Iranian honey bee.

Table 3 .
Composite genotypes (haplotypes), haplotype diversity and sample size of all the populations studied.

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.