Journal of Hymenoptera Research 38: 1–9, doi: 10.3897/JHR.38.7302
New molecular evidence for fragmentation between two distant populations of the threatened stingless bee Melipona subnitida Ducke (Hymenoptera, Apidae, Meliponini)
Geice R. Silva 1, Bruno A. Souza 2, Fábia M. Pereira 2, Maria T. R. Lopes 2, Sérgio E. S. Valente 1, Fábio M. Diniz 2
1 Department of Biological Sciences, Federal University of Piauí, CEP: 64049-550, Teresina, PI, Brazil
2 EMBRAPA Meio-Norte, Laboratory of Molecular Biology & Biotechnology, CP: 01, CEP: 64.006-220, Teresina, PI, Brazil

Corresponding author: Fábio M. Diniz (

Academic editor: Jack Neff

received 18 February 2014 | accepted 3 March 2014 | Published 12 June 2014
(C) 2014 Geice R. Silva. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
For reference, use of the paginated PDF or printed version of this article is recommended.

Citation: Silva GR, Souza BA, Pereira FM, Lopes MTR, Valente SES, Diniz FM (2014) New molecular evidence for fragmentation between two distant populations of the threatened stingless bee Melipona subnitida Ducke (Hymenoptera, Apidae, Meliponini). Journal of Hymenoptera Research 38: 1–9. doi: 10.3897/JHR.38.7302


For a snapshot assessment of the genetic diversity present within Melipona subnitida, an endemic stingless bee distributed in the semi-arid region of northeastern Brazil, populations separated by over 1, 000 km distance were analyzed by ISSR genotyping. This is a prerequisite for the establishment of efficient management and conservation practices. From 21 ISSR primers tested, only nine revealed consistent and polymorphic bands (loci). PCR reactions resulted in 165 loci, of which 92 were polymorphic (57.5%). Both ΦST (ARLEQUIN) and θB (HICKORY) presented high values of similar magnitude (0.34, p<0.0001 and 0.33, p<0.0001, respectively), showing that these two groups were highly structured. The dendrogram obtained by the cluster analysis and the scatter-plot of the PCoA corroborate with the data presented by the AMOVA and θB tests. Clear evidence of subdivision among sampling sites was also observed by the Bayesian grouping model analysis (STRUCTURE) of the ISSR data. It is clear from this study that conservation strategies should take into account the heterogeneity of these two separate populations, and address actions towards their sustainability by integrating our findings with ecological tools.


Population differentiation, Hymenoptera, Jandaíra, genetic diversity, ISSR markers


Melipona subnitida Ducke (Hymenoptera, Apidae, Meliponini) is an endemic stingless bee distributed in the Caatinga, the semi-arid region of northeastern Brazil (Michener 2007, Camargo and Pedro 2013). The species, called Jandaíra by natives, has great ecological significance as a pollinator of the local native flora and cultivated crops and is of economic importance in honey production, which is valued for its alleged medicinal properties and antibacterial activity (Cruz et al. 2004; Alves et al. 2008; Silva and Paz 2012).

The stingless bees are currently threatened by the increasing destruction of native semi-arid vegetation and by the intensification of agriculture in the Caatinga (Roulston and Goodell 2011; Pereira et al. 2011). In response, small populations of stingless bees may gradually decline, resulting in local extinction. The assessment of the genetic diversity present within Melipona subnitida populations is therefore a prerequisite for the establishment of efficient management and conservation practices.

Molecular markers have proven to be decisive in elucidating diversity and differences at the DNA level in microorganisms, plants and animals (Panwar et al. 2010; Sebastien et al. 2012; Rana et al. 2012; Bonatti et al. 2014). Among several markers, Inter Simple Sequence Repeats (ISSR) have been used due to their low cost and high level of polymorphism, and as an alternative to overcome reproducibility problems commonly found in other PCR-based markers (Abbot 2001; Lima-Brito et al. 2011). ISSRs are semiarbitrary markers based on DNA amplification by PCR in the presence of single 15- to 20-bp long primer complementary to a target short sequence repeat (Zietkiewicz et al. 1994). Despite the existence of few genetic studies related to stingless bees, it has recently been shown that ISSR markers can be useful in the analysis of natural bee populations, contributing to the development of management strategies of these important genetic resources (Nascimento et al. 2010; Tavares et al. 2013).

The present study uses ISSR analysis to investigate the degree of genetic differentiation between two Melipona subnitida populations separated by over 1, 000 km in the Brazilian Caatinga.

Worker bees were randomly collected from natural colonies (one bee from each of 30 colonies) distributed in 2 locations only: (1) 15 nests in Natal (NAT; coordinates: 5°48'04"S, 35°11'08"W; State of Rio Grande do Norte) and (2) 15 nests in Ilha das Canárias, Parnaíba River Delta (PAR; coordinates: 2°46'39"S, 41°51'59"W; on the border of the states of Piauí and Maranhão) in Brazil (Figure 1A). All the samples were taken to the laboratory and stored at –20 °C until further use. Total genomic DNA was extracted from each adult worker thorax using a phenol/chloroformalcohol isoamyl (25:24:1, v:v:v) extraction of SDS/proteinase-K digested tissue of each individual (Sambrook et al. 1989). High molecular weight DNA was isolated by ethanol precipitation and visualized by agarose gel electrophoresis.

Figure 1.

A Sampling sites of Melipona subnitida: NAT (coordinates: 5°48'04"S, 35°11'08"W; State of Rio Grande do Norte) and PAR (coordinates: 2°46'39"S, 41°51'59"W; on the border of the States of Piauí and Maranhão) B Clustering analysis using UPGMA for Melipona subnitida genotypes included in this study based on DICE similarity coefficient values. Numbers indicate bootstrap values for nodes retained by more than 50% of bootstrap replicates (1000 replications) C Scatter-plot of the principal coordinate analysis (PCoA) using ISSR loci. ■ PAR genotypes; ● NAT genotypes D Bar plot from Inferred population structure of using the Bayesian grouping admixture model-based program STRUCTURE (K = 2).

The extracted DNA was then amplified by PCR using twenty-one primers developed by the University of British Columbia (primer set #9). PCRs were carried out in 20 µL-reaction volumes containing approximately 10–50 ng of DNA, 1× PCR buffer (40 mM Tris-HCl; 100 mM KCl), 1.5 µM of primer, 6.25 mM MgCl2, 1.5 µM of each dNTP, 1.25 U of Invitrogen Taq DNA polymerase. All amplifications were carried out on a VERITITM Gradient Thermalcycler (APPLIED BIOSYSTEMS). The following PCR conditions were used: an initial denaturation at 94 °C for 1.5 min, followed by 45 cycles of 94 °C/40 s, 36 °C/1 min. and 72 °C/2 min., and a final extension of 72 °C/5 min. ISSR markers were screened by silver nitrate detection on denatured 6% polyacrylamide gels, which were scored for absence (0) and presence (1) of bands across genotypes and entered into a binary matrix.

Sample polymorphism was estimated as the percentage of polymorphic bands (PPB) in the total number of bands. The program HICKORY v.1.1, which implements the Bayesian method (Holsinger et al. 2002), was used for estimating θB (FST analogue), heterozygosity (HS), and the inbreeding coefficient (f), an FIS analogue for dominant markers. Analysis of Molecular Variance (AMOVA) was conducted using all the amplified loci to check the occurrence of population structure among sampling localities using ARLEQUIN v.3.11 (Excoffier et al. 2006). The use of different algorithms for the calculation of FST analogues was an additional effort to check for the reliability of the data presented by ISSR markers.

Further, a Bayesian grouping admixture model (burn-in length 100, 000 iteractions; run length 100, 000; K = 1 to 8 subpopulations tested) was used to infer the number of subpopulations (software STRUCTURE 2.3.3; Pritchard et al. 2000). These results were analyzed using STRUCTURE HARVESTERWeb v.0.6.9 (Evanno et al. 2005).

The similarity among samples within populations was estimated using PAST v1.34 (Hammer et al. 2001) according to Dice’s coefficient. Cluster analysis using the unweighted pair-group method arithmetic average (UPGMA) and multidimensional principal coordinate analysis (PCoA) were performed on the data set to reveal the degree of genetic differentiation between sites.

From 21 primers initially screened for their ability to generate ISSR loci, only nine revealed consistent and polymorphic bands (loci) with 30 Jandaíra worker bees. The other 12 ISSR markers were monomorphic or had unreliable amplification and therefore are not included in the genetic diversity analysis. Polymorphic ISSR primers were also considered reproducible after repeated PCRs, under the same reaction conditions and, therefore, selected for genotyping (Table 1). PCR reactions involving these nine primers resulted in 165 loci (bands) or 18.3 bands/primer, of which 92 were polymorphic (10.2 polymorphic bands/primer) ranging in size from 250 to 1636 bp, corresponding to an average polymorphism of 57.5%. Genotyping showed that most of the detected loci were polymorphic. Overall ISSR polymorphism in Melipona subnitida was similar to that of Melipona quadrifasciata Lepeletier (67%) (Nascimento et al. 2010). The number of polymorphic bands per primer ranged from 5 (UBC-884 and UBC-888) to 23 (UBC-899).

Table 1.

Primer names and sequences used in the ISSR analysis, number of polymorphic bands per primer and range of molecular weight in base pairs (bp) amplified by PCR-ISSR for 30 Melipona subnitida worker bees.

Primer code Primer sequence (5’–3’)* Total number of bands/loci Number of polymorphic loci Size range of bands (bp)
Total Overall PPL PAR Polymorphism (%) NAT Polymorphism (%)
UBC-827 ACA ACA ACA ACA ACA CG 15 6 0.400 3 50.0 6 100.0 500-1010
UBC-834 AGA GAG AGA GAG AGA GYT 25 8 0.320 7 87.5 2 25.0 300-1000
UBC-841 GAG AGA GAG AGA GAG AYC 21 16 0.762 11 68.8 13 81.3 500-1010
UBC-845 CTC TCT CTC TCT CTC TRG 9 9 1.000 8 88.9 8 88.9 506-1636
UBC-884 HBH AGA GAG AGA GAG AG 18 5 0.278 4 80.0 3 60.0 500-1010
UBC-886 VDV CTC TCT CTC TCT CT 16 11 0.688 6 54.5 11 100.0 500-1600
UBC-888 BDB CAC ACA CAC ACA CA 14 5 0.357 3 60.0 4 80.0 500-1010
UBC-899 CAT GGT GTT GGT CAT TGT TCC A 34 23 0.676 14 60.9 21 91.3 250-1636
UBC-852 TCT CTC TCT CTC TCT CRA 13 9 0.692 3 33.3 9 100.0 950-1636
Total 165 92 - 59 - 77 - -
Average - - 0.575 - 64.9 - 80.7 -

PPL: Proportion of polymorphic loci; *The following designations were used for degenerate sites: Y (C or T), R (A or G), H (A, T, or C), B (G, T, or C), V (C, G, or A) and D (G, A, or T). PAR: Ilha das Canarias, Parnaiba River Delta; NAT: Natal.

ISSR genotyping revealed differences in genetic diversity based on the percentage of polymorphic bands (% PPB), which was also estimated separately for each population. Result suggests that the NAT population (80.7%) is characterized by a higher genetic diversity than the PAR population (64.9%), which in theory might give the NAT population an increased ability to adapt to selective pressures.

The θB, f and HS values obtained from four different models of population structure using the Bayesian analysis are shown in Table 2. Of the models applied to the ISSR dataset, the full model, in which θ and f are estimated simultaneously, was preferred primarily because of its smaller deviant information criterion (DIC) value (657.47), with a difference of more than six units required to indicate strong support over all the other models (Holsinger et al. 2002). In the Bayesian approach θB (analogue to Wright’s FST), f (analogue to FIS), and HS (average panmictic heterozygosity across populations) were estimated to be 0.33, 0.31 and 0.29, respectively, indicating a pronounced genetic differentiation between populations, possibly caused by restricted gene flow and random genetic drift (Epperson 1995).

Table 2.

Summary of genetic variability, partitioning of diversity and limits for credible interval by Bayesian statistical procedures.

Models f θB HS DIC
Mean 2.5 % 97.5 % Mean 2.5 % 97.5 % Mean 2.5 % 97.5 %
Full 0.3096 0.0455 0.67 0.3275 0.2425 0.4188 0.2925 0.2565 0.3254 657.465
f = 0 0 - - 0.2737 0.2099 0.3451 0.3272 0.3122 0.3418 665.837
θB 0.3733 0.1702 0.6037 0 - - 0.3546 0.3319 0.3761 1107.250
f-free 0.4997 0.0224 0.9742 0.3578 0.2807 0.4389 0.2771 0.2467 0.3098 700.287

θB is analogous to Wright’s FST, f is analogous to FIS, HS is the average panmictic heterozygosity across populations, and DIC is deviance information criterion.

The analysis of molecular variance also provided additional support to the evidence of population differentiation in Melipona subnitida. The AMOVA analysis indicated that most of the genetic variation found in Melipona subnitida samples could be attributed to differences among individuals within populations (approximately 66% of the variance), but also a large part of the variation (34.35%) was due to differences among localities. The inbreeding coefficient (f = 0.31) provided by HICKORY agrees with the results obtained by the AMOVA, as moderate endogamy might be expected for strongly structured populations.

Both θB (HICKORY) and ΦST (ARLEQUIN) presented high values, but of similar magnitude and significance (0.33, p<0.0001 and 0.34, p<0.0001, respectively), showing that these two populations are highly structured (Tables 2–3). These two different approaches showed a general agreement among the results.

Table 3.

Hierarchical analysis of molecular variance (AMOVA) detected by the ISSR genotyping.

Source of variation df Sum of squares Variance components Percentage of total variance (%) P-value*
Among populations 1 96.333 5.697 34.35 < 0.0001
Among individuals within populations 28 304.800 10.886 65.65 0.014
TOTAL 29 401.133 16.583
Fixation index FST index = 0.3435

*p-values calculated from a random permutation test (10, 000 permutations).

A high degree of population differentiation has also been observed for Melipona quadrifasciata populations based on ISSR patterns (Nascimento et al. 2010) and Melipona rufiventris Lepeletier (Tavares et al. 2007) based on allozyme, microsatellite and random amplified polymorphic DNA (RAPD) molecular markers. Both studies were conducted within the State of Minas Gerais (Brazil). However, no correlation was found between the first internal transcribed spacer (ITS1) sequence divergence of Melipona subnitida populations and geographical distances in northeastern Brazil, which might be explained by the extremely high mutation rates of the ITS region in Melipona subnitida (Cruz et al. 2006).

Genetic differentiation within Melipona subnitida populations was probably because of low gene flow, caused by limited dispersal ability (Engels and Imperatriz-Fonseca 1998; Tavares et al. 2007), the large distance separating the NAT and PAR populations, and extensive disturbances of population dynamics due to anthropogenic habitat degradation and fragmentation (Quezada-Euan et al. 2007). More recently, mtDNA data has also pointed to genetic differentiation between these Melipona subnitida populations from Rio Grande do Norte (Mossoró city) and Piauí (Parnaíba city, on the border of Maranhão state) (Bonatti et al. 2014).

Furthermore, the dendrogram obtained by the cluster analysis (Figure 1B) and the scatter-plot of the PCoA (Figure 1C) revealed a clear separation of the species in two main clusters confirming a significant molecular genetic difference between the two populations. This topology corroborates the data presented by the Bayesian θB and AMOVA. A clear evidence of subdivision among sampling sites was also observed by the Bayesian grouping model analysis of the ISSR data (Figure 1D).

This study provides additional evidence that ISSR markers can be useful tools in defining population genetic substructuring in Melipona species. More importantly, the distinctiveness of populations in these two regions suggests that the NAT and PAR populations of Melipona subnitida have separate evolutionary histories. It is clear from this study that conservation strategies should take into account the heterogeneity of these two separate populations, and that actions should be addressed towards their sustainability by integrating our findings with ecological tools. Failing to do so would risk decimating the entire bee population by uncontrolled human activities in the region.


The authors are grateful to Mr. José Maria Vieira-Neto for help in field collection. This work was supported by grants from the Northeast Bank of Brazil (BNB/ETENE/FUNDECI) and Embrapa (Macroprograma 2, SEG Code:

Abbot P (2001) Individual and population variation in invertebrates revealed by Inter-simple Sequence Repeats (ISSRs). Journal of Insect Science 1.8, 3 pp.
Alves DFS, Cabral-Junior FC, Câmara-Cabral PPA, Oliveira-Junior RM, Rego ACM, Medeiros AC (2008) Effects of topical application of the honey of Melipona subnitida in infected wounds of rats. Revista Colégio Brasileiro de Cirurgiões 35: 88–193.
Bonatti V, Simões ZLP, Franco FF, Francou TM (2014) Evidence of at least two evolutionary lineages in Melipona subnitida (apidae, Meliponini) suggested by mtDNA variability and geometrc morphometris of forewings. Naturwissenschaten 101: 17–24. doi: 10.1007/s00114-013-1123-5
Camargo JMF, Pedro SEM (2013) Meliponini Lepeletier, 1836. In: Moure JS, Urban D, Melo GAR (Orgs) Catalogue of Bees (Hymenoptera, Apoidea) in the Neotropical Region - online version. [accessed Dec/20/2013]
Cruz DO, Freitas BM, Silva LA, Silva SEM, Bomfim IGA (2004) Adaptação e comportamento de pastejo da abelha Jandaíra (Melipona subnitida Ducke) em ambiente protegido. Acta Scientiarum Animal Sciences 26: 293–298.
Cruz DO, Jorge DMM, Pereira JOP, Torres DC, Soares CEA, Freitas BM, Grangeiro TB (2006) Intraspecific variation in the first internal transcribed spacer (ITS1) of the nuclear ribosomal DNA in Melipona subnitida (Hymenoptera, Apidae), an endemic stingless bee from northeastern Brazil. Apidologie 37: 376–386. doi: 10.1051/apido:2006003
Engels W, Imperatriz-Fonseca VL (1998) Caste development, reproductive strategies and control of fertility in honeybees and stingless bees. In: Engels W (Ed) Social Insects. Springer-Verlag, Berlin, Germany, 167–230.
Epperson BK (1995) Spatial distributions of genotypes under isolation by distance. Genetics 140: 1431–1440.
Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology 14: 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x
Excoffier L, Laval G, Schneider S (2006) Computational and molecular population genetics lab (CMPG). Institute of Zoology University of Berne Baltzerstrasse 6 3012 Bern, Switzerland.
Hammer O, Harper DAT, Ryan PD (2001) Past: Paleontological Statistics Software Package for Education and Data Analysis. Paleontology Electronic 4: 1–9.
Holsinger KE, Lewis PO, Dey DK (2002) A Bayesian approach to inferring population structure from dominant markers. Molecular Ecology 11: 1157–1164. doi: 10.1046/j.1365-294X.2002.01512.x
Lima-Brito J, Castro L, Coutinho J, Morais F, Gomes L, Guedes-Pinto H, Carvalho A (2011) Genetic variability in Sambucus nigra L. clones: a preliminary molecular approach. Journal of Genetics 90: 47–52.
Michener CD (2007) The Bees of The World, 2nd ed. Johns Hopkins Univ. Press, Baltimore, MD.
Nascimento MA, Batalha-Filho H, Waldschmidt AM, Tavares MG, Campos LAO, Salomão TMF (2010) Variation and genetic structure of Melipona quadrifasciata Lepeletier. Genetic and Molecular Biology 33: 394–397. doi: 10.1590/S1415-47572010005000052
Panwar P, Nath M, Yadav VK, Kumar A (2010) Comparative evaluation of genetic diversity using RAPD, SSR and cytochrome P450 gene based markers with respect to calcium content in finger millet (Eleusine coracana L. Gaertn.). Journal of Genetics 89: 121–133. doi: 10.1007/s12041-010-0052-8
Pereira DS, Menezes PR, Filho VB, Sousa AH, Maracajá PB (2011) Abelhas indígenas criadas no Rio Grande do Norte. Acta Veterinária Brasílica 5(1): 81–91.
Pritchard JK,  Stephens M,  Donnelly P (2000)  Inference of population structure using multilocus genotype data. Genetics 155: 945–959.
Quezada-Euan JJG, Paxton RJ, Palmer KA, Itza WJM, Tay WT, Oldroyd BP (2007) Morphological and molecular characters reveal differentiation in a Neotropical social bee, Melipona beecheii (Apidae: Meliponini). Apidologie 38: 1–12. doi: 10.1051/apido:2007006
Rana TS, Narzary D, Ohri D (2012) Molecular differentiation of Chenopodium album complex and some related species using ISSR profiles and ITS sequences. Gene 495: 29–35. doi: 10.1016/j.gene.2011.12.031
Roulston TH, Goodell K (2011) The role of resources and risks in regulating wild bee populations. Annual Review of Entomology 56: 293–312. doi: 10.1146/annurev-ento-120709-144802
Sambrook J, Fritsch EFE, Maniatis T (1989) Molecular Cloning: A Laboratory Manual, 2nd ed. Cold Spring Harbor Lab. Press, New York, USA.
Sebastien A, Gruber MAM, Lester PJ (2012) Prevalence and genetic diversity of three bacterial endosymbionts (Wolbachia, Arsenophonus, and Rhizobiales) associated with the invasive yellow crazy ant (Anoplolepis gracilipes). Insect Sociaux 59: 33–40. doi: 10.1007/s00040-011-0184-8
Silva WP, Paz JRL (2012) Abelhas sem ferrão: muito mais do que uma importância econômica. Natureza online 10: 146–152.
Tavares MG, Almeida BS, Passamani PZ, Paiva SR, Resende HC, Campos AO, Alves RMO, Waldschmidt AM (2013) Genetic variability and population structure in Melipona scutellaris (Hymenoptera: Apidae) from Bahia, Brazil, based on molecular markers. Apidologie 44: 720–728. doi: 10.1007/s13592-013-0220-y
Tavares MG, Dias LAS, Borges AA, Lopes DM, Busse AHP, Costa RG, Salomão TMF, Campos LAO (2007) Genetic divergence between populations of the stingless bee uruçu amarela (Melipona rufiventris group, Hymenoptera, Meliponini): Is there a new Melipona species in the Brazilian state of Minas Gerais? Genetics and Molecular Biology 30: 667–675. doi: 10.1590/S1415-47572007000400027
Zietkiewicz E, Rafalski A, Labuda D (1994) Genome fingerprinting by simple sequence repeat (SSR) anchored polymerase chain reaction amplification. Genomics 20: 176–183. doi: 10.1006/geno.1994.1151