Molecular genetic diversity in populations of the stingless bee Plebeia remota: A case study

Genetic diversity is a major component of the biological diversity of an ecosystem. The survival of a population may be seriously threatened if its genetic diversity values are low. In this work, we measured the genetic diversity of the stingless bee Plebeia remota based on molecular data obtained by analyzing 15 microsatellite loci and sequencing two mitochondrial genes. Population structure and genetic diversity differed depending on the molecular marker analyzed: microsatellites showed low population structure and moderate to high genetic diversity, while mitochondrial DNA (mtDNA) showed high population structure and low diversity in three populations. Queen philopatry and male dispersal behavior are discussed as the main reasons for these findings.

Most of the Brazilian tropical flora is pollinated by bees, especially by those belonging to the tribe Meliponini (stingless bees) (Kerr et al., 1996;Nogueira-Neto, 1997). The increase in habitat loss can lead to severe consequences for bee populations and species diversity (Foley et al., 2005;Brown and Paxton, 2009). Studies focusing on Meliponini general biology, including genetic diversity, are still scarce in the literature. Genetic data are essential for a better understanding of macro-and micro-evolutionary processes and patterns in organisms, and provide support for conservation and managing programs (Moritz, 2002;Frankham et al., 2004).
Recent studies have indicated a low genetic diversity in feral populations of Brazilian stingless bees (Costa et al., 2005;Arias et al., 2006;Tavares et al., 2007;Borges et al., 2010;Brito and Arias, 2010;Francisco and Arias, 2010). This low genetic diversity may have negative consequences for the long-term population survival rate and raises important questions related to conservation programs for the Meliponini.
An adequate understanding of how genetic diversity is distributed and maintained among stingless bee populations requires a consideration of behavioral components such as philopatry. Philopatry may restrict individual dispersion, leading to inbreeding, which consequently reduces heterozygosity; philopatry also increases the effects of genetic drift due to population subdivision and isolation.
However, studies in a variety of organisms have shown that if one gender is philopatric then the other one normally mediates gene flow through dispersion (Whitehead, 1998;Nyakaana and Arctander, 1999;Kappeler et al., 2002;Apio et al., 2010). This behavioral mechanism minimizes the negative effects of philopatry.
The stingless bee Plebeia remota occurs in Bolivia and southeastern and southern Brazil (Camargo and Pedro, 2012). This species generally builds its nests in tree cavities, with colonies of up to 5,000 bees (van Benthem et al., 1995); the workers are small (~0.5 cm in length) (Hilário et al., 2007). In a previous investigation, Francisco and Arias (2010) described low intrapopulation mitochondrial polymorphism for this species. In the present work, we reanalyzed most of those samples to measure nuclear genetic diversity based on the amplification of microsatellite loci with specific primers and also by sequencing two mitochondrial genes.
One worker bee from each of 65 nests was analyzed for nuclear and mitochondrial loci. The samples originated from four localities (referred to from here on as "populations"): Cunha in São Paulo state (n = 13), Curitiba (n = 6) and Prudentópolis (n = 34) in Paraná state and Blumenau (n = 12) in Santa Catarina state ( Figure 1). Total DNA was extracted using Chelex ® 100 (Bio-Rad) according to a protocol described by Walsh et al. (1991). All individuals were genotyped for 15 microsatellite loci (Francisco et al., 2011): Prem03, Prem07, Prem57, Prem58, Prem70, Prem75a, Prem78, Prem79, Prem81a, Prem82, Prem83, Prem84, Prem87, Prem93 and Prem94. Microsatellite amplification and visualization were done as described by Francisco et al. (2011). Allelic richness (A), observed and expected heterozygosities (H O and H E , respectively) from Hardy-Weinberg proportions, percentage of polymorphic loci and allele frequencies were calculated for each population using Genalex v.6.41 (Peakall and Smouse, 2006). Due to differences in sample size, rarefaction was applied to allelic richness (Ar) by using the program HP-Rare 1.0 (Kalinowski, 2005). Log likelihood ratio statistics for linkage disequilibrium were computed using Genepop v.4.1.4 (Rousset, 2008). The Bonferroni correction (Rice, 1989) was applied when multiple comparisons were done. Population structure was analyzed with the program Structure v.2.3.3 (Pritchard et al., 2000). The program was set up for 500,000 Markov chain Monte Carlo repetitions after an initial burn-in of 20,000 repetitions. The number of structured populations (K) was estimated based on 10 replications for each K (from 1 to 4). The estimate of the best K was calcu-lated as described by Evanno et al. (2005) using Structure Harvester v.0.6.92 (Earl and VonHoldt, 2012). The program Clumpp v.1.1.2 (Jakobsson and Rosenberg, 2007) was used to align the 10 repetitions of the best K. The program Distruct v.1.1 (Rosenberg, 2004) was used to graphically display the results produced by Clumpp. Population structure was also analyzed using the D est estimator (Jost, 2008) which was calculated for each population pair by the program SMOGD v.1.2.5 (Crawford, 2010).
Two mitochondrial genes, cytochrome c oxidase subunit I (COI) and cytochrome b (Cytb), were partially amplified by using the primers mtD06 + mtD09 (Simon et al., 1994) and mtD26 (Simon et al., 1994) + AMB16 (Arias et al., 2008), respectively. PCR assays were done with 1 mL of DNA, 1x PCR buffer, 200 mM of each dNTP, 3 mM of MgCl 2 , 0.8 mM of each primer, 1 M of betaine anhydrous (USB Corporation) and 1 U of Taq DNA polymerase (Invitrogen) in a final volume of 10 mL. The amplification conditions consisted of an initial denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 60 s, annealing at 42°C for 60 s and elongation at 64°C for 80 s, and a final elongation step at 64°C for 10 min. PCR products (2 mL aliquots) were analyzed by electrophoresis in 0.8% agarose gels stained with GelRed (Biotium) and visualized under UV light. About 100-200 ng of each product was purified with 0.5 mL of ExoSAP-IT(USB Corporation) and used for sequencing reactions according to the manufacturer's recommended protocols (BigDye Terminator v.3.1 Cycle sequencing kit, Applied Biosystems). The sam- Francisco et al. 119 ples were analyzed in an automatic sequencer ABI PRISM 3100 Genetic Analyzer (Applied Biosystems). DNA sequences were edited using the Geneious v.5.1.6 software package (Drummond et al., 2010). The alignment was done using the algorithm Muscle (Edgar, 2004) from Geneious, with a maximum number of eight iterations. DnaSP v.5.10.01 software (Librado and Rozas, 2009) was used to identify individual haplotypes and their frequencies. A haplotype network was generated with the software Network v.4.6.1.0. Exact tests between pairs of populations were done using Arlequin v.3.5.1.3 (Excoffier and Lischer, 2010). All microsatellite loci analyzed were polymorphic. The allele frequencies for each locus and for each population are included in the Supplementary material to this paper (Table S1). The intrapopulation genetic diversity indices ranged from moderate to high ( A total of 794 bp (415 from COI and 379 bp from Cytb) was obtained for all individuals (GenBank accession numbers JQ517144-JQ517273). Ten haplotypes were identified, most of which were exclusive to a specific population, except for one (h04) ( Table 2). Table 2 also shows the haplotype and nucleotide diversity indices; they were not correlated to sample size. Figure 2 shows the network built to represent the associations between haplotypes and the genetic differentiation among populations. The maximum number of nucleotide differences between two haplotypes was 14 (1.8%). Exact tests based on haplotype frequencies showed differentiation between all population pairs (all p < 0.0005).
Population structure and genetic diversity varied, depending on the molecular marker analyzed. The microsatellite data showed a low population structure and moderate to high genetic diversity, whereas the mtDNA data showed a high population structure and low diversity in three populations. The mtDNA data suggested an absence of female gene flow among the populations, and reinforced the philopatric behavior of queens and its strong influence on the genetic differentiation observed. These findings agree with previous data obtained by RFPL of mtDNA that also showed no gene flow through females (Francisco and Arias, 2010).
The level of intrapopulation nuclear diversity was moderate to high. The Prudentópolis population had the lowest genetic diversity index, which suggested genetic isolation. The nuclear data also indicated a low population structure among the Cunha, Curitiba, and Blumenau populations. An absent or low genetic structure can be attributed to homoplasy in microsatellite size but should be accompanied by a decrease in genetic variability (Estoup et al., 2002), which was not the case here.
Since female philopatry was detected in these populations, the absence of genetic structure can be explained by male dispersal. The few genetic studies of Meliponini male congregations have demonstrated the presence of males from distant areas, with more than 100 colonies acting as male donors (Paxton, 2000;Cameron et al., 2004;Kraus et al., 2008;Mueller et al., 2012). 120 Genetic diversity in Plebeia remota  Thus, the divergence between mitochondrial and nuclear data is a consequence of the reproductive behavior of P. remota. The low mtDNA variability indicates a low dispersal capability of females, i.e., queen philopatry. In contrast, the high nuclear genetic variability is maintained by male dispersal. This male behavior is crucial to avoid inbreeding and to keep the population genetically healthy. Since Meliponini species show queen philopatry (Nogueira-Neto, 1954;Engels and Imperatriz-Fonseca, 1990), we expect a similar genetic scenario in other species.

Acknowledgments
We thank Susy Coelho and Ana Carolina Lima Novelli for technical assistance and anonymous reviewers for their comments and suggestions on an earlier version of this manuscript. FOF was supported by scholarships (99/11190-6 and 08/08546-4) and grants (04/15801-0 and 10/50597-5) from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP). LRS was supported by a scholarship from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). This work was developed in the Research Center on Biodiversity and Computing (BioComp) of the Universidade de São Paulo (USP), supported by the USP Provost's Office for Research.