EUROPEAN JOURNAL OF ENTOMOLOGY EUROPEAN JOURNAL OF ENTOMOLOGY Genome-wide discovery and characterization of microsatellite markers from Melipona fasciculata (Hymenoptera: Apidae), cross-ampli ﬁ cation and a snapshot assessment of the genetic diversity in two stingless bee populations

. Brazilian native meliponines are currently threatened by increased human impacts. The assessment of their genetic variation by microsatellite DNA markers can assist in the conservation of populations and help in the planning and establishment of ef ﬁ cient management strategies. The purpose of this study was to develop the ﬁ rst set of microsatellite markers for Melipona fasciculata , selected from partial genome assembly of Illumina paired-end reads. Primer pairs were designed for each detected locus at their ﬂ anking regions. Bee samples were genotyped from two different populations of Northeastern Brazil for marker characterization and validation. A total of 17 microsatellite loci displayed polymorphism. Mean H E and H O heterozygosities were 0.453 and 0.536, respectively. PIC across all loci ranged from 0.108 to 0.714. A genetic diversity analysis revealed high values for population differentiation estimates ( F ST = 0.194, R ST = 0.230, and D est = 0.162) within the investigated region. PCoA and Bayesian clustering showed a separation of the species into two distinct clusters. These microsatellite markers have demonstrated strong potential for population-level genetic studies. Moreover, the preliminary analysis of the genetic diversity in M. fasciculata provides provisional evidence of signi ﬁ cant population differentiation between the two studied populations.


INTRODUCTION
Stingless bees (Hymenoptera: Apidae: Meliponini) are a diverse group of bees regarded for their great economic and ecological importance. For instance, beekeeping provides a sustainable source of income under a low-cost investment for smallholder farming communities, and these native bees provide an effi cient pollination service in both natural and agricultural systems (Garibaldi et al., 2013).
Currently, native meliponines are threatened by increased human impacts such as destruction of native vegetation and consequent landscape transformation, and indiscriminate use of pesticides for agriculture (Winfree et al., 2009;Potts et al., 2010;Silva et al., 2015). Anthropogenic disturbances or intervention may negatively affect the existence of small populations of native Eur. J. Entomol. 115: 614-619, 2018 doi: 10.14411/eje.2018.058 NOTE tion kit (Illumina, San Diego, CA, USA). DNA was tagged and fragmented by the Nextera XT transposome, followed by limitedcycle PCR amplifi cation, AMPure XP magnetic-bead purifi cation (Agencourt Bioscience Corporation, Beverly, MA, USA) and the Illumina Nextera XT bead-based normalization protocol. The DNA library was sequenced using a MiSeq Benchtop Sequencer (Illumina). Contigs were created from the resulting paired-end sequence data (reads) using CLC Genomics Workbench 7.0.4 (Qiagen, Carlsbad, CA, USA).

SSR loci search and primer design
All these contigs were subsequently added directly into Msatcommander 0.8.2 (Faircloth, 2008) for detection of possible microsatellite loci with at least four repeats, except for dinucleotides (six repeats), and designing of primer pairs for each detected locus at their fl anking regions. Long mononucleotide repeats were ignored for marker development. Primer design was performed with the Primer3 (Rozen & Skaletsky, 2000).

SSR-PCR amplifi cation for primer validation and genotyping
Genomic DNA from fi ve individuals, each from different colonies, were initially used to validate all designed primer pairs using polymerase chain reactions (PCRs). Reactions were performed in a 10-μL total volume containing at least 10 ng of genomic DNA, with 1.25 to 1.5 × buffer, 2 to 2.5 mM MgCl 2 , 10 mM dNTP mix, 0.25 mM of each primer and 0.25 units of Taq DNA polymerase (Thermo Fisher Scientifi c, Waltham, MA, USA) or HotStar Taq DNA Polymerase (Qiagen). All amplifi cations were run in a Veriti 96-well Thermal Cycler (Applied Biosystems, Foster City, CA, USA) using the PCR temperature profi le indicated in Table 1. We also tested the cross-species amplifi cation success of these loci with three samples from each of four additional species, Melipona marginata, M. subnitida, Nannotrigona testaceicornis, and Frieseomelitta varia, as the DNA extraction described above. All amplifi cation products, from source and non-source bee species, were screened by silver nitrate detection on denatured 6% polyacrylamide gels. Subsequently, additional M. fasciculata samples were genotyped from two different locations, Piauí (20 inds) and Maranhão (20 inds), each from different colonies, to obtain baseline allele frequency information.

Data analysis
The genotyped data were initially analyzed using Microchecker ver. 2.2.3 (van Oosterhout et al., 2004) to test for the presence of null alleles, large allele dropout and scoring errors by stuttering. Observed and expected heterozygosities (H O and H E ), the number of alleles (Na), and the polymorphic information content (PIC) were determined using Cervus ver. 3.0 (Kalinowski et al., 2007). Allelic richness (A R ) as a measure of the number of alleles per locus independent of sample size was calculated by Fstat ver. 2.9.3.2 (Goudet, 1995). Deviations from Hardy-Weinberg equilibrium (HWE) and tests for linkage disequilibrium were conducted using Genepop software (Raymond & Rousset, 1995). The Bonferroni correction was applied when multiple pairwise tests were performed to assess the signifi cance (P < 0.05).
The genetic diversity for each locus was evaluated by Arlequin ver. 3.5 (Excoffi er & Lischer, 2010), which determined the value of θ (F ST ) for the whole sample set.
A Bayesian grouping admixture model was used to infer possible population structuring using the software Structure ver. 2.3.3 (Pritchard et al., 2000). The program was set up for 1,000,000 Markov chain Monte Carlo repetitions after an initial burn-in of 500,000 steps. The estimate of the best K was calculated based on fi ve replications for each K (from 1 to 6) using Structure harvester ver. 0.6.92 (Earl & vonHoldt, 2012). The program Clumpp fl owering plants are available, and (ii) its production of honey and geopropolis with antioxidant potential (Dutra et al., 2014) and anti-infl ammatory effect (Liberio et al., 2011). The species has also an important place in a rapidly growing market with wide possibilities for economic exploitation, mainly by family farmers (Holanda et al., 2012;Venturieri et al., 2012;Alves, 2013;Gostinski, 2018). As such, beekeepers have been continuously trying to improve honey production by nest transportation and trading among themselves (Kerr et al., 1996;Jaffé et al., 2015). In fact, this unregulated and uncontrolled practice of transporting and trading colonies from one region to another affects the genetic structure of stingless bee populations by altering factors such as geographic distance between populations or physical barriers to gene fl ow (Nogueira et al., 2014;Byatt et al., 2016;Jaffé et al., 2016).
Microsatellites, stretches of short DNA sequences tandemly repeated, have become the markers of choice for high-resolution assessment of genetic variation and population structure studies, most importantly, due to their abundance throughout the eukaryote genome and their hypervariability (Goldstein & Schlötterer, 1999). Emerging technologies in DNA sequencing (i.e. next generation sequencing -NGS) have proven to be useful for identifying microsatellite loci from the large amounts of sequence data they generate with much less effort and low cost, therefore, challenging traditional approaches for their development (Zalapa et al., 2012;Souza et al., 2015). Microsatellite markers developed for a studied species can also be transferred among related species, of the same genus or different genera in the same family, which considerably reduces their development costs and shortens development times in a non-source species (Selkoe & Toonen, 2006).
In this paper, we describe the fi rst set of microsatellite markers developed for Melipona fasciculata, selected from partial genome assembly of Illumina paired-end reads, and test the crossamplifi cation of all microsatellite loci in four non-source species. A preliminary analysis of its genetic variation within a relatively short geographical range is also performed to characterize and validate these polymorphic markers.

Bee materials and genomic DNA isolation
Genomic DNA from Melipona fasciculata was extracted from each adult worker thorax (n = 50) using the Wizard Genomic DNA Purifi cation Kit (Promega, Madison, WI, USA) according to the manufacturer's instruction. Bees were collected from hives originally from the Northeast region of Brazil, in the states of Piauí (Elesbão Veloso city; 6°11´56.2˝S and 42°11´43.8˝W) and Maranhão (São Bento city; 2°42´30.6˝S and 44°50´18.9˝W). DNA was also isolated from four additional species, M. marginata, M. subnitida, Nannotrigona testaceicornis and Frieseomelitta varia (Table 1) for cross-species amplifi cation testing. Three worker bee samples from each species were collected from hives located at Embrapa Meio-Norte in Teresina, Piauí (5°02´22.8˝S and 42°48´12.9˝W). The extracted DNA samples were electrophoresed on 0.8% agarose gel to test for overall quantity and quality of the DNA yield.

Library preparation and next generation sequencing
A single individual with the highest quality DNA yield was selected for sequencing. DNA was quantifi ed using a PicoGreen protocol and was run using a Perkin Elmer Fusion DNA Quantifi er (Perkin Elmer, Waltham, MA, USA). An Illumina paired-end library was created using 1 ng of genomic DNA, following the standard protocol of the Illumina Nextera XT Library Prepara-ver.1.12 (Jakobsson & Rosenberg, 2007) was used to align the fi ve repetitions of the best K. The program Distruct ver. 1.1 (Noah, 2004) was used to graphically display the results produced by Clumpp. Population structure was also analyzed using principal coordinate analysis (PCoA), R ST , a measure of genetic differentiation analogous to F ST but incorporating allele size information, and the D est estimator of actual differentiation as implemented in Genalex v. 6.5 (Peakall & Smouse, 2012).

Sequence assembly, SSR identifi cation and primer design
The genomic library, which was previously loaded as 16% of a MiSeq Reagent Kit v2 300 cycle sequencing run, produced 2,669,884 reads, which were assembled into 47,087 contigs. The program msatcommander 0.8.2 (Faircloth, 2008) identifi ed 9,954 contigs (11.3% of total contigs) containing 11,869 microsatellite loci, being in the majority mononucleotide (6,444) and dinucleotide (3,225) repeats. For ease of imaging and scoring, we chose to examine only tri-(734) and tetranucleotide (574) loci. From these, 37 loci were chosen for primer designing and validation in M. fasciculata, based on the presence of long, uninterrupted repeat units (≥ 10 repeats), fl anking regions of at least 50 bp in length containing no more than four monobase repeats [(A) n , (G) n , (C) n , (T) n ], and amplicon length between 100 and 270 bp.

Polymorphism and validation of Melipona fasciculata SSR markers
The Micro-Checker analysis of the entire dataset revealed null alleles for loci Mfsc3 and Mfsc11, at lower frequencies than 0.2 (Table 2). Null allele frequencies below 0.2 are acceptable in most microsatellite data sets (Dakin & Avise, 2004). When the dataset was divided into two populations (Piauí and Maranhão) only locus Mfsc3 indicated the presence of nulls, which may be a possible cause of its deviation from HWE in the Piauí population, even after Bonferroni correction for multiple comparisons at the 5% signifi cance level (critical value for P > 0.0029). No loci showed signifi cant linkage disequilibrium after Bonferroni correction.
PCR products of expected size with clear and consistent bands were obtained from 18 primer pairs, out of 37 sets tested on 50 individual bees from the two surveyed populations in Northeastern Brazil ( Table 2). The proportion of markers that generated consistent amplicons within their expected sizes was 48.6%. Expected product sizes for each microsatellite locus were based on sequence data from the partial genome assembly process. Seventeen microsatellite loci were polymorphic across the entire  M. subnitida, Nannotrigona testaceicornis, and Frieseomelitta varia). However, heterologous amplifi cation proved more successful with closely related species, with 10 loci (out of 18 tested loci, i.e. 55.6%) in M. marginata and 6 loci (33.3%) in M. subnitida. Cross-species amplifi cation tests showed lower success rate for Frieseomelitta varia (22.2%) and Nannotrigona testaceicornis (11.1%), presumably due to the larger evolutionary distance between the target and the source species (Barbara et al., 2007).
The genotyping of the entire dataset revealed 70 alleles, ranging from one, for locus Mfsc34 to 10, for locus Mfsc17, with an average of 3.9 ± 2.7 alleles per locus (Table 2) Table 2, the level of polymorphism of each locus was also evaluated by the allelic richness (A R ) and the polymorphic information content (PIC). The values of allelic richness varied from 2 to 9.1 (average of 3.9 ± 2.4), while PIC values ranged between 0.108 and 0.714. Mean PIC (0.372 ± 0.198) characterizes all microsatellite loci as reasonably informative markers (Botstein et al., 1980). Overall mean observed and expected heterozygosity was estimated to be 0.536 and 0.453, respectively. These estimates were higher when compared to most heterozygosities found for Melipona species, exception made for M. subnitida (Souza et al., 2015). It is noteworthy that low levels of heterozygosity are known to occur in social Hymenoptera compared to other insects (Graur, 1985). Nine microsatellite loci exhibited signifi cant probabilities (P < 0.05) of departure from Hardy-Weinberg equilibrium, likely due to the mixing of individuals from populations of different allelic frequencies (Templeton, 2006). However, when populations were considered separately, the number of loci that deviated signifi cantly from HWE were fi ve in each population (Table  2), mostly caused by excess of heterozygotes. F IS estimates for most loci were negative, also indicating a trend towards an excess of observed heterozygosity. In a population this feature could be a signature of population bottleneck events (Cornuet & Luikart, 1996), which might be associated to destruction of native vegetation and consequent landscape transformation by human activities (Winfree et al., 2009;Potts et al., 2010).

Assessment of the genetic diversity in Melipona fasciculata
Genetic diversity between M. fasciculata populations, as measured by the mean number of alleles per microsatellite locus, mean allelic richness, heterozygosity and PIC, was characterized by a slightly higher level of genetic variability from samples collected in Maranhão when compared to those sampled in Piauí ( Table 2).
The high F ST (0.194) and R ST (0.230) estimates found in M. fasciculata suggest the existence of genetic differentiation. However, additional genetic surveys should be carried out to confi rm this observation. Similarly high F ST estimates were reported in wild populations of M. rufi ventris (Tavares et al., 2007) and M. beecheii (Quezada-Euan et al., 2007), with values of 0.250 and 0.280, respectively. D est , which is a measure based on the proportion of alleles that are unique to a subpopulation, provided further evidence of population differentiation (D est = 0.162). Low rates of dispersion and short fl ight distance, less than 2000 m, might have contributed to the levels of population differentiation (Silveira et al., 2002;Araújo et al., 2004;Duarte et al., 2014).
The scatter-plot of the principal coordinates analysis (PCoA) showed a clear separation of the species in two distinct clusters of stingless bees ultimately defi ned by the origin of each individual population, therefore, confi rming a signifi cant molecular genetic difference between the two populations (Fig. 1A). The analysis of microsatellite variation using the admixture model of STRUC-TURE, at the fi rst level of sub-population separation (K = 2), have also revealed two distinct clusters (Fig. 1B). These clusters rep-  resent each sampling population separately, except for very few individuals mostly located in Piauí that appears to be admixed. All these estimates indicate that there is little contemporary gene fl ow between these two M. fasciculata populations. This suggests questions about whether the destruction of native semiarid vegetation is increasing genetic drift by reducing genetic connectivity among M. fasciculata populations currently restricted to the remaining fragments of native forests.
Overall, analyses provide provisional evidence of signifi cant population differentiation between Maranhão and Piauí, in M. fasciculata. However, the data generated by this study should be further investigated using the same microsatellites markers, but larger sample size and more widespread sampling throughout the distribution of the species. Given all these considerations, the 18 isolated microsatellite loci have demonstrated strong potential for population-level genetic studies and can be used effectively as molecular tools to aid in the conservation of the species. Crossspecies amplifi cation further indicates that most of these loci can be useful across a wider range of species. Moreover, the results obtained from this snapshot assessment of the genetic diversity in M. fasciculata support their use for conducting population genetics and landscape genetics studies.

ACKNOWLEDGEMENTS. Embrapa (The Brazilian Agricultural
Research Corporation) funded this study through the project grant no. MP 06.14.01.001.00.00. This study was also fi nanced in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES). Sequencing was conducted on a MiSeq DNA sequencer purchased with a bequest from E.A. Nielsen to the Marine Gene Probe Laboratory. The authors are grateful to J.M. Vieira-Neto (in memoriam) for help in fi eld collection. We also thank two anonymous reviewers for comments that considerably improved the manuscript.