Genetic Diversity of Meliponamandacaia SMITH 1863 ( Hymenoptera , Apidae ) , an Endemic Bee Species from Brazilian Caatinga , Using ISSR

1 Departamento de Ciências Biológicas, Universidade Estadual do Sudoeste da Bahia, Avenida José Moreira Sobrinho s/n, 45206-190 Jequié, BA, Brazil 2 Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, USP, Rua do Matão 277, Cidade Universitária, 05508-090 São Paulo, SP, Brazil 3 Escola Agrotécnica Federal de Catu, Rua Barão de Camacari 118, 48110-000 Catu, BA, Brazil 4 Departamento de Biologia Geral, Universidade Federal de Viçosa, Avenida Peter Henry Rolfs s/n Campus Universitário, 36571-000 Viçosa, MG, Brazil


Introduction
The bee species Melipona mandacaia SMITH 1863 belongs to the subtribe Meliponina, which comprises the stingless bees.This species is endemic to Caatinga biome, being widespread in the Brazilian states of Piauí, Ceará, Bahia, Paraíba, and Pernambuco, usually close to São Francisco River [1,2].It plays an important role in Caatinga, acting as specific pollinators of this biome, besides presenting a great potential in meliponiculture [3].
The Caatinga is an exclusively Brazilian biome, comprising a wide drought belt in South America.This biome encompasses about 800,000 km 2 (8.6% of Brazilian territory) being surrounded by Atlantic Forest to the east, Amazon Forest to the west, and Brazilian savannah (Cerrado) to the south [4].
The genetic resources from Caatinga suffer accentuated pressure by continuous deforestation and expansion of agricultural frontiers.Consequently, the number of trees used for nesting of stingless bees is becoming more and more scarce [5], leading to population decreases or even local extinction of some bee species [6].Studies in Caatinga bees have shown several peculiarities, such as endemism and specific interactions between bees and local flora [7].Therefore, knowledge about this ecosystem and its potentialities, as well as the conservation of original covers, is essential.However, little is known about the genetic population structure of endemic species from Caatinga, and the conservation and the importance of its biota to local biodiversity remain overlooked [8].Studies of genetic diversity in natural populations usually refer to quantification of levels of variation within and among populations [9].When combined to other biological aspects such as reproductive behavior and dispersion process, these studies are able to provide valuable insights for defining conservation programs of a given species [10].
Amongst the several molecular markers available for genetic studies, the Inter Single Sequence Repeats (ISSRs) stands out as an efficient technique in the genome characterization of plants, fungi, vertebrates, and insects [11].This methodology allows detecting polymorphism in DNA regions flanked by microsatellites without requiring previous isolation and sequencing of specific DNA fragments.
Moreover, the estimates of genetic diversity of M. mandacaia might contribute for defining further management and conservation strategies.Therefore, the goal of the present study was to estimate the genetic diversity and structure of M. mandacaia throughout Bahia state, focusing on three main questions.(i) How much genetic diversity is there in M. mandacaia?(ii) How this diversity is structured?(iii) Is the genetic variation correlated to geographic distribution?

Material and Methods
2.1.Sampling and DNA Extraction.Samples of 104 colonies of M. mandacaia were collected in 12 localities in Caatinga along Bahia state between 2003 and 2010 (Figure 1 and Table 1).The collected specimens were stored in absolute ethanol and kept at −20 • C prior analyses.The total DNA was extracted from one worker bee per colony, based on the methodology proposed by Waldschmidt et al. [12].The DNA samples were quantified in 0.8% agarose gel to verify their concentration and integrity.

ISSR-PCR.
Firstly, 70 ISSR primers (manufactured by UBC) were tested and 15 were selected based on their reproducibility, definition, and number of bands.Afterwards, the effects of concentration of primer (0.05; 0.10, and 0.15 μM), DNA (10 and 50 ng), and annealing temperature (48 to 60 • C) were analyzed.After optimizing the PCRs, 10 primers were selected for population analyses (Table 2).The amplification reactions were adjusted based on the methodology proposed by Eiadthong et al. [13].
The PCRs comprised a final volume of 25 μL including 10 ng of template DNA, 2.5 μL of 10x buffer (Biotools), 2.0 μL of dNTPs at 200 μM each, 0.5 U of Taq polymerase (Biotools) and 0.4 pmoles of primer.The PCR conditions involved an Psyche 3

Data Analysis.
The amplification products were transformed into binary data according to the presence (1) or absence (0) of bands.The genetic diversity (H e ) [14] and the percentage of polymorphic loci were estimated using the software TFPGA v1.3 [15].The analysis of molecular variance (AMOVA) [16] and the estimate of structuring index among localities (Φ ST ) were performed using the software Arlequin v3.5.1.2[17].The test significance was based on 1,000 permutations.
Genetic diversity and population structure were also estimated using the Bayesian approach [18] available in the software HICKORY 1.1 [19].The H B and θ B indexes, analog with H e and Φ ST , respectively, were estimated as well.The analyses were carried out using four a priori models: full model, f = 0 model, θ = 0 model, and f free model.The best model was determined based on the deviance information criterion (DIC) according to software's authors, in which the lowest DIC values indicate the best adjusted model.A total of 100,000 MCMC (Markov Chain Monte Carlo) generations and a burn-in of 20% were implemented.The analyses were run twice in order to check the convergence of parameters.
Mantel's test was also performed to check the correlation between geographic distance and Φ ST values amongst localities.This test was applied in order to verify whether distance isolation between localities was present or not.The analyses were carried out using the software Arlequin v3.5.1.2[16] and the statistical significance was based on 1,000 permutations.

Results
The 10 selected ISSR primers yielded 109 bands with a mean polymorphism of 72.47% (Table 2).The primers UBC-889 and UBC-835 presented the highest number of bands (14) while UBC-841 had the lowest number (8).The primer UBC 807 presented the highest polymorphism (100%).The mean number of bands per primer was 10.9.The direct estimate of gene diversity (H e ) was equal to 0.2616.
The genetic structure based on AMOVA partitioned in two hierarchical levels showed that 70.39% of genetic variation was within localities whereas 29.61% of variation was related to variation among localities.The Φ ST value was 0.2961 (P < 0.000001; Table 3).When three hierarchical levels were considered (assuming localities between left and right bank of the São Francisco River as groups in the third hierarchical level) in AMOVA, no genetic variation among groups was revealed, resulting in 70.43% of variation observed among localities within groups and 29.67% within localities (Table 3).The Φ ST obtained using three hierarchical levels was equal to 0.29568.The pairwise Φ ST estimate among localities showed sites with moderate genetic structure and others with extremely low Φ ST values (Table 4).
When genetic diversity was estimated using Bayesian analysis, the full model presented the best adjustment based on DIC (Table 5).The estimates of genetic diversity (H B ) and genetic structure (θ B ) in this analysis were 0.2573 (highest posterior density (HPD) 97.5%: lower 0.2397, upper 0.2759) and 0.3289 (HPD 97.5%: lower 0.2746, upper 0.3809), respectively.
Mantel's test showed a positive correlation between genetic and geographic distances (r = 0.4542; P < 0.01), indicating isolation by distance among localities.

Discussion
The significant percentage of polymorphism and the mean number of bands obtained were 72.47% and 10.9, respectively.These values differ from previous reports in other species of Hymenoptera.Nascimento et al. [20]   with 59.18% of polymorphism in M. quadrifasciata from Minas Gerais state.Paplauskiené et al. [21] observed a mean number of six bands per primer using 11 ISSR primers for a subspecies of Apis mellifera.Borba et al. [22], studying lineages of Trichogramma, a group of small parasitic wasps, using 11 ISSR primers reported a mean number of 16 bands per primer and a high polymorphism (96%).The low value of genetic diversity (H e = 0.2616 and H B = 0.2573) herein observed has been reported in other studies based on molecular markers in stingless bees.Tavares et al. [23], using RAPD markers, observed an H e = 0.21 for M. mondury and H e = 0.23 for M. rufiventris.Borges et al. [24], based on studies of microsatellite markers in Partamona helleri that reported a genetic diversity (H e ) of 0.18.Nascimento et al. [20], analyzing the genetic variability of M. quadrifasciata with ISSR markers, reported a genetic diversity of H e = 0.20.Such reduced values of genetic diversity might be related to the swarming behavior for formation of new nests.In this process, the new nest can only be funded within short distances since it will depend on both food and workers from the mother colony up to its full establishment, which favors endogamy [25].
The two-hierarchical level AMOVA revealed a higher percentage of variation within localities (70.39%).The high concentration of genetic variability within localities is usually associated with lack or restrictions to gene flow among individuals from different localities, which might potentially lead to increased inbreeding within localities.When we partitioned the AMOVA into three hierarchical levels, placing São Francisco River as a geographic barrier among groups, no variation was observed, while the variation within localities remained high (70.43%).
The variation percentage detected by AMOVA and both Φ ST (0.2961) and θ B (0.3289) values showed that M. mandacaia is moderately structured through sampled localities when compared to other studies in stingless bees, where higher Φ ST indexes suggest high genetic structuring ( [26], Φ ST = 0.90; [27], Φ ST = 0.76 and 0.77; [20], Φ ST = 0.59), even though these studies were based on distinct molecular markers.Moreover, the pairwise Φ ST matrix among localities also confirms the genetic differentiation of M. mandacaia in Bahia state, revealing localities with high values of genetic structure, close to those reported by Nascimento et al. [20], while others presented low levels of structure (Table 5).
Mantel's test (r = 0.4542; P < 0.01) indicated a positive correlation among geographic and genetic distances, revealing an isolation by distance among sampled localities.Nascimento et al. [20] reported r = 0.3998 (P < 0.05) and also inferred that isolation by distance was present in M. quadrifasciata.This evidence corroborates the results of genetic structure obtained by AMOVA, Φ ST , and Bayesian analysis.
Mantel's test values and the pairwise Φ ST matrix among localities showed that the longer the geographic distance between them, the higher the genetic differentiation is.According to the Φ ST matrix (Table 4), three groups can be distinguished: group 1 (1 (Paulo Afonso), 2 (Macururé), 3 (Uauá) and 4 (Juazeiro)); group 2 (5 (Lapão), 6 (Irecê), 7 (São Gabriel), 8 (Central) and 9 (Uibaí)); group 3 (10 (Morpará), 11 (Muquém do S. Francisco) and 12 (Serra do Ramalho)).The second group presented the lowest genetic differentiation among localities due to the increased geographical distance, when compared to groups 1 and 3, it presents intermediary Φ ST values, as expected according to geographical distances.The groups 1 and 3 presented high levels of genetic differentiation when compared to group 2. However, the differentiation within groups 1 and 3 is as high as the differentiation between groups.This result might be related either to a preexisting differentiation among analyzed samples or to the geographic distance among the localities in this group.
The moderate genetic structure observed in this work could be explained due to the short period of separation among analyzed localities or else to a restricted gene flow.The restriction or lack of gene flow might be a consequence historical habitat fragmentation of biome Caatinga.Phylogeographic analyses using DNA sequences can help elucidate this scenario.
The environmental degradation has threatened bee populations worldwide [28], once woody and large trees, one of the main resources for nesting of bees, have decreased [29].In spite of its continuous deforestation, there are only 11 protection areas (including national parks, ecological stations, and biological reserves) in Caatinga that, together, represent less than 1% of the original biome [30].Taking into account the relevance of M. mandacaia as pollinators to the maintenance of Caatinga and the increased environmental impacts on this biome, our data suggest that further genetic studies, such as phylogeographic ones, are necessary to depict how the shape of genetic variation in this species across its geographic distribution is organized.

Table 1 :
Sampled localities of Melipona mandacaia, number of colonies per locality (N), altitude, and geographic coordinates.

Table 2 :
Selected ISSR primers with their respective sequences, number of bands, and percentage of polymorphism.
• C for 1 minute, 53 • C for 2 minutes, and 72 • C for 2 minutes; plus a final extension step at 72 • C for 7 minutes.All reactions included a negative control with all PCR components but without DNA.The PCR products were separated by electrophoresis in 1.2% agarose gel with 0.2 μg/ml of ethidium bromide, immersed in 1X TBE buffer (Tris-Borate 90 mM, EDTA 1 mM, and pH 8.0).The DNA fragments (bands) were visualized under ultraviolet light and photodocumented.

Table 3 :
Analysis of molecular variance (AMOVA) with two hierarchical levels and with three hierarchical levels, testing the São Francisco River as a geographic barrier in M. mandacaia.
using 11 ISSR primers observed 13.36 bands per primer, on average, Psyche

Table 4 :
Matrix of Φ ST values for each pairwise combination among specimens from 12 localities based on 109 ISSR loci.
The pairs with highest structuring levels are shown in bold.

Table 5 :
Model comparison of posterior means for Bayesian H B and θ B by Hickory.Comma values for H B and θ B represent lower and upper limits for 97.5% of high posterior density.Best fit model chosen based on DIC is indicated in bold.