Evidence for Gradients of Melipona rufiventris (Hymenoptera: Apidae) Genetic Diversity within the Brazilian Semiarid

INTRODUCTION

Stingless bees belong to the family Apidae in the order Hymenoptera and are members of the tribe Meliponini, which includes all genera of extant stingless bees [1, 2]. Meliponines, in turn, are pollinators of the local native flora and cultivated crops [3-7], and therefore, of great ecological and economic significance, and adapted to survive in a variety of habitats, such as different natural forest environments, savanna, wetlands, protected areas, farmlands, wooden houses, among others [8]. The intensification of agriculture aiming solely at increasing food/forage production could lead to habitat loss and fragmentation, which have been recognized as crucial factors for the decline of stingless bee populations, mainly due to the inappropriate use of fertilizers and pesticides [9-12].

The species Melipona rufiventris Lepeletier, 1836 is part of the group of stingless bees native to Brazil, commonly known as uruçu-amarela, whose populations have suffered high losses in recent years by the increasing destruction of native semi-arid vegetation and the predatory collection of honey. It is currently included on the red list of threatened Brazilian fauna as an Endangered Apidae [13].

In addition to these constraints imposed by anthropogenic disturbances of the natural habitat, drastically reducing many populations, studies have shown that the virgin queens of M. rufiventris mate with only one male (monandry) during a nuptial flight [14]. So, all females will share the same father, inheriting the male's chromosomes. Genetic variability among their offspring could be, therefore, low. Although single mating is an ancestral condition in this group, low colony-level genetic diversity can not be taken for granted, since alternative mechanisms to increase intra-colonial genetic diversity can be developed to surpass natural biological constraints [15, 16].

Knowledge of the genetic variability and diversity in M. rufiventris and its populations, and how they are structured, is essential to a correct interpretation of population dynamics. Information on the population structure of the species can also be used to predict the potential impact of disturbances on the habitat and to formulate appropriate management strategies for species conservation [17]. Molecular markers can be used to measure genetic diversity and variability within and among populations of a single species from different geographical origins [18]. Hence, important questions of conservation relevance can be addressed to this threatened stingless bee species.

This study provides a snapshot of the genetic diversity and structure of Melipona rufiventris populations in the Brazilian semiarid using microsatellite markers, a starting point to understand the spatial and temporal genetic changes of M. rufiventris populations in Brazil.

MATERIALS AND METHODS

Bee materials and genomic DNA isolation

Worker bees were randomly collected from natural colonies distributed in 3 locations as follows: 25 nests in Campo Maior (CAM; coordinates: 4°49'19"S, 42°09'52"W), 7 nests in Castelo do Piauí (CAP; coordinates: 5°23'15"S, 41°31'17"W), and 6 nests in Guadalupe (GUA; coordinates: 6°47'30"S, 43°34'14"W), all in the states of Piauí-Brazil. All samples were taken to the laboratory and stored at -20°C until further use. Genomic DNA was extracted from each adult worker's thorax using the standard HotSHOT protocol [19]. Alkaline lysis buffer was heated to 95°C for 60 min; then, samples were cooled to 4°C and pH adjusted to 5 with 40 mM Tris-HCl. DNA was quantified using a NanoDrop 2000 Spectrophotometer (Thermo Scientific) and quality was checked using electrophoresis on 0.8% agarose gels.

Microsatellite markers testing and genotyping

Initially, cross transferability of heterologous microsatellite primers was tested in the M. rufiventris genome to find polymorphic loci (Table 1). PCR amplification was carried out using 19 microsatellite primer pairs developed from M. subnitida [20] and 18 pairs from M. fasciculata [8].

PCR reactions were carried out in a total volume of 10 μL containing approximately 10-50 ng of DNA, 1× PCR buffer (40 mM Tris-HCl; 100 mM KCl), 0.2-0.25 μM of primers, 2.5-3.0 mM MgCl₂, 0.5-1.0 μM of each dNTP, 0.25-0.75 U of Invitrogen Taq DNA polymerase. All amplifications were carried out using a VERITI^TM Gradient Thermal cycler (Life Technologies). PCR conditions (PCR₁) were as follows: 94°C for 5 min, 30 × (94°C for 40 s, Ta (50-60)°C for 30 s, 72°C for 40 s) and 72°C for 7 min. An alternate PCR profile (PCR₂) was used to optimize primers Mfsc11, based on the following conditions: 94 °C for 1 min, 40 × (94 °C for 30 s, Ta°C for 30 s, 72 °C for 30 s) and finally 72 °C for 3 min. SSR markers were screened by silver nitrate detection on denatured 6% polyacrylamide gels. A 10-bp ladder (Life Technologies) was used as a size marker. Amplifications of microsatellite loci were considered successful when the gel presented one or two clear and consistent bands, with product sizes close to those of the original species.

Table 1. Markers, primer sequence, experimental parameters used for the amplification of microsatellite loci in Melipona rufiventris, and their GenBank accession number.

Ta, Annealing temperature; PCR profiles: (PCR₁ = [94°C-5 min; 30 cycles × (94°C-40 seg; Ta-30 seg; 72°C-40 seg); 72°C-7 min], PCR₂ = [94°C-1 min; 40 cycles × (94°C-30 seg; Ta-30 seg; 72°C-30 seg); 72°C-3 min].

Data analysis

The genotyped data were analyzed using Micro-Checker 2.2.3 [21] to test for the presence of null alleles or possible scoring inconsistencies. The number of alleles (A), observed and expected heterozygosities (H_O and H_E), and the polymorphic information content (PIC) were determined using cervus 3.0.3 [22]. Allelic richness (A_R) was calculated by fstat version 2.9.3.2 [23]. Tests for Hardy–Weinberg Equilibrium (HWE) and linkage disequilibrium were conducted using the web-based software genepop [24]. Bonferroni-corrected P-values were used to assess the significance (P<0.05). A Bayesian grouping admixture model was used to infer possible population structuring using the software structure v2.3.3 [25]. The program was set up for 1,000,000 Markov chain Monte Carlo repetitions after an initial burn-in of 100,000 steps. The estimate of the best K was calculated based on 10 replications for each K (from 1 to 4) as described by [26] using structure harvester v.0.6.92 [27]. The program clumpp v.1.1.2 [28] was used to align the five repetitions of the best K. The program distruct v.1.1 [29] was used to graphically display the results produced by clumpp. Population structure was also analyzed using principal coordinate analysis (PCoA) as implemented in GenAlEx v.6.5 [30].

RESULTS AND DISCUSSION

Several genetic diversity studies have been conducted on bee populations using transferred microsatellite markers [31-35].

After PCR testing of 37 potential heterologous microsatellite primer pairs, only nine markers (24.3%) revealed polymorphism in M. rufiventris. A monomorphic primer (Mfsc13) amplified clear and consistent bands, however, only one band was detected within the studied populations. All the other primer pairs (73%) produced unsatisfactory results with strong non-specific banding patterns or showed no amplification, even though M. subnitida and M. fasciculata are congeners of M. rufiventris. According to Silva et al. [8] amplification success declines according to the genetic distance between the taxa, that is, the phylogenetic proximity is the main success factor in the transferability of heterologous primers. Other factors such as the size and complexity of the genome and the location of the microsatellite, whether in the coding region or not, may also influence the transferability of microsatellite primers.

Four loci, Mfsc11, Msub31, Msub38, and Msub51, significantly (P < 0.05) departed from HWE when the data set was examined as a whole, generally due to the presence of null alleles or missing data. These loci showed heterozygote excess and possible causes are the artifactual scoring of nonspecific and stuttering bands [21], as well as the small effective population size [36, 37]. Allelic richness (A_R) ranged from 2 to 6.6. with a mean of 3.2, while PIC values ​​vary between 0.12 and 0.67, with an average of 0.37, therefore, moderately informative. PIC values ​​above 0.5 are considered very informative; moderately informative between 0.25 and 0.5 and less informative lower than 0.25 [38]. The observed heterozygosity (Ho) varied from 0.00 to 0.85 with a mean of 0.47, whereas the expected heterozygosity (He) presented values ​​between 0.14 and 0.72, with a mean of 0.43. The Msub31 locus was the most polymorphic, and Msub51 was the only locus that presented some evidence of null alleles, but at a lower frequency than 0.200 (Table 2). Null allele frequencies below 0.200 are acceptable in most microsatellite data sets [39].

The mean He was the same as that found by Lopes et al. [40] for M. rufiventris (He = 0.43) and higher than those found in studies with other meliponine species such as 0.38 for M. mondury and 0.35 for M. mandacaia [41], and 0.105 for M. mondury and 0.189 for M. quadrifasciata [31]. The genetic diversity found in M. rufiventris although superior to most of the aforementioned studies is considered low, and this reduced value is associated with genetic, biological, and environmental aspects, as is the case of anthropogenic interferences. These interferences contribute to the reduction of genetic variability through habitat fragmentation and predatory actions of honey collection, which reduce the number of colonies in a given locality [40-44].

When the data set was divided into 3 collection sites (Campo Maior, Castelo do Piauí, and Guadalupe) it was observed that Campo Maior presented the highest average in relation to the number of alleles/population (3.0), ranging from 1 to 7, when compared with Castelo do Piauí and Guadalupe. Departure from Hardy-Weinberg equilibrium was observed at Mfsc11, Msub18, Msub31, and Msub38 loci in the same collection site, after Bonferroni correction. No loci in samples from Castelo do Piauí and Guadalupe deviated significantly from HWE, P > 0.05 (Table 2).

Table 2. Variability of 9 microsatellite loci and estimates of genetic diversity in Melipona rufiventris within the Brazilian semiarid.

A: Number of alleles in the population; Ho: Observed heterozygosity; He: Expected heterozygosity; PIC: Polymorphic information content; pHWE: Hardy-Weinberg equilibrium probability; *: 5% significance (Bonferroni correction < 0.005).

When estimating the values ​​of the F statistic, based on these molecular markers, the Fis was -0.177, which indicates a low level of inbreeding, whereas, the overall Fst and Rst were 0.151 and 0.288, respectively. The value found for Fst was of similar magnitude compared to that estimated in a previous report for Melipona asilvai populations using microsatellite markers (Fst = 0.166) [45]. Although the FST value (0.151) obtained in this study, based on heterologous markers, is not very high, it indicates some level of structure in the populations, which becomes particularly relevant for a threatened species such as M. rufiventris. Also according to Nei [46], Fst values ​< 0.05 are considered low, between 0.05 and 0.15 moderate, and Fst values ​> 0.15 are high, thus indicating a high population structure within the study area.

The principal coordinate analysis (PCoA) showed a separation of the species into two main groups with some degree of intersection, thus confirming a significant genetic differentiation among the three populations (Figure 1a). In the Bayesian analysis, the optimal K-value was determined to be 2, where the analysis revealed two distinct groupings (Figures 1b-1c).