Analysis of the population structure of Tetragonisca (hymenoptera, meliponini) by microsatellite markers and network interactions

Tetragonisca angustula is the stingless bee more managed by beekeepers, whose honey has many medicinal properties and is distributed throughout Latin America. They have variation in color of mesepisternum and can be black, yellow or mixed. They might be classified into two species (T. angustula and T. fiebrigi) or two subspecies (T. angustula angustula and T. angustula fiebrigi). This study aimed study was to analyze the genetic variability of populations T. angustula sampled using the technique simple repeated sequences and to analyze the interaction of molecular markers used in other studies to verify the possible existence of a marker that distinguishes genetically. A total of 60 individuals were sampled in the cities of Maringá, Cianorte and Terra Boa, Paraná state, Southern Brazil. The results showed that all loci are polymorphic. The value of FST=0.1173 shows that populations are moderately differentiated and analysis of molecular variance indicates that 76% of the change occurs within the analyzed populations. The largest value of delta K obtained by Bayesian inference estimated the actual number of people equal to three. The analysis of interaction networks has shown that there are more interactions with isozymes, but none of them allowed separate T. angustula in two species.


Introduction
Bees of the Meliponini tribe (Hymenoptera, Apidae) are considered highly eusocial (Hrncir et al., 2016), being known as native bees, indigenous or stingless bees, exhibit distribution in the tropical and subtropical regions of the world (Camargo & Pedro, 2008). These pollinators showing flower visits (Giannini et al., 2015), are deeply related to human well-being, performing different activities that promote the maintenance of ecosystem health and function, breeding of wild plants, crop production, and food security, as well as producing honey and other beekeeping products (Potts et al., 2016), being ecologically, economically and culturally important (Ayala et al., 2013).
Several studies involving biochemical and molecular markers were conducted with T. angustula aiming to analyze the genetic structure of the population of this bee. However, there is no consensus on their classification, with some authors classifying them into two subspecies (T. angustula and T. angustula fiebrigi) (Schwarz, 1938), and others in two species (T. angustula and T. fiebrigi) (Camargo & Pedro, 2008). Both are based on morphological differences, separating them by the pigmentation on the thorax called mesepisternum (Castanheira & Contel, 2005).
Other studies support the specific name described by Camargo & Pedro (2008). Among them, Stuchi et al. (2012), using esterase isoenzymes found in T. fiebrigi workers esterases one, two and four and in T. angustula esterases three and four.
Also, using analyzes involving the use of PCR-RFLP (Koling & Moretto, 2010), and analyzing karyotypic (Barth et al., 2011), divergence between T. angustula and T. fiebrigi, with the presence of B chromosomes occurring in T. fiebrigi. Moreover, analysis of the sequencing mitochondrial DNA of microsatellite regions by Francisco et al. (2014), reported the occurrence of hybridization between T. angustula and T. fiebrigi, and that the divergence between these species is recent and with asymmetric introgression, suggesting the subspecies classification. Santos et al. (2014) analyzed the genetic variability in T. angustula samples, T. fiebrigi and hybrids of three Brazilian states (Paraná, Rondônia, and São Paulo), utilizing PCR-RFLP. They found that T. angustula and T. fiebrigi are not completely separated, there is a hybridization between them and this was observed between T. angustula females and T. fiebrigi males. In studies using the RAPD marker it was not possible to differentiate at the specific level, only populational (Baitala et al., 2006).
Obtaining genetic markers that may distinguish these bees is an important tool, as these will contribute to studies on the distribution of these stingless bees in neotropical areas, and serve as a basis for planning conservation strategies both in the ecosystem and in meliponaries.
Interaction networks is a methodology that uses graphical means to facilitate understanding of the dynamics of mutualistic relationships because they allow a representation of its complexity and an assessment of the entire structure of interactions (Bascompte & Jordano, 2007). Its widely used in management studies, environmental impact, ecological (plantpollinator interaction, host-parasite) in silico and in social networks and metabolic pathways within cells (Albert & Barabasi, 2002). This study aimed to analyze the genetic variability of the population of T. angustula sampled using the technique SSR and verify the interaction of molecular markers used in other works with the morphological variants of T. angustula to check the possible existence of a marker that differentiates them genetically.

Sampling location
Were used adult workers of T. angustula and T. fiebrigi collected in four Paraná state locations: Maringá (six nests), Cianorte (one nest), Terra Boa (four nests) and Malú (one nest), district of Terra Boa (Figure 1). The bees were collected randomly in pet bottles positioned in the straw nest. They were sacrificed by freezing and stored in vials with 70% alcohol. Souce: Authors.
The morphologic analysis was performed according to the color of the mesepisternum proposed by Camargo & Pedro (2008) (Figure 2). Of the twelve nests, six were natural, from tree holes or exposed holes and pipes in walls. The other six nests Research, Society and Development, v. 10, n. 4, e4711424811, 2022 (CC BY 4.0) | ISSN 2525-3409 | DOI: http://dx.doi.org/10.33448/rsd-v10i4.24811 4 are from rational boxes that farmers reported that they had not been divided, only used baits for the bees to install themselves in them.

DNA extraction and Genotyping
The extraction of nuclear DNA methodology was performed as described by Yu et al. (1993), but using the chest of five workers of each nest instead of the head. Quantitation of DNA was assessed by spectrophotometer Picodrop (PICOPET01 model) following the manufacturer's recommendations. Polymerase Chain Reaction (PCR) amplification was performed on all individuals genotyped at four microsatellite loci (Table 1).  Amplifications via PCR were performed in a thermocycler (Biosystems, Veriti 384) as the original protocol described by Brito et al. (2009). The amplification process was started with an initial denaturation at 94 °C for 4 minutes followed by 35 cycles of 30 seconds at 94 °C, annealing for 30 seconds at a specific temperature for each primer, the extension for 30 seconds at 72 °C and 5 minutes at 72 °C for the final extension.
The products of the amplifications were subjected to electrophoresis on 4% agarose gel (prepared with 50% MS8 agarose and 50% regular agarose) at 60 V using the total volume of the reaction (20 µL) added to two µL load (bromophenol blue [0.25%] and glycerol [30%]). To determine the size of amplified fragments was used a molecular weight marker of 100pb DNA ladder (Invitrogen). After electrophoresis, the gels were stained in ethidium bromide solution (0.5/ mL) for 60 minutes, followed by the revelation of the bands by UV light and photo-documented in L-Pix HE Loccus biotechnology system.
The molecular analysis of data involved the interpretation of nuclear DNA fragments obtained by agarose gel electrophoresis. It was calculated measures of proportion polymorphic loci, the mean number of alleles observed heterozygosity (HO) and expected heterozygosity by the Hardy-Weinberg equilibrium (HE). To determine the significance of genetic differentiation between the population was performed the analysis of multi-loci, which were calculated the fixation indexes FIS and FST using POPGENE software version 1.32 (Yeh et al., 1999). The analysis of molecular variance (AMOVA) was obtained by GenAlex software version 6501 (Peakall & Smouse, 2012).
The estimate of the actual number of populations was performed according to Evanno et al. (2005), using delta statistical K = ln (P (X/K) between successive values of K (number of population) and where X = probability of observing the data. The estimation was made with K values ranging from two to eight and the value of K used was the one that showed higher delta K. This estimate was performed using the Structure 2.3.4 software (Pritchard et al., 2000) and structure Harvester (http://taylor0.biology.ucla edu/struct_harvest/). The burn-in was set at 3000 repetitions and the Monte Carlo method Markov chains (MCMC) in 30 thousand whit set number of iterations was ten.

Analysis of interaction networks
The analyses were performed using data from six studies. These works include molecular analyzes with a different population of Tetragonisca sp. with microsatellite (Bronzato, 2011, this study), isoenzymes (Ruiz, 2006, Silva, 2007, Stuchi et al., 2012 and Polymerase chain reaction-restriction fragment length polymorphism (PCR-RFLP) (Santos et al., 2015).
Data were transformed into a binary matrix to build the network qualitative interactions, that when an individual has given allele, it has received the number one (presence of the allele), and when the allele was not present in an individual, it has received the number zero (absence of the allele). The calculations and graphs of interaction networks were obtained from the bipartite package (Dormann et al., 2008) R software version 2.14.1 (R Development Core Team, 2011).

Results
This section may be also divided by subheadings. It should provide a concise and accurate description of the experimental results, their interpretation as well as the experimental conclusion that can be drawn.

Morphological identification
Bees collected were separated into two populations based on the geographical distance between the nests ( Table 2).
The first population is made by bees that were collected in the city of Maringá and the second by the bees collected in the cities of Cianorte and Terra Boa, which are distant 20 km apart. The city of Maringá is 74 and 64 km away from Cianorte and Terra Boa, respectively.

Genotypic and allelic distribution
All four loci analyzed showed polymorphism at the genotyping of 60 individuals, presenting 15 genotypes. The twopopulation showed a different number of genotypes with twelve genotypes in the population one and ten in population two.
Besides, the population one exhibited one heterozygous genotype more than the second one (three and two heterozygote genotypes, respectively). The number of alleles per locus detected varied from two (Tang11 and Tang77) to four (Tang12 and Tang78), with an average of three per locus (Table 3). At the Tang11 locus, the 180 bp allele was the most frequent in both populations, ranging from 43.33% to 95%. The 180 bp allele of locus Tang12 appears only in population one, at a frequency of 6.67%. This allele is 18 presents in homozygosis in two individuals from the same nest, classified as a hybrid. however, the 210 bp allele of this same locus is present only in population two, at a frequency of 10%.
It is also present in homozygosity, but in three individuals from a natural nest, classified as T. angustula. In this locus, the most frequent allele is 200 bp, ranging from 56.67% to 83.33%. At the Tang77 locus, the 220 bp allele is the most frequent in both populations, appearing with 70% in population one and with 73.33% in population two. At loco Tang78, the most frequent allele is 250 bp, varying from 46.67%, in population one, to 51.67%, in population two. however, the 310 bp allele of these 19 loci appears exclusively in population two, with a frequency equal to 1.67% in the heterozygous form of the Cianorte nest.
In both populations observed heterozygosity is lower than expected, with an excess of homozygotes that caused deviation from Hardy-Weinberg equilibrium (p <0.05). These results were confirmed with the FIS values (fixation index) positive and close to one (Table 4), demonstrating that there is an excess of homozygotes in the two populations. Research, Society andDevelopment, v. 10, n. 4, e4711424811, 2022 (CC BY 4.0) | ISSN 2525-3409 | DOI: http://dx.doi.org/10.33448/rsd-v10i4.24811 8 AMOVA Table 5 shows the results of the AMOVA. It is observed that most of the genetic variability is derived from its intrapopulation component (76%), which may suggest gene flow between upcoming groups of samples, and only 18% of the estimated variance stems from its interpopulation component. According to the best estimate of K, the largest structure found divided the analyzed population into three groups (K = 3) (Figure 3). One group is made of four nests, the nest of bees with hybrid staining from Malú district, and three from Maringá, one of T. angustula and two of T. fiebrigi (light gray area).
A second group is formed by a nest of T. fiebrigi from Maringá and the three nests of bees with hybrid staining, two from Maringá and one from Terra Boa (black area). The last group is formed by the nest of Cianorte city and three nests of T.

Interaction networks
From the data obtained from six works was generated a binary matrix of 894 rows, comprising individuals (555 classified as being T. angustula, 319 as T. fiebrigi and 20 with hybrid color) by 93 columns, which were the molecular markers' allele used in each study. Thus, it was possible to account 7874 interactions for this network ( Figure 4A). The alleles Tang11b, Tang12a, Tang12b, and Tang78c has little interaction with the group identified as hybrids in the study performed ( Figure 4B).

Discussion
The results demonstrate that the evaluated Tetragonisca populations have reduced genetic variability, suggesting that there is a homogeneity of the evaluated individuals. This can occur due to the short distance between the sampling sites, as well as the populations of this region being on the same evolutionary line (Santiago et al., 2016). Also, this genetic variability can occur due to the species being easily managed and adapted to the urban environment and the meliponicultores exchange colonies, divisions and capture via bait, causing a genetic exchange to occur, resulting in a homogenization of these populations (Francisco et al., 2017).
Microsatellite analyzes show similarity between the different species of T. angustula, hybrid and T. fiebrigi evaluated.
Despite being considered taxonomically different species, T. fiebrigi bees have a more restricted geographical distribution (Argentina, Bolivia, Paraguay, and in Brazil, in the states of Mato Grosso, Mato Grosso do Sul, Paraná, Rio Grande do Sul and São Paulo). T. angustula, found from southern Mexico to Argentina and occurring throughout Brazil (Camargo & Pedro, 2008). The molecular studies of mitochondrial DNA and analysis of microsatellites demonstrate that hybridization and introgression are taking place between T. angustula and T. Fiebrigi, being considered as subspecies (Francisco et al., 2014).
The analysis of the four loci in these populations revealed only 12 alleles fewer to the 53 detected by Brito et al. (2009) with the same four loci. The authors using the same markers obtained an average of 13.25 alleles per locus, ranging from ten to 16 alleles by analyzing 19 nests. It is believed that the smaller number of alleles found in this study is due by the geographical proximity of the nests which bees were collected when compared to the work of Brito et al. (2009), that used bees collected from nests in several regions of Brazil.
Through the FST value can be observed that the population is moderately differentiated since the estimated value was 0.1173 and according to Wright (1969) FST values between 0.05 and 0.15 indicate moderate differentiation. This moderate degree of differentiation of populations that are close may be the result of two processes: fertilization by a single male and the philopathy of the female. In most meliponines, queens are fertilized by only one male, a condition called monoandria (Peters et al., 1999).
This behavior generates less genetic variability within a colony, when compared to the polyandric behavior present in A. mellifera, for example. Low genetic variability is harmful to populations because it is related to low adaptive capacity and to the increase in their probability of extinction (Francisco et al., 2017). However, there is another negative consequence, the effects of which are practically immediate.
In bees, low genetic variability has been found to decrease immunocompetence in Bombus muscorum (Whitehorn et al., 2011) and increases the occurrence of infections in Apis mellifera and Bombus sp. (Baer & Schmid-Hempel, 2001, Tarpy, 2003, Cameron et al., 2011, Whitehorn et al., 2011. Therefore, this is yet another process by which populations with low genetic diversity can decrease their population size and even become extinct. Thus, it is possible that the number of related genotypes causes the present analysis to result in a deviation from the Hardy-Weinberg equilibrium. The small number of nests and loci analyzed may be another factor that contributed to the small variability detected, due to the low allelic polymorphism observed. As for the restricted dispersion of the females, this is due to their nesting habits, because when a new nest is about to form, a virgin queen is accompanied by a group of workers from the old nest. This new nest is dependent on the old nest until it becomes independent, and may take a few days or even months depending on the species (Nogueira-Neto, 1997).
In this way, the distance from the new nest is restricted to that of the workers' foraging. It is known that in T.
angustula the distance of foraging covers a radius of about 500m (Nogueira-Neto, 1997). Since the female's philopatry was detected in these bees, the absence of genetic structure can be explained by the dispersion of the males.
The few genetic studies of male Meliponini congregations have demonstrated the presence of males from distant areas, with more than 100 colonies present (Paxton, 2000, Cameron et al., 2004, Kraus et al., 2008, Mueller et al., 2012. These factors together (monandry, nesting distance and dependence on the old nest) result in a high degree of kinship between the nests, with no great genetic differentiation between populations. Thus, the division of these groups may not be related to distance or geographical barriers. For example, the Cianorte nest of the dark gray group and the Terra Boa nest of the black group are separated by only 15.34 km. Thus, it is believed that this division of groups is caused by behavioral differences in migration or gene flow from males.
The PCR-RFLP technique used in the population of Rondônia, Paraná and São Paulo and obtained values of molecular variance that indicate high genetic differentiation among the population with interpopulation variation equal to 83% and intra-population equal to 17% (Santos et al., 2015). However, this marker was efficient in identify the matrilineage in the analyzed nests and observe the occurrence of hybridization between T. angustula and T. fiebrigi, suggesting incomplete specific separation.
Results obtained by Bronzato (2011), with bees of T. angustula and T. fiebrigi indicated that it is not possible to separate these two groups into Tetragonisca species only in population, which suggests that hybridization can occur between these groups. Francisco et al. (2014) found reciprocal monophyly in mtDNA of T. angustula and T. fiebrigi but with evidence of hybridization with the microsatellite data, suggesting that the divergence time between them is insufficient to prevent gene flow, similar results to the present study.
Disagreement on identification using of the evaluated molecular and morphological markers is probably related to speciation that may have occurred between the two species and a secondary contact before a complete separation among them.
The secondary contact may have occurred by human activities on these bees, such as deforestation and the use of these bees by beekeepers because of the ease of management, honey production with good acceptance and high prices in the market.
It is believed that the greater number of interactions exhibited by isozymes (greater thickness on the upper bars corresponding to these markers), is because more than 54% of the subjects included in this analysis came from the study of these markers, while 29% comes from analysis with microsatellites and 17% by PCR-RFLP. Besides, the study by Santos et al. (2015), all haplotypes of PCR-RFLP marker used showed interaction with individuals analyzed and classified as hybrids.
The network of interactions between bees of T. angustula and T. fiebrigi showed that it is not possible to separate the two species and even the bees identified as hybrids. Therefore, these results corroborate the findings of studies by Baitala et al. (2006) and Francisco et al. (2014), but Santos et al. (2015) showing that these bees have not formed two species yet. The differences can be observed in their morphology, geographical distribution (more restricted in T. fiebrigi) and various molecular markers already analyzed.

Conclusion
The morphology of the coloration of the mesetristernum of Tetragonisca bees is not a good marker for the taxonomic identification of these bees, as they are not genetically differentiated according to the coloration, which is believed to be a form of polymorphism. None of the molecular markers analyzed was efficient to characterize and definitively separate T. angustula and T. fiebrigi as two species. The nuclear genetic variability in the populations analyzed was low.