Genetic and morphometric variation in Schwarziana quadripunctata and Schwarziana mourei (Hymenoptera: Apidae: Meliponini)

Schwarziana bees are a ground-nesting stingless bee distributed in the Neotropical region. Schwarziana quadripunctata was the first described and the most studied species of this genus. Now, there are four valid species of Schwarziana bees, but it has been suggested that the diversity of this taxon may be higher, due to undescribed cryptic species. In this study, we investigated the populational diversity of S. quadripunctata using workers collected at 11 localities in Brazil (from the Northeast to South region). We also included one population of S. mourei (collected in São Paulo state, 2 nests). We analysed the bees using geometric morphometrics and molecular analyses amplifying mtDNA cytochrome oxidase I (COI) and 16S to access the diversity among the populations. From the results of geometric morphometrics, the Mahalanobis distances between S. mourei and S. quadripunctata are greater than those distances among S. quadripunctata populations. A similar scenario can also be observed looking to the phylogenetic tree generated by the molecular markers. Morphometry and molecular markers data showed significant association with geographic distance, indicating the existence of intrapopulation variation in S. quadripunctata. Our hypothesis was supported, that the populations of S. quadripunctata showed differences in haplotypic diversity. Overall, these analyses revealed a moderate level of intraspecific variation among S. quadripunctata populations and discriminated well the species S. quadripunctata from S. mourei.


Introduction
Schwarziana bees are a ground-nesting stingless bee distributed in the Neotropical region. Schwarziana quadripunctata was the first species described and the most studied species of this genus. Now, there are currently four valid species of Schwarziana bees (Melo 2015), but it has been suggested that the diversity of this taxon may be higher, due to undescribed cryptic species. Schwarziana is a stingless bee genus which is distributed in Brazil, Argentina and Paraguay (Camargo and Pedro 2013;Grüter 2020). These bees are small, having a body length of c.a. 6-7.5 mm, and mostly species nest in the ground (Fig. 1). This genus was originally described by Lepeletier (1836) as Melipona quadripunctata. However, Moure (1943) suggested a separation in a new subgenus (Schwarziana) within Trigona and later this subgenus was recognized as a different genus by Schwarz (1948). Schwarz (1948) described two subspecies, S. quadripunctata and S. bipartite, but later Camargo (1974) suggested that these two subspecies were synonyms. According to Melo (2003), the 1 3 bee diversity within Schwarziana is likely greater than has so far been recognized. The Moure catalogue published in 2013 (Camargo and Pedro 2013) recognized the two species S. quadripunctata and S. mourei, with the new species S. mourei distributed in the central region of Brazil and Paraguay. In 2015, two new species of Schwarziana were identified in Brazil, S. bocainensis and S. chapadensis (Melo 2015;Grüter 2020). Although very similar to S. quadripunctata, S. mourei shows several distinguishing features, including a completely yellow scutellum  (variable in S. quadripunctata,Fig. 1a,d), normal corbicular depression (for S. quadripunctata it extends to more than half of the tibia, Fig. 1b, ), and yellowish hairs, mainly in the region of mesepisternum (Camargo 1974).
Schwarziana bees are stingless bees (Meliponini, Apidae); this group is considered important pollinators with more than 552 species described (Grüter 2020). Despite this high number of species, it is expected that in general many more stingless bees exist in nature, but are yet to be described (Rasmussen and Cameron 2010;Grüter 2020). Among the corbiculate tribes Apini, Bombini, Euglossini and Meliponini, the stingless bees are easily recognized by their small body size and reduced wing venation (Melo 2020). Depending on the genera, stingless bees present minor morphological differences among the species, generating uncertainty over the recognition of new species (e.g., Tetragonula, Rasmussen et al. 2017). To solve this problem, some studies have used auxiliary tools to properly identify bees using morphometric or molecular studies (e.g., morphometric in Nannotrigona testaceicornis, Mendes et al. 2007; in Euglossine bees Grassi-Sella et al. 2018;both markers in Melipona beecheii Quezada-Euán et al. 2007;Francoy et al. 2011;in Plebeia remota, Francisco et al. 2008;in honeybees, Oleksa and Tofilski 2015;Ferreira et al. 2020;in Tetragonilla collina, Rattanawannee et al. 2017;in Mourella caerulea, Galaschi-Teixeira et al. 2018).
Traditional morphometry allows access of biodiversity through identifications of species using parts of the insect body, particularly parts of legs (femur, tibia, tarsus), tergite, sternite, proboscis and wing venation patterns. However, more recent studies have used geometric morphometrics due to its higher accuracy (Baylac et al. 2003). The forewing represents an ideal structure to study geometric morphometrics, because the wing landmarks are genetically inherited, and the wings have two dimensions with several homologous landmarks. Although the forewing venation patterns in Meliponini bees are discontinued (Michener 2007;Melo 2020), it is possible to mark several landmarks that can be used in geometric morphometrics. This technique, which is mainly based on the comparison of the shape obtained from the plotted landmarks, when applied to studies in bees allowed to access biodiversity and differentiate species and populations (Francisco et al. 2008;Halcroft et al. 2016;Hurtado-Burillo et al. 2016;Galaschi-Teixeira et al. 2018;Francisco et al. 2008;Francoy et al. 2008;Grassi-Sella et al. 2018) and was effective to confirm species recognition (Quezada-Euán et al. 2007;Francisco et al. 2008;Francoy et al. 2011;Galaschi-Teixeira et al. 2018). Geometric morphometric analyses have also shown good efficiency in group differentiation (Mendes et al. 2007;Francoy et al. 2009). For example, Francisco et al. (2008) analysed populations of Plebeia for mitochondrial DNA, morphometric geometrics and cuticular hydrocarbons, reinforcing the existence of two distinct taxonomic unities.
Considering the establishment of the identification of the two Schwarziana species, we aimed to study morphological and genetic differences among populations of S. quadripunctata and included one population of S. mourei (2 nests) to quantify the variation level among S. quadrifasciata populations in relation to the interspecies variation between S. quadrifasciata and S. mourei. Due to the colony-founding strategy (Dependent Colony Foundation-DFC) adopted by this stingless bee, we expect high population structuring for the mtDNA (high haplotypic diversity) in S. quadripunctata, as the sampling was conducted from different regions of Brazil. The general objective of this study was to assess the variation within the group through geometric morphometrics and mtDNA, using both techniques as auxiliary tools to characterize S. quadripunctata populations at intraspecific differentiation level.

Sampling
The present study was performed considering two described species: Schwarziana quadripunctata and Schwarziana mourei (Fig. 1). Samples were collected in 12 localities from Brazil (Table 1, Fig. 2) between the years of 2007 and 2009, using nets at the nest entrance, except for the population of Santa Catarina (SC), in which the workers were collected at two different locations during foraging for nest material (faeces, although no information about distance was recorded). The sample of S. mourei was collected from one locality (Botucatu, São Paulo state) in two nests (locality 12 in Fig. 2). We defined three different regions, Northeast, Southeast and South to be analysed based in the latitude coordinate differences (Table 1). All samples were stored in 70% aqueous alcohol solution. The identification of the species was done using the key to identify the species published by Melo (2003). Before the analysis, all samples were identified morphologically by a specialist to certify the correct identification.

Geometric morphometrics analysis
The right forewings were removed with tweezers and mounted between slide and coverslip. We aimed to use at least ten individuals per colony/location, although for some sites we had fewer samples (Table 1). These forewings were photographed with a digital camera (Leica DFC500) attached to a stereomicroscope (Leica MZ16) and then captured by the program Leica IM50 Software. Ten landmarks were plotted using the software tpsDig (version 2.04) (Rohlf 2015) (Figs. S1, S2). The tps file containing the cartesian coordinates of each plotted landmark was used as input in the MorphoJ software v1.06a (Klingenberg 2011). The landmarks were aligned using Procrustes superimposition and the data were analysed using canonical variate analysis (CVA). In short, the Procrustes alignment scales all wings to the same size, superimposing them at the mass centre of the plotted landmarks and rotating them for the best adjustment. The differences among individual wings in this final landmark configuration is then considered for the calculations of partial and relative warps, which were used as parameters for calculating distances among the studied groups or individuals, as determined by the researcher (for more details, please see Marcus 1993 andRohlf 1999). We performed tests with region (geographic region) and species establishment. The centroid sizes were obtained using MorphoJ and wing size differences among region and species were analysed using Mann-Whitney U test for the significance level (0.05) in R (version 4.0.2). Mahalanobis distances between centroids of each population origin (region) of the samples were analysed and used for pairwise discriminant analysis. Mahalanobis distances were used to construct a neighbourjoining dendrogram of the morphological similarity using MEGA 6.0 (Tamura et al. 2013).

Molecular analysis
One individual from each nest of the different regions were selected for molecular analysis. We did not perform the molecular analysis from all the sampled colonies. DNA extraction from legs was performed using the phenol-chloroform method (Sheppard and McPheron 1991). Two mitochondrial genes 16S and COI were partially amplified using heterologous primers designed for orchid bees 16SWb (Dowton and Austin 1994) and 874-16S (Cameron et al. 1992) for 16S gene and the COI gene primers (CO1-F and CO1-R) designed by Dick et al. (2004). Fragments were amplified using the polymerase chain reaction (PCR) in a total volume of 25 μL, containing 250 μM of each dNTP, 2.5 mM of MgCl 2 , 2.5 μl buffer Invitrogen 10×, 1 μM of  Table S1). For each of the gene loci and both concatenated, we did phylogenetic analysis using maximum likelihood, following the Tamura-Nei model with 500 bootstraps in the software MEGA 6.0 (Tamura et al. 2013). To assess genetic differentiation among populations, we estimated Fst in the software Arlequin 3.1 (Excoffier and Schneider 2005) using the F-statistics according to Weir and Cockerham (1984) and the genetic variation partitioning within and among populations by Analysis of Molecular Variance (AMOVA) using the distance method proposed by Kimura (1980). We correlated morphological distances using Mahalanobis distances, molecular differentiation using Fst (concatenated genes) and geographic distances from different populations averaged in kilometres. Mantel tests were performed in R 1.4.1093 with package vegan (Oksanen et al. 2020).

Results
We obtained a total of 276 wings of Schwarziana for the discrimination of the S. mourei species and populations of S. quadripunctata from distinct localities and regions (Northeast, Southeast and South) using wing geometric morphometrics. All the landmarks were found in the samples. The number of colonies sampled in the three regions were 2 for the Northeast, 19 for Southeast and 7 for the South region (Table 1, Fig. 2). Additional vouchers were deposited in the Collection Camargo at USP FFCLRP-Ribeirão Preto (registry numbers: 100479-100539). The CVA showed accuracy in separating the two studied species (Fig. 3). The two first canonical axes of the CVA explained approximately 90.74% of the variation among the individual samples (Fig. 3). This result showed that the S. quadripunctata centroids of the populations from the three regions are different, although there is a slight overlapping of the individuals, and that there is a clear separation between S. quadripunctata and S. mourei. The centroid size of the wings and the dendrogram based on the Mahalanobis distances also support CVA results (Table S2A, transformation grids at Fig. S2). Wing size, accessed by the centroid size, showed significantly higher variation among populations than wing shape (Table S2, Fig. S3). For the molecular analyses, we compared the individuals collected from 2 colonies of the Northeast region, 15 colonies from Southeast and 4 colonies from the South region (Table 1). COI and 16S sequences were aligned to assess intrapopulational variation of S. quadripunctata samples and interspecific differentiation with S. mourei based on haplotype diversity and nucleotide diversity. We amplified 611 base pairs of gene COI, with a total of 53 polymorphic sites and 9 haplotypes (Haplotype diversity Hd = 0.885; Nucleotide diversity π = 0.0218; Average number of nucleotide differences k = 13.320, and Fst = 0.918, P < 0.001). However, when analysing the alignment of amino acid sequences only 3 substitutions are not synonymous (Fig. S4). Results from AMOVA analysis with three hierarchical levels showed that 48.45% of the variation was explained by variation among regions (Northeast, South and Southeast), 43.33% among populations within the states and 8.23% within populations. For 16S, we amplified 568 pb, with a total of 14 polymorphic sites and 7 haplotypes (Haplotype diversity Hd = 0.842; Nucleotide diversity π = 0.00649; Average number of nucleotide differences k = 3.668, and Fst = 0.9581, P < 0.001). AMOVA analysis for 16S gene showed that 62.38% of the variation was explained by the variation among regions (Northeast, South and Southeast), 30.43% among populations within the states and 6.19% within populations. For the concatenated haplotypes, 1179 nucleotides were aligned and generated 9 haplotypes (Haplotype diversity Hd = 0.885; Nucleotide diversity π = 0.01445; Average number of nucleotide differences k = 16.988, and Fst = 0.9293, P < 0.001) (Tables S3, S4). Results from AMOVA showed that 51.88% of the variation represents variation among groups, 40.15% among populations within the states and 7.97% variation within populations. We constructed a Maximum likelihood tree showing the relationships of concatenated (COI plus 16S) sequences characterized in S. quadripunctata and S. mourei populations (Fig. 4). Based in the tree, there is a high relationship among the haplotypes characterized in each sampling region and as expected, haplotypes of S. quadripunctata and S. mourei haplotypes were grouped in different clusters.
For the Mantel tests, we correlated the geographic distances with morphological (Mahalanobis distances) and molecular (pairwise Fst) data. Nonsignificant pairwise correlations were seen between molecular and geographic distances (r = 0.254, p = 0.178) and between morphological and geographic distances (r = − 0.18, p = 0.63). Significant correlations were seen when comparing molecular and morphological distances (r = 0.459, p = 0.01**) and overall

Discussion
The two auxiliary tools, geometric morphometrics and molecular markers, showed some level of differentiation between S. quadripunctata populations associated to geographic distances. Although we had few samples of S. mourei, the tools discriminated well this species from S. quadripunctata. The Mahalanobis distances between S. mourei and S. quadripunctata are greater than those distances among S. quadripunctata populations, which can also be observed in the dendrogram of molecular markers. The two markers showed a high correlation with geographic distance, indicating the existence of significant intrapopulation variation in S. quadripunctata. Our hypothesis was supported, that the populations of S. quadripunctata showed differences in haplotypic diversity, although some caution is required, that we did not cover its total distribution occurrence.
Wing size, accessed here by the centroid size, showed significantly higher variation among populations compared to wing shape, as measured by Procrustes distances. Considering that wing size can be an estimator of body size (Dellicour et al. 2017), this is an expected finding as size variation in insects is dependent to the quantity and quality of food provided to the larvae (Peruquetti 2003;Campos et al. 2018). Alternatively, this variation may be associated with a differential response of these populations to specific features of their local environments (Grassi-Sella et al. 2018). If environmental factors explain most of the observed differences in size, the complexity in shape of the wing involves environmental and genetic factors subjected to stronger developmental constraints (Klingenberg 2016). In addition, shape variation may be linked to local adaptive responses to foraging behaviour and flight dynamics, to face the diversity of ecotypes that the populations are exposed, as well as differences in factors, such as climate, altitude, food availability between areas (Benítez et al. 2013).
The population of S. mourei from Botucatu (2 nests) had the highest rate of exclusive mutations (mainly substitutions) compared to the samples of S. quadripunctata, both in the COI and the 16S genes, even considering that only two samples were analysed. The populations of S. quadripunctata were quite different, although the Mahalanobis distances were smaller compared to the interspecific differentiation. For the molecular markers, the COI gene showed a greater number of haplotypes compared to 16S. Nine COI haplotypes and seven 16S haplotypes were identified, with 53 and 14 variable sites, respectively. A high haplotype diversity was observed for the COI and 16S haplotypes distribution (Hd = 0.885 and 0.842, respectively), but a low and moderate average number of nucleotide substitutions between COI and 16S haplotypes (k = 13.32 and 3.668, respectively), considering the analysed samples. The average number of nucleotide substitutions was higher when compared to those previously reported in other stingless bees, such as Partamona rustica and P. seridoensis (Miranda et al. 2016;Miranda et al. 2017). This maybe could be associated to an older origin of the genus Schwarziana relative to Partamona (Rasmussen and Cameron 2010). The populations of Bahia (Northeast) and Rio Grande do Sul (RS and Santa Catarina (Sc) (region South) were more closely related than the populations of São Paulo and Minas Gerais haplotypes. In another example, the samples from a city in São Paulo state, Cunha colony (Cunha_SP_V-1), differed from the other samples because of some substitutions (COI: bp 33; 16S: bp 235), being considered a different haplotype even compared with samples collected in the same location. In the COI gene, the samples from Carandaí city from the state of Minas Gerais (MG) presented base substitutions at sites 126 and 162, differentiating themselves from specimens collected close to them, such as the samples from Ouro Preto city (MG), both from MG state. Therefore, the nine haplotypes are shown to be exclusively located, which could explain the high haplotype diversity. Another population that showed some differences such as those mentioned above was that of Espírito Santo (ES). The Mahalanobis distances showed the samples of S. mourei and the ES population are well-separated compared to the populations of S. quadripunctata. We also noticed that the latter are quite distant from the samples of S. mourei, a fact corroborated by the distances generated in the morphometric analyses. Interestingly, for the gene 16S, the population of ES share a haplotype with the population from MG, probably due to the locality proximity. The Fst estimates were significant and AMOVA indicated a greater variation among the groups (states of origin), suggesting population structuring. We also found a positive correlation in the Mantel tests among the three analysed factors, geographical, morphological and molecular distances. By associating these results with the significant Fst value, we can suggest, according to the Migration-Drift model (Balding and Nichols 1995), an isolation by distance of these populations and consequently, a limited dispersion ability of Schwarziana.
A possible caveat that needs to be considered is that mitochondrial markers may not be a good marker to determine the genetic structure of stingless bees' populations due to two main reasons. First, colony reproduction in stingless bees is a dependent colony foundation (DFC) process, in which the production of new colonies has a temporary link between the maternal and daughter colonies (Nogueira-Neto 1954;Cronin et al. 2013). Reproductive females (new queens) are philopatric and do not disperse widely to start a new colony (Engels and Imperatriz-Fonseca 1990;Wille and Orozco 1975). The second reason is that maternal inheritance of mitochondrial markers does not allow to estimate the gene flow promoted by the dispersal of males. This sex-asymmetrical dispersal was previously reported in some social insects (Cronin et al. 2013, Dessi et al. 2022, and in euglossine bees (Cerântola et al. 2011;López-Uribe et al. 2014), which leads to an expected biogeographic consequence, higher populational genetic structuring when accessed by mitochondrial DNA compared to nuclear DNA.
As seen, the high haplotype diversity found in stingless bee populations is associated to the occurrence of unique and/or distinctive haplotypes in each population. The presence of unique and/or distinctive exclusive haplotypes in S. quadripunctata populations, as well as in other stingless bees (Miranda et al. 2016;Miranda et al. 2017;Dessi et al. 2022), suggests that an area is colonized by a low number of founder females (maternal lineages). Furthermore, the usual occurrence of unique or distinctive haplotypes in the colonies from distinct areas suggests that the progress of the occupation does not occur through the arrival of new females to these sites; instead, the growth of the populations occur through the reproduction of the few original founder colonies. Our analyses reveal that Northeast and South populations seem to be very close as was showed by the CVA analyses, unfortunately we do not have a clear explanation for this finding, given the geographic distances involved between regions (Northeast, Southeast and South). There is the possibility that microevolutionary factors, particular to some geographic areas, may influence this pattern, but also not only a single factor, such as vegetation or altitude.
Interestingly, comparing the variation observed in S. quadripunctata colonies from Cunha, a city in São Paulo State, with the variation found in colonies of Plebeia remota from the same locality (Francisco et al. 2008), it was also seen that the P. remota population from Cunha (SP) showed distinctive haplotypes. It was also suggested that during the diapause, which is shorter in Cunha, P. remota build an extra layer of cerumen in the nests (Ribeiro et al. 2003). This observation associated with differences in external morphology suggest that the Plebeia population of Cunha represents a distinct taxonomic entity compared to the other studied populations (Francisco et al. 2008). Those facts point out the need for further studies to verify if populations of other Meliponini species were under influence of some specific microevolutionary factors associated with geographic area.
In conclusion, the combination of molecular and morphological tools to assess bee biodiversity is important, attesting the variation in S. quadripunctata. We included one population of S. mourei, which is not a common studied bee. Increasing sample size of S. mourei and including samples of the two other recently described species of S. bocainensis and S. chapadensis (Melo 2015;Grüter 2020) will further highlight the diversity and evolutionary history of these Schwarziana bees.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s00040-022-00878-0. permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.