Genetic Diversity in Nannotrigona testaceicornis (Hymenoptera: Apidae) Aggregations in Southeastern Brazil.

The Meliponini, also known as stingless bees, are distributed in tropical and subtropical areas of the world and plays an essential role in pollinating many wild plants and crops These bees can build nests in cavities of trees or walls, underground or in associations with ants or termites; interestingly, these nests are sometimes found in aggregations. In order to assess the genetic diversity and structure in aggregates of Nannotrigona testaceicornis (Lepeletier), samples of this species were collected from six aggregations and genetically analyzed for eight specific microsatellite loci. We observed in this analysis that the mean genetic diversity value among aggregations was 0.354, and the mean expected and observed heterozygosity values was 0.414 and 0.283, respectively. The statistically significant Fis value indicated an observed heterozygosity lower than the expected heterozygosity in all loci studied resulting in high homozygosis level in these populations. In addition, the low number of private alleles observed reinforces the absence of structuring that is seen in the aggregates. These results can provide relevant information about genetic diversity in aggregations of N. testaceicornis and contribute to the management and conservation of these bees' species that are critical for the pollination process.

Meliponini are popularly called indigenous stingless bees (Michener 2007). Today, over 600 species in 56 named genera live in tropical and subtropical areas of the world (Cortopassi-Laurino et al. 2006). They are common visitors to flowering plants in the tropics and they are known to visit the flowers of ∼90 crop species (Heard 1999, Slaa et al. 2006). Depending on the biome, the need for pollination by stingless bees can range from 30 to 90% of savanna species, including in some rainforest patches and Amazon (Kerr et al 2001). Thus, these bees play an important ecological and economic role that makes them a very important group for pollination (Biesmeijer and Sla 2006, Klein et al. 2007).

A number of species of stingless bees can form their nests in a highly aggregated spatial distribution (Starr and Sakagami 1987; Eltz et al. 2003). The queens with philopatric behavior, in which the daughter colony is supplied by the mother colony for a long period and thus limiting the distance that the swarm could move in relation to the mother nest, could have critical influence in the formation of these aggregations (Wille 1983, Inoue et al. 1984). This behavior presents a strong impact on population structuring, and may have had an important evolutionary role for the great species diversity reported for Meliponini (Brito and Arias 2010).

Few studies have been conducted on genetic diversity in stingless bees and even fewer in bees that form nest aggregations. Therefore, genetic data in this area is essential for a better understanding of the dynamic of these organisms and especially for the development of conservation strategies (Moritz 2002, Cortopassi-Laurino et al. 2006).

In studies carried out with native stingless bees, molecular techniques have been increasingly used (Arias et al. 2006, Quezada-Euán et al. 2007, Francisco et al. 2008). In the last decade, different molecular markers have been characterized, most of which are derived from different sequences of mitochondrial (mtDNA) and microsatellite loci. Microsatellites have been extensively used in parentage testing, intraspecific variation, hybridization between species, population dynamics, genetic mapping, and phylogeographic studies (Moritz and Hillis 1996, Chakraborty and Kimmel 1999).

The aim of this study was to evaluate the genetic diversity and structure of Nannotrigona testaceicornis aggregations by specific microsatellite markers in six nests aggregations localized in Southeastern Brazil.

Material and Methods

Collection of Bees

About 302 individuals were collected at the entrance of 32 natural nests distributed among six N. testaceicornis aggregations in Southeast Brazil, two at the University of São Paulo Campus of Ribeirão Preto (SP) (SUL and OFI), one in Campinas (SP) (CAM), one in Bonfim Paulista (SP) (BON), one in Uberlândia (MG) (UBE) and one in Caratinga (MG) (CAR) (Fig. 1). The geographic coordinates were also collected from the six aggregations. The distance between nests of each aggregate did not exceed five meters. Bees were placed in 70% ethanol and stored at −20 °C until DNA isolation.

Fig. 1.

Location of the aggregations in the states of São Paulo and Minas Gerais, Brazil. 1-SUL; 2-OFI; 3-CAM; 4-UBE; 5-BON; 6-CAR.

Genetic Analysis

For DNA isolation, Wizard Genomic DNA Purification Kit (Promega, Madison, WI) was used according to manufacturer’s protocol, and then stored at −20 °C. Eight pairs of microsatellite markers (primers) of N. testaceicornis (Ntes02, Ntes34, Ntes37, Ntes41, Ntes50, Ntes59, Ntes70, and Ntes73; Oliveira et al. 2009) (Supp Table 1 [online only]) were evaluated in ten individuals of each nest. DNA amplification was performed by polymerase chain reaction (PCR) according to Oliveira et al. (2009). The electrophoresis was performed using 12% denaturing polyacrylamide gel (acrylamide:bis-acrylamide = 29:1) containing 7 M urea in 1× TBE buffer. The DNA bands were visualized using a DNA silver staining (Sanguinetti et al. 1994). Scoring of the PCR bands from the silver-stained gels determined the sizes of the alleles. The alleles were arbitrarily named according to the size of the PCR bands and then compared with Oliveira et al. (2009).

Statistical Analyses

Statistical analyses were performed to estimate the allelic and genotypic frequencies, adherence to Hardy-Weinberg (H-W) equilibrium and U test using GENEPOP software, v3.0 (Raymond and Rousset 1995). The Genetic Data Analysis Package GDA software (Lewis and Zaykin 1997) was used to estimate heterozygosity and the F statistics. The ADZE software (Szpiech et al. 2008) was used to calculate the allelic richness and the number of private alleles by aggregation and in all possible combination of aggregations. Because of the small sample size of the CAM aggregation, were also carried out analyzes of private alleles in the absence of that aggregation.

The population structure was assessed by analysis of molecular variance (AMOVA), using the ARLEQUIN software, v3.1 (Excoffier et al. 2005) and the Bayesian clustering approach implemented in STRUCTURE v.2.3.3 (Pritchard et al. 2000). The program was set to 1.000.000 interactions and burnin of 500.000 under the admixture model. Models with K from 1 to 6 were repeated 10 times. The best K was estimated as described by Evanno et al (2005) using STRUCTURE HARVESTER software v.0.6.94 (Earl and vonHoldt 2012).

GENELAND (Guillot et al. 2012), available in the R statistical package, was used to perform a spatial genetic analysis by integrating geographic and genetic information, assuming the genotype matrix of each individual, the range of possible clusters and a maximum of 10 populations. Ten independent runs with 106 interactions Markov Chain Monte Carlo (MCMC) were performed. Each one with the correlated allele frequencies model and geographical coordinates were used to analyze the datasets.

Correlation of geographic and genetic distances was determined by Mantel regression of Fst values on log10-transformed geographic distances among populations using isolation by distance (IBD) 1.52 (Bohonak 2002).

Results

Among the six aggregations, the total number of alleles was 38, ranging from one to 11 per locus: the locus Ntes02 was the only one to show only one allele and the Ntes 50 was the locus with the largest number of alleles (n = 11). The CAR aggregation showed 30 of the 38 observed alleles and the number in other aggregations ranged from 15 to 22 alleles (Supp Table 2 [online only]). The allelic richness ranged from 1.87 to 3.75, with average richness among the six aggregations of 2.76 alleles. Rarefaction analysis showed that the aggregations presented a very small number of private alleles, with numbers ranging from 0.03 (CAM) to 0.9 (CAR) (Fig. 2). Despite this small number of private alleles (<1), we analyzed their distribution in all combinations of aggregations and did not observed any correlation between private alleles and geographic distance (Supp Fig. 1 [online only]).

Fig. 2.

Number of private alleles per aggregation.

The heterozygosity observed in all aggregations was lower than expected, indicating low gene diversity. The average diversity between aggregates was 0.354. A summary of the genetic variation between aggregates is shown in Table 1.

Table 1.

Genetic variation among aggregations

	SUL	OFI	CAM	UBE	BON	CAR
N	36	88	18	54	56	50
Nn	4	9	2	6	6	5
Na	30.8	78.4	15.0	51.9	52.9	49.9
T	17	22	15	19	22	30
A	2.24	3.37	1.87	2.38	2.94	3.75
PLP (95%/99%)	50.0	62.5	62.5	50.0	75.0	87.5
Hea	0.295	0.366	0.326	0.307	0.389	0.435
Ho	0.144	0.218	0.234	0.272	0.349	0.434
DH	0.1222	0.0147*	0.2958	0.1692	0.3679	0.3965
Fis	0.525	0.407	0.294	0.115	0.102	0.002

N, number of individuals per aggregation; Nn, number of nests per aggregation; Na, average number of individuals per locus; T, total alleles per aggregation; A, average allelic richness; PLP, polymorphic loci; He, average expected heterozygosity; Ho, average observed heterozygosity; DH, average deficit of heterozygotes (*P < 0.05); Fis, average inbreeding coefficient (all values were significant at CI = 95%).

aThe average He among populations was 0.353 and Ho was 0.275.

Considering the values of polymorphic loci (PLP), the highest percentage was observed for the CAR aggregation and the lowest for the UBE and SUL aggregations. No aggregation showed difference in PLP, regardless of the criteria value used (95 or 99%).

The H-W equilibrium test showed a highly significant P value (P < 0.05) for all aggregates, showing that they are not in H-W equilibrium (data not shown). To verify if the H-W equilibrium deviation in N. testaceicornis aggregations was related to the excess or deficit of heterozygotes, the “U” test was performed (Rousset and Raymond 1995), and showed that all aggregates presented deficit of heterozygotes (DH) for at least one marker (Supp Table 3 [online only]). However, when for the average DH, only OFI aggregation presented statistically significant results.

The inbreeding coefficient (Fis) was calculated per locus and aggregation. All values were statistically significant for confidence interval (CI) of 95% (Table 1).

The population differentiation parameter Fst showed significant values between pairs of aggregates, demonstrating that they are genetically different from each other at CI = 95% (Table 2).

Table 2.

Estimates of Fst between pairs of aggregates

Aggregations	SUL	OFI	CAM	UBE	BON
OFI	0.1389
CAM	0.2050	0.1347
UBE	0.2390	0.1159	0.1972
BON	0.0734	0.0772	0.1308	0.1204
CAR	0.2446	0.1459	0.2400	0.2191	0.1646

The pairwise combinations were all significant at CI = 95%.

The genetic variation was assessed by AMOVA. First, all comparisons among the six aggregates were tested. The percentage of variation for this analysis was 14%, with P < 0.001 (SD = 0), indicating a low but significant value. Therefore, each aggregation was evaluated individually in order to verify whether there was genetic variation within the aggregations. This analysis however, did not show statistically significant P value.

The aggregations structure analysis using the STRUCTURE, identified five clusters (K = 5) according to the STRUCTURE HARVESTER. However, there was no structuring among aggregations. Mixture was present in all populations (Fig. 3).

Fig. 3.

Structure analysis among the six aggregations. STRUCTURE showing no structure between the aggregations (K = 5).

The GENELAND model also identified five spatially coherent clusters with the following groups and probability respectively: SUL (0.45), OFI (0.40), CAM (0.45), UBE (0.45), BON (0.48), and CAR (0.60). However, this structure was undetected under a model of uncorrelated frequencies, possibly due to low levels of observed genetic differentiation between the aggregations.

The correlation between geographical distances and genetic differentiation (IBD) was not statistically significant between aggregations (r = 0.5043; P = 0.0537).

Discussion

In this study, we analyzed the genetic diversity in aggregations of N. testaceicornis. The allelic diversity is one of the parameters used to study the genetic diversity of populations. Our analyses showed values ranging from 1.87 to 3.75. Although the CAM aggregation showed low allelic richness (1.87), this value is justified due to the limited number of individuals (18).

The exact test showed a highly significant deviation for the H-W equilibrium with a large imbalance between homozygotes and heterozygotes, and in all aggregations, there was a deficit of heterozygotes for at least one marker.

The low to moderate Fst values indicated a small structure in the aggregations also identified by GENELAND (two genetic groups are observed: group I (OFI, SUL, and CAM) and group II (UBE and CAR) but not identified by the STRUCTURE. The results observed in GENELAND were probably because this software successfully clusters individuals into groups, even with low or moderate Fst values. We believe that this variation, in spite of statistical significance, is too small to structure a population. Similar results were observed by Viana et al (2013) in a study about the genetic variability in M. scutellaris.

Private allelic richness, is a convenient measure of how distinct a population is from other populations (Kalinowski 2004). The low number of private alleles identified by aggregation or in combination was against the results obtained with the Bayesian clustering methods. The mixture of aggregates observed in the structure analysis could be related to the lack of correlation between genetic and geographical distance (IBD) that we have identified between these aggregates, thereby enabling the existence of a gene flow, e.g., by drone migration.

Many Hymenoptera species are found in aggregations (groups of nearby but unrelated nests) and the spatial distribution of these nests may be influenced by biotic and abiotic factors (Potts and Willmer 1998). Oldroyd et al. (1995) raised some hypotheses to explain this phenomenon in Apis mellifera: 1) a result of the preference of swarms to travel short distances from their original nest, 2) a result of an environmental imbalance, 3) increase in nest defense, 4) improved mating efficiency. Unlike Apis bees, who leave their colonies at once and never return, swarming among stingless bees is slow, with high dependence on the daughter colony in relation to the mother colony. For many days, stingless bees supply the new nest with food and building materials from the mother colony (Nogueira-Neto 1997). This relationship could limit the displacement of the reproductive swarm in relation to the parental nest (Cameron et al. 2004), increasing the inbreeding index. On the other hand, studies by Cameron et al. (2004) with nest aggregation and male Trigona collina congregations suggest that this species minimizes the chance of brother–sister mating and maximizes the probability of mating of queens with unrelated males. These authors found that each aggregation nest showed no tendency to be more related than isolated colonies, probably because males move through the aggregates to reach a farther congregation. As filial nests in stingless bees swarming can only be founded in the proximity of the mother colony, panmixis is effected mainly by drone migration (Engels and Engels 1984).

The panmixia by drone migration could explain why the mixture of alleles between aggregations. However, we still observed a low number of alleles and a high rate of homozygosity.

The low levels of genetic diversity detected could suggests that stochastic population processes such as genetic drift and small and fluctuation population sizes that results in population bottlenecks, have had a strong effect on the maintenance of genetic variation (Hinten et al. 2003).

Thus, we suggest two situations for these aggregations: 1) one in which they were separated by barriers that over the time become transposable; or 2) these bees formed a continuous population and over the time started to isolate themselves, still being in this process.

If geographically separated populations become progressively more genetically isolated over time, gene flow is expected to yield to drift as the predominant mechanism shaping population structure (Clegg and Phillimore 2010).

When the habitat is fragmented, dispersion and colonization potential is often reduced. Many species of birds, mammals and insects inside forests are unable to cross short distances of open areas (Lovejoy et al. 1986, Bierregaard et al. 1992). Araújo et al. (2004) found a correlation between body size and flight radius in stingless bees, that is, the bigger the bee, the greater the distance traveled during the flight, and concluded that body size could act as a limiting factor in the flight capacity. According to the authors, the maximum flight distance constrained by body size can have a direct influence on the dispersal capacity of the population. For N. testaceicornis bees, Araújo et al. (2004) observed flight capacity of 951 m. Therefore, bees of this species could not overcome open areas greater than this distance to establish new nests.

Despite the apparent high adaptability of several species of stingless bees, including N. testaceicornis to urban environments (Souza et al. 2002), our data revealed that the aggregations collected at six different locations in two Brazilian states showed low genetic diversity. Regardless of the possible reasons that may be causing this phenomenon, such as deforestation and fragmentation of forests, (and others not cited in this study such as the introduction of exotic species, the predatory collection of honey, and the intensive use of herbicides and pesticides) (Freitas et al. 2009), it is essential the maintenance of the genetic variability for conservation of these populations due to their essential role in pollination.

Supplementary data

Supplementary data are available at Journal of Insect Science online.