DNA Barcoding of Stingless Bees (Hymenoptera: Meliponini) in Northern Peruvian Forests: A Plea for Integrative Taxonomy

: Stingless bees (Hymenoptera: Meliponini) are among the most important pollinators of tropical forests. Peru is considered a hotspot of biodiversity of Meliponini, but many areas of this country (e.g., Peruvian Amazon) remain unexplored. We aimed to produce a ﬁrst inventory of stingless bee species dwelling in humid and seasonally dry forests of northern Peru by combining traditional (morphologically-based) taxonomy and DNA barcoding. Specimens were collected in 2020 at ﬁve sites located in San Martin and Piura regions. We identiﬁed 12 genera of Meliponini. Among those, Trigona and Plebeia were the most abundant (45.9% and 12.8% respectively), whereas Nannotrigona and Scaura were the least represented ones (2.3%). We assigned a reliable species identiﬁcation to about 30% of specimens ( Trigona amazonensis , T. muzoensis , T. williana , Partamona testacea , Scaura tenuis , Tetragona goettei, and Tetragonisca angustula ). Yet, more than a half of the specimens received a provisional identiﬁcation (e.g., Geotrigona cf. fulvohirta , T. cf. amalthea , T. cf. fuscipennis, T. cf. hypogea , Melipona cf. cramptoni , Partamona cf. epiphytophila , Ptilotrigona cf. perenae , Scaura cf. latitarsis , Tetragona cf. clavipes , Trigonisca cf. atomaria ). We also highlighted an extensive polyphyly that affected a number of currently recognized species (e.g., T. fulviventris , T. guianae , Plebeia franki , P. frontalis , M. eburnea , M. illota ), whose members were split into various clades. Finally, 16% of individuals failed to be identiﬁed at the species level ( Trigona sp. 1, T. sp. 2, Nannotrigona sp., Partamona sp., Scaptotrigona sp. 1, S . sp. 2, Trigonisca sp. 1, and Trigonisca sp. 2). We discuss our ﬁndings according to the current faunistic and biogeographic knowledge of Meliponini in Peru and the Neotropical region. We also remark on the importance of conducting a taxonomic revision of stingless bees and improving both their morphology-based identiﬁcation keys and BOLD repository. Finally, we claim that integrative taxonomy shall be strongly implemented to truly assess the biodiversity of Neotropical stingless bees, allowing conserving these important pollinators and the associated traditional meliponiculture in an effective manner.


Introduction
The pantropical tribe Meliponini (Hymenoptera: Apidae) counts over 500 species and represents the most diverse group among corbiculate bees [1]. These eusocial bees, also known as stingless bees for their atrophied stings, are the main pollinators of tropical forests [2]. They are numerically dominant within bee communities in the Neotropical region, dwelling in forests, savannahs, or other habitats with open vegetation [3][4][5].
Before the arrival of Apis mellifera in the New World, stingless bees played an important role in pre-Columbian civilizations, providing honey, wax, and pollen [6]. Nowadays, Meliponini are traditionally reared by native populations of South America to exploit the biological and therapeutic properties of their honeys and related products [7]. In addition, modern breeding techniques are being developed to enhance the pollination of several crops by stingless bees [8][9][10][11].
With about 175 species, Peru is considered a hotspot of biodiversity of Meliponini in the Neotropics [12]. However, this is likely a conservative estimate of the real biodiversity of Meliponini in the country, because a comprehensive faunistic database of Peruvian stingless bees is far from being complete. In fact, although some inventories have been compiled in the past [12][13][14][15][16][17][18], many areas of the country have been scarcely explored and many aspects of their biology and ecology remain to be understood [19]. Moreover, species recognition of stingless bees is often hindered by the lack of dichotomous keys or any other valuable sources to attain a correct species identification (e.g., illustrated keys with high-quality images) in understudied genera and species complexes.
In the last decades, DNA barcoding [20,21] has represented a genetic tool of paramount importance to taxonomically identify those bee species [22][23][24][25] for which dichotomous keys were unavailable and/or whose morphologically-based identification proved challenging [26,27]. Thus far, most of the sequences from the barcode fragment of cytochrome oxidase I (COI) of Meliponini have been acquired to support morphometric data, resolve some specific taxonomic issues [28][29][30][31][32][33][34][35][36][37], or conduct phylogenetic and phylogeographic studies [38,39]. These sequence data, together with those acquired from museum stingless bee specimens, currently populate the barcode reference database BoldSystems, which counts over 2500 DNA sequences from Meliponini specimens, for a total of about 350 recognized species. This referenced sequence collection may offer a valuable background to help species identification of stingless bees by non-experts, but the robustness of this repository and its efficacy in truly helping the taxonomic identification of newly-collected specimens has never been thoroughly evaluated. In fact, no taxonomic survey has been carried out so far in the Neotropical region by using DNA barcoding to assess the biodiversity of Meliponini locally.
Through comparisons of morphologically-and DNA barcoding-based identification, we aimed to produce a first inventory of stingless bee species dwelling in humid and seasonal dry forests of San Martin and Piura regions in Northern Peru. In doing so, we evaluated the state-of-the-art and usefulness of combining traditional morphologicallybased taxonomy and DNA barcoding in diagnosing field-collected individuals. We also highlighted an extensive COI polyphyly in a number of morphologically-recognized species of stingless bees, which encourages future research to solve the taxonomy of some groups within this tribe in the Neotropics and deepen their biogeographic history.

Study Sites and Stingless Bee Sampling
Specimens of stingless bees were collected in July-November 2020 in 5 [41]. Three sites were located in rainforests: JP (800-110 m a.s.l. and −6 • [42,43]. The greatest drought of MA is due to the influence of subtropical atmospheric subsidence and the rain-shadow effect of the Andean Cordillera, which prevents the Amazon humid air from reaching the Pacific coast [43]. At each site, we walked 5 random transects of 1 km each, catching stingless bees both with direct and indirect methods. The direct method consisted of capturing flying bees with sweep nets on flowers, near water sources, or from wild colonies encountered 5 m right or left of the transect line [44]. For the indirect method, we installed 4 baited stations every 200 m along the transect, consisting of plates containing rotten fish, fermented fruits, and dog feces [12]. Furthermore, at the baited stations we sprayed on vegetation (30 to 100 cm above the ground; 15 sprays) a solution of honey and water in the ratio of 1:2 (v:v) honey: water with salt (NaCl) [45]. Stingless bees were captured with sweep nets at the baited stations 60 min after their installation. Sampling, conducted on sunny and windless days by 4 operators, started at 8:00 a.m. and ended up at 4:00 p.m. All specimens were preserved in 96% ethanol for further analysis.

Morphological and DNA Barcoding Data Collection
Stingless bees were identified based on their morphology according to dichotomous keys, when available, and/or by relying on diagnostic photographs [18,[46][47][48][49][50][51][52][53][54][55][56][57][58]. A Zeiss Axio Zoom V16 microscope was used for the observation of diagnostic morphological characters. Multifocal digital images of the specimens were acquired by the software ZEN 2.3 pro. The software HELICON FOCUS 6.7.1 was used to combine stacking images taken at different focus distances.
COI barcode fragments were generated from each collected specimen by extracting the total genomic DNA from one middle leg following the "salting out" procedure [59]. Extracted DNA was then eluted in 100 µL H2O milliq and stored at −20 °C. PCR was performed to amplify the 650 bp barcode fragment of COI by using primers LCO1490 (5′-GGT CAA CAA ATC ATA AAG ATA TTGG-3′) and HCO2198 (5′-TAA ACT TCA GGG TGA CCA AAAAAT CA-3′) [60]. PCR was conducted in a total reaction volume of 25 µL, containing 0.5 pmol of each primer, 10 mM Tris-Cl, pH 8.3 and 50 mM KCl, 1.5 mM MgCl2, 2.5 mM dNTPs, 2 µL of the DNA template, and 1 unit of Taq DNA polymerase (Meridian, Douglas, CO, USA). PCR cycling conditions consisted of an initial denaturation of 3 min at 94 °C followed by 35 cycles of 30 s at 94 °C, 30 s at 50 °C, and 1 min at 72 °C, with a final elongation step of 10 min at 72 °C. Products were visualized on a 1% agarose gel stained At each site, we walked 5 random transects of 1 km each, catching stingless bees both with direct and indirect methods. The direct method consisted of capturing flying bees with sweep nets on flowers, near water sources, or from wild colonies encountered 5 m right or left of the transect line [44]. For the indirect method, we installed 4 baited stations every 200 m along the transect, consisting of plates containing rotten fish, fermented fruits, and dog feces [12]. Furthermore, at the baited stations we sprayed on vegetation (30 to 100 cm above the ground; 15 sprays) a solution of honey and water in the ratio of 1:2 (v:v) honey: water with salt (NaCl) [45]. Stingless bees were captured with sweep nets at the baited stations 60 min after their installation. Sampling, conducted on sunny and windless days by 4 operators, started at 8:00 a.m. and ended up at 4:00 p.m. All specimens were preserved in 96% ethanol for further analysis.

Morphological and DNA Barcoding Data Collection
Stingless bees were identified based on their morphology according to dichotomous keys, when available, and/or by relying on diagnostic photographs [18,[46][47][48][49][50][51][52][53][54][55][56][57][58]. A Zeiss Axio Zoom V16 microscope was used for the observation of diagnostic morphological characters. Multifocal digital images of the specimens were acquired by the software ZEN 2.3 pro. The software HELICON FOCUS 6.7.1 was used to combine stacking images taken at different focus distances.
COI barcode fragments were generated from each collected specimen by extracting the total genomic DNA from one middle leg following the "salting out" procedure [59]. Extracted DNA was then eluted in 100 µL H2O milliq and stored at −20 • C. PCR was performed to amplify the 650 bp barcode fragment of COI by using primers LCO1490 (5 -GGT CAA CAA ATC ATA AAG ATA TTGG-3 ) and HCO2198 (5 -TAA ACT TCA GGG TGA CCA AAAAAT CA-3 ) [60]. PCR was conducted in a total reaction volume of 25 µL, containing 0.5 pmol of each primer, 10 mM Tris-Cl, pH 8.3 and 50 mM KCl, 1.5 mM MgCl2, 2.5 mM dNTPs, 2 µL of the DNA template, and 1 unit of Taq DNA polymerase (Meridian, Douglas, CO, USA). PCR cycling conditions consisted of an initial denaturation of 3 min at 94 • C followed by 35 cycles of 30 s at 94 • C, 30 s at 50 • C, and 1 min at 72 • C, with a final elongation step of 10 min at 72 • C. Products were visualized on a 1% agarose gel stained using Midori Green Advance dye (Nippon Genetics Europe, Düren, Germany). PCR products were purified using the ExoSAP-IT™ PCR Product Cleanup Reagent (Applied Biosystem) and sent to the sequencing facility of Microsynth AG (Balgach, Switzerland).
Sequences were deposited in GenBank under Acc. n • OP093417-OP093549. All specimen vouchers were deposited in Estudios Amazonicos Biological Material Depositary Center (Tarapoto, Peru).

Species Identification
COI sequences were edited and aligned with Staden Package 2.0.0b11-2016 (http: //staden.sourceforge.net/ accessed on 1 February 2021. The identification tool in BOLD (http://www.boldsystems.org/index.php/IDS_OpenIdEngine, accessed on 1 February 2021) was used to taxonomically identify our specimens. Their attribution to species was based on Kimura 2-parameter (K2P) corrected genetic distances [61] with a threshold of 3% as a cut-off value [20]. Refined single linkage analysis (RESL) in BOLD [62] was used to obtain the BINs corresponding to the taxonomic categories. For each identified stingless bee genus, we aligned the COI sequences obtained from our sample with those available in BOLD and NCBI repositories. To do so, we generated a non-redundant database of COI sequences archived in public repositories of GenBank and BOLD for each genus of Meliponini (Geotrigona n = 1; Melipona n = 182; Nannotrigona n = 1; Scaptotrigona n = 136; Partamona n = 101; Plebeia n = 102; Ptilotrigona n = 5; Scaura n = 11; Tetragona n = 20; Tetragonisca n = 1469; Trigona n = 100; Trigonisca n = 3). To assess the relationships among our query sequences and their neighboring reference sequences, we built consensus (distancebased) neighbor joining (NJ) phylogenetic trees (500 bootstrap replicates, K2P model, midpoint rooted) using MEGA 11 [63]. Barcode-phylogenetic trees generated for the various stingless bee genera were only built to evaluate the clustering of sequences from the same (predicted) species in monophyletic clades (not to trace evolutionary histories).
The final attribution of the species name to each individual was achieved by relying on its identification obtained through BOLD (BOLD ID) and morphology (MORPHO ID), and by considering the relative position of analyzed specimens within the comprehensive genus-level phylogenetic trees (see Table 1 caption for further details). Table 1. List of collected Peruvian stingless bees with their final identification obtained by integrating DNA barcoding, morphology, and inspection of COI phylogenetic trees. Final taxonomic IDs were established as follows: a Specimens ascribed to a species with BOLD ID % similarity > 97, confirmed by morphology and with sequences included in a monophyletic clade with all conspecifics (with node score > 70) received the "BOLD ID attribution"; b Specimens ascribed to a species with BOLD ID % similarity > 97, confirmed by morphology, but whose sequences were interspersed but grouped in various monophyletic clades (node score > 70) embedding conspecifics were named as "cf. + BOLD ID attribution" or "cf. + BOLD attribution + clade ID" when referring to clades grouping only Peruvian specimens; c Specimens unidentified at the species level in BOLD, but belonging to a clade (node score > 70) including at least one identified ( a, b ) specimen at the species level, and confirmed by morphology, were named as "cf. + BOLD ID attribution of the sister taxon"; d Specimens with BOLD and/or morphologically-based species attributions, not examined on reconstructed phylogenetic trees due to lack or strong underrepresentation of BOLD data, received a final (but uncertain) taxonomic ID provided by authors; * Unidentified specimens at species level in BOLD and by morphology, forming an independent clade in the phylogenetic tree were named as "genus + sp.".

SAMPLE CODE
BOLD ID (% Similarity)

Results
COI barcode fragments were obtained from a total of 133 specimens of Peruvian stingless bees. Except for five specimens (JP29, JPA0015-1, JPA0016-1, POA44, and MA51), molecular data attributed everyone to the correct stingless bee genus (Table 1). Among the 12 identified genera, Trigona and Plebeia were the most abundant, constituting 45.9% and 12.8% of the total sample, respectively, whereas Nannotrigona and Scaura were the least represented (2.3%) (Supplementary Material Figure S1).
The relative presence of Meliponini genera varied across sampling sites, with some of them found exclusively in a single site (Nannotrigona in MA and Partamona in JP), while some others were specific to certain forest ecosystems, e.g., rainforests (Melipona in JP, POA, and PY, and Scaura in JP and POA) ( Figure S1).
By integrating molecular and morphological data, we assigned a reliable speciesspecific identification to 40 specimens (30.1% of the whole sample). The identification was provisional for 72 individuals (54.1%), but likely to be definitive after comparing our specimens with reference material or consulting a specialist of the taxon (=cf.). Finally, the ascription of 21 individuals (15.8%) at the species level remained uncertain (=sp.) ( Table 1). The final identification obtained for individuals identified within each stingless bee genus is detailed in Table 1 (see also results of SDMs in Table S1).

Genus Trigona
Approximately one-third of specimens (n = 20; 33%) of Trigona were determined unambiguously at the species level (Table 1). Based on genetic (BOLD % similarity = 98.5) and morphological data, we consistently ascribed 10 individuals (clustering in a monophyletic clade; bootstrap = 69) to T. muzoensis. In addition, six and four individuals were consistently assigned to T. amazonensis and T. williana, respectively, as they showed (i) high matching scores with conspecific entries in BOLD (% range of similarity = 98.6-99.3 and 97.9-98.6, respectively) ( Table 1), (ii) an intra-specific phylogenetic cohesion (bootstrap = 100) (Figure 2), and (iii) a morphological identification (52) coherent with barcode determination ( Figure S2). Unfortunately, public sequence data were not available for T. amazonensis and T. williana in BOLD and, thus, an assessment of their monophyletic condition remained unexplored.
Most Trigona individuals (n = 37; 61%) were attributed with a certain degree of uncertainty (cf.) to previously described species (Table 1). First, we cautiously assigned eight specimens to T. cf. fuscipennis because of an apparent polyphyly for this taxon: in fact, our specimens grouped in a monophyletic clade (bootstrap = 92) including all Costa Rican BOLD conspecifics, except ASINH876-12 (Costa Rica), which, however, appeared distantly related (Figure 2). In any case, all specimens showed diagnostic features of T. fuscipennis [46] (Figure S2).
Seven other individuals were ascribed to T. cf. hypogea: sequences from these specimens showed a phylogenetic cohesion in the phylogenetic tree (bootstrap = 94) (Figure 2), but four individuals were unidentified by BOLD, whereas the other three received a 97% similarity score with T. hypogea (Table 1). All these seven individuals are morphologically similar, showing the diagnostic characters [48], except for the bristles on the surface of scape that appeared slightly longer ( Figure S2).  which were further split into two distinct sub-groups, not significantly related to each other (T. cf. fulviventris Clade A and Clade B; Table 1, Figure 2). All specimens from MA showed jaws with four teeth and a predominantly black body with the typical rust/yellowish coloration of T. fulviventris abdomen (46) (Figure 3). On the other hand, the single specimen from UT (UT11) showed the typical elongated black-colored abdomen of T. guianae, hyaline wing membrane, brown C and R veins, light brown pterostigma, and brown microtrichia [52,58] (Figure 3). Coherently, UT11 was identified as T. guianae in BOLD (98.9%) (and also appeared as a well-delimited OTU based on the distance-based SDM; Table S1). All individuals from POA with a black abdomen formed the sister-group of UT11 ( Figure 2) and, thus, were assigned to T. cf. guianae Clade A (Table 1; Figure 3). PY individuals also showed a blackish abdomen and formed an additional, but unrelated, group of T. guianae, and hence were assigned to T. cf. guianae Clade B (Figure 2). Apparently, the specimens belonging to the two clades do not show evident morphological differences (Figure 3).

Genus Plebeia
All specimens (n = 17) ascribed to Plebeia received an uncertain taxonomic identification. However, we provisionally ascribed our Peruvian specimens to previously described species (= cf.) waiting for confirmation (Table 1).
In particular, seven specimens (unidentified at the species level by BOLD; Table 1) both from JP (JPA002, JP5, JP16, and JP45) and POA (POA8, POA32, and POA34) clustered in a monophyletic group (bootstrap = 99%), including the BOLD entries SICOB7718-11 (Panama) and GBAH6196-09 (unknown origin) of P. franki. Within this group, our specimens were split into two well-supported but distinct clades (bootstrap = 99 and 83, respectively) ( Figure 4). Since morphology did not help identifying these individuals, we provisionally ascribed these specimens to P. cf. franki Clade A or Clade B, respectively. The specimens belonging to Clade A and Clade B do not show distinctive morphological characters ( Figure 5). We only observed some polymorphism (but not correlated to the Finally, four Trigona individuals identified by BOLD at the genus level only (Table 1) formed two lineages that were unrelated to all other recognized groups ( Figure 2) within Trigona: one included three specimens from JP (JPA009-1, JP10, and JP21), which, hence, were attributed to Trigona sp. 1; the other was formed by a single individual from PY (PY8), which, according to morphology [52,58], could be attributable to T. dallatorreana, but here provisionally named Trigona sp. 2 ( Figure 3). This latter, in particular, based on ASAP (Table S1), results in an OTU separated from all other BOLD entries.

Genus Plebeia
All specimens (n = 17) ascribed to Plebeia received an uncertain taxonomic identification. However, we provisionally ascribed our Peruvian specimens to previously described species (= cf.) waiting for confirmation (Table 1).
In particular, seven specimens (unidentified at the species level by BOLD; Table 1) both from JP (JPA002, JP5, JP16, and JP45) and POA (POA8, POA32, and POA34) clustered in a monophyletic group (bootstrap = 99%), including the BOLD entries SICOB7718-11 (Panama) and GBAH6196-09 (unknown origin) of P. franki. Within this group, our specimens were split into two well-supported but distinct clades (bootstrap = 99 and 83, respectively) ( Figure 4). Since morphology did not help identifying these individuals, we provisionally ascribed these specimens to P. cf. franki Clade A or Clade B, respectively. The specimens belonging to Clade A and Clade B do not show distinctive morphological characters ( Figure 5). We only observed some polymorphism (but not correlated to the two groupings), as JP5, POA32, and POA34 differed from JP16 and POA8 for the coloring of the abdomen and wing veins, as well as for the shape of abdomen and the pterostigma ( Figure 5). frontalis Clade B, including samples from JP, POA, and PY, which, however, were somehow morphologically similar; (iii) P. cf. frontalis Clade C, including a single specimen from POA (POA40); and (iv) P. cf. frontalis Clade D, including a single specimen from UT (UT58) (Figure 4). Although all examined Peruvian specimens showed the diagnostic morphological features of P. frontalis [47], we observed in POA40 and UT58 some evident phenotypic differences, such as a yellow abdomen with dark basal bands of the tergites, a yellow coloration of the first two pairs of legs, and some yellow spots in the mesepisterna. Furthermore, the yellow spotting on the posterior tibias extends well beyond the basal part ( Figure 5).
Peruvian specimens (n = 5) of Melipona (PY23, PY56, JPA0012, JP48, and POA42) were scattered in an overall poorly resolved phylogenetic tree ( Figure S3). JP48, POA42, and PY23 were highly similar (99.84-100% similarity) to BOLD private entries of M. eburnea, and, consistently, these owned morphological features typical of the species [18] ( Figure  S2). Noteworthy, the first two individuals (JP48 and POA42; M. cf. eburnea Clade A) were phylogenetically split apart from the latter (PY23; M. cf. eburnea Clade B) ( Figure S3), questioning the actual monophyly of M. eburnea. On the other hand, the morphological inspection of JPA0012 and PY56 did not lead to a reliable species diagnosis ( Figure S2), and, therefore, these individuals were provisionally attributed to M. cf. cramptoni (for a 99.33% of similarity with a private entry) and M. cf. illota (for a 99.83% of similarity with a BOLD specimen, though placed elsewhere in the phylogenetic tree; Figure S3), respectively. All other Plebeia were identified by BOLD as P. frontalis (% similarity = 98.84-99.84) ( Table 1). However, both BOLD entries and sequences obtained from our specimens attributed to this species appeared scattered throughout the phylogenetic tree, revealing an extensive polyphyly for P. frontalis (Figure 4). Our Peruvian sample was split into four main different lineages: (i) P. cf. frontalis Clade A, grouping all specimens from MA (bootstrap = 99); (ii) a poorly supported group (bootstrap = 59), provisionally named P. cf. frontalis Clade B, including samples from JP, POA, and PY, which, however, were somehow morphologically similar; (iii) P. cf. frontalis Clade C, including a single specimen from POA (POA40); and (iv) P. cf. frontalis Clade D, including a single specimen from UT (UT58) (Figure 4). Although all examined Peruvian specimens showed the diagnostic morphological features of P. frontalis [47], we observed in POA40 and UT58 some evident phenotypic differences, such as a yellow abdomen with dark basal bands of the tergites, a yellow coloration of the first two pairs of legs, and some yellow spots in the mesepisterna. Furthermore, the yellow spotting on the posterior tibias extends well beyond the basal part ( Figure 5).
Peruvian specimens (n = 5) of Melipona (PY23, PY56, JPA0012, JP48, and POA42) were scattered in an overall poorly resolved phylogenetic tree ( Figure S3). JP48, POA42, and PY23 were highly similar (99.84-100% similarity) to BOLD private entries of M. eburnea, and, consistently, these owned morphological features typical of the species [18] (Figure S2). Noteworthy, the first two individuals (JP48 and POA42; M. cf. eburnea Clade A) were phylogenetically split apart from the latter (PY23; M. cf. eburnea Clade B) ( Figure S3), questioning the actual monophyly of M. eburnea. On the other hand, the morphological inspection of JPA0012 and PY56 did not lead to a reliable species diagnosis ( Figure S2), and, therefore, these individuals were provisionally attributed to M. cf. cramptoni (for a 99.33% of similarity with a private entry) and M. cf. illota (for a 99.83% of similarity with a BOLD specimen, though placed elsewhere in the phylogenetic tree; Figure S3), respectively.
Morphology [15], BOLD ID (99.52% similarity), and tree reconstruction (bootstrap = 91) were coherent in assigning nine specimens (JP) to Partamona testacea (Table 1, Figures S2 and S3). On the other hand, two specimens from the same locality (JP4 and JPA005) were genetically related to a private BOLD entry (99.84% similarity) of P. epiphytophila. Although this identification was confirmed based on morphology [15] (Figure S2), the lack of public sequences in BOLD obliged us to cautiously attribute these individuals to P. cf. epiphytophila. Finally, a single individual, JPA006, showed a peculiar position within the phylogenetic tree (placed as a sister group of P. testacea + P. mulata) ( Figure S3), with no significant affinities with any particular taxon (neither based on BOLD ID). Therefore, this individual was preliminarily attributed to Partamona sp. ( Figure S2). Noteworthily, all SDMs pointed to JPA006 as an independently evolving lineage/species (Table S1).
Ptilotrigona individuals collected in JP (n = 2) and POA (n = 3) were genetically attributed to P. pereneae (99.17% similarity). However, too few sequences of Ptilotrigona were available in BOLD to provide a satisfactory phylogenetic reconstruction, and since the morphological distinction between P. pereneae and P. lurida can be somewhat ambiguous [49] ( Figure S2), we ascribed these specimens to Ptilotrigona cf. perenae (Table 1).
Four individuals, unambiguously identified through morphology as Scaptotrigona sp., were erroneously attributed to Oxytrigona sp. by BOLD (Table 1). Indeed, these specimens clustered within the phylogenetic tree reconstructed for Scaptotrigona ( Figure S3), but separately with respect to all BOLD entries available for this genus. In particular, we ascribed JPA0015-1, JPA0016-1, and POA44 to Scaptotrigona sp. 1 and the slightly diverging specimen MA51 to Scaptotrigona sp. 2 ( Figure S2).
Among the three specimens of Scaura, only JP8-1 could be clearly ascribed to Scaura tenuis based on morphology [52] (Figure S2), BOLD identification (98.53% similarity), and tree reconstruction ( Figure S3). The remaining two individuals (JPA004 and POA31) were attributed to Scaura cf. latitarsis, as they clustered with other BOLD entries of the same species (although this same clade also embedded BOLD entries of another taxon, i.e., S. argyrea; Figure S3). Based on the species description [66], morphological analyses also confirmed this taxonomic attribution ( Figure S2).
Tetragona specimens were split into three groups ( Figure S3). Morphology [58] ( Figure S2) and BOLD identification were coherent in considering JP12 and POA25 as members of T. clavipes (Table 1), but, due to the low bootstrap support (63%), we provisionally ascribed these specimens to T. cf. clavipes. Specimens identified by BOLD as T. mayarum (95.46-95.67% similarity with BOLD private entries) comprised a well-supported clade (bootstrap = 100), sister to T. ziegleri. Since morphological analyses did not allow confirming their genetic identification, we decided not to assign species-specific naming to these specimens (Tetragona sp.) ( Figure S2). Finally, specimens ascribed to T. goettei (99.51-100% BOLD similarity) formed a monophyletic cluster ( Figure S3), and the attribution at the species rank was also confirmed on a morphological basis [52] (Figure S2).
As for Trigonisca, JP29 unequivocally belonged to this genus (Trigonisca sp. 1) but was erroneously attributed by BOLD to Plebeia (Table 1, Figure S2). On the other hand, MA9-1 (99.64% in BOLD) was clearly attributed to Trigonisca atomaria: however, this "whollyblack" individual of T. atomaria ( Figure S2) diverged morphologically from the partially reddish-colored of a (seemingly) conspecific specimen from Tumbes dry forest [16]. Since dichotomous keys to recognize this species are not available, we provisionally considered this individual as T. cf. atomaria (Table 1, Figure S2). Finally, it was not possible to attribute a taxonomic identification for specimen MA2, which, hence, was named as Trigonisca sp. 2 (Table 1, Figure S2).

Discussion
By comparing morphologically-and molecularly-based species diagnoses, we provided the first faunistic inventory of stingless bees inhabiting humid and seasonally dry forests of San Martin and Piura regions, which have been, thus far, poorly explored for Meliponini biodiversity.
Overall, with few exceptions, both morphological and DNA barcoding approaches were coherent at identifying stingless bee genera. In fact, DNA barcode-based misidentifications at the genus level only affected the recognition of five specimens, which, based on morphology, were unambiguously ascribed to Nannotrigona, Scaptotrigona, and Trigonisca (Table 1). Although morphology alone allows an easy discrimination of genera of Neotropical Meliponini, these extrinsic errors of the DNA barcode tool due to the inclusion of wrongly identified sequences into the reference database [67] may accumulate over time and impede a correct genetic identification of stingless bees. These errors shall be corrected by double-checking deposited sequences and associated metadata (e.g., type of voucher, geographical coordinates, images, etc.) to avoid major mistakes in taxonomic identifications.
Beyond the few inconsistencies observed, all identified genera were largely expected based on the faunistic checklists available for Peruvian Meliponini [12][13][14][15][16][17][18]. Among them, Trigona and Plebeia were the most abundant, ubiquitously found in all five localities (Table 1), which is consistent with the high speciosity and abundance of these two Meliponini genera in the Neotropics [68]. On the other hand, despite being the largest Neotropical stingless bee genus-counting about 74 species [68]-Melipona was not abundant in our collection (Table 1). Although disconcerting, the scarcity of Melipona spp. could be due to the susceptibility of these species to deforestation and its restriction to well-preserved areas and particular habitat types of relatively high quality [12,69,70].
Other genera occurred exclusively at a single site, such as Partamona and Nannotrigona. Their limited localization seems odd as these two genera are very common and widely distributed in the Neotropics [15], even in urbanized and disturbed areas [54]. Other genera (e.g., Melipona, Scaura, Ptilotrigona) were recorded in humid ecosystems only, but this might simply reflect undersampling in seasonally dry forests during a period when flowering is limited. In addition, given that nesting biology is highly diverse and unknown in many species of Meliponini, the scarce availability of nesting substrates (e.g., some Scaura species need both termite mounds and pre-existing cavities for nesting [58]) could also have limited the presence of certain genera to some areas.
These inconsistencies could be explained by a difficulty of sampling elusive genera, either because of their limited geographic distribution in the Neotropics (e.g., Noguerapis), or a low local abundance despite their extensive geographic ranges (e.g., Celetrigona, Cephalotrigona, Aparatrigona, Leurotrigona, and Schwarzula) [51]. Additional faunistic surveys will be needed to assess the abundance and actual distributional patterns of the observed (or so-far-undetected) Meliponini genera and confirm their possible peculiar localization in different biomes of the northern Peruvian forests.
In general, morphology-based species identification matched with those provided by BOLD (Table 1), which, hence, appears as an affordable tool to recognize stingless bee taxa accepted by the current taxonomy. However, many collected specimens were not recognized at the species level (these were just ascribed to genera Trigona, Tetragona, Partamona, Trigonisca, and Scaptotrigona) and, hence, their identity shall be ascertained.
Most of the identified species were previously reported in Peru [12,[16][17][18], except for Melipona cf. cramptoni, Plebeia cf. frontalis, and P. cf. franki (Table 1). In fact, M. cramptoni (a "dark form" of M. fulva Lepeletier, 1836 [51]) is known to be distributed in Colombia, Guyana, Venezuela, and Brazil; P. franki in Costa Rica, Panama, and Colombia, whereas P. frontalis in Central America (e.g., in Mexico [47]) and to a lesser extent in South America, where it seems southward limited to Colombia [51]. It is worth noting, however, that a comparison of P. cf. frontalis from MA ( Figure 5) with a photographed, and phenotypically similar, individual of Plebeia sp. of the Tumbesian dry forest [16] would allow confirming if the occurrence in Peru of this taxon was already recorded.
Comprehensive phylogenetic reconstructions at the genus level were mostly characterized by weakly supported nodes and/or extensive polyphyly (Figure 2, Figure 4 and Figure S3). In particular, Trigona and Plebeia (Figures 2 and 4), but also Melipona and Scaura ( Figure S3), showed incompletely resolved COI phylogenetic trees. Consequently, the accuracy of DNA barcoding to correctly identify Neotropical Meliponini appears hindered by a lack of reliable "barcoding gaps" (i.e., the separation between intraspecific variation and interspecific divergence estimated on the basis of the mitochondrial COI fragment) delimiting species boundaries (an intrinsic error of the barcoding method [67]). This raises concerns about the robustness of diagnostic morphological characters currently used to distinguish Neotropical Meliponini at the species rank.
In addition, the lack of public sequences in BOLD for some species (e.g., T. amazonensis, T. muzoensis, and T. williana) do not allow evaluating their genetic cohesion, which, then, should be thoroughly assessed in the next future.
Concurrently to a poorly-established taxonomy, the pervasive intraspecific gene pool fragmentation could be the result of paleoclimatic and geological (e.g., Andean uplift; see [50,51]) events that are known to have profoundly contributed to the origin and evolution of Neotropical biodiversity [71]. The rise and fall of biogeographic barriers and, in some cases, bee-plant specialization [72], likely triggered isolation at both large (e.g., continental) and small (e.g., regional) geographical scales [73], favoring a rapid genetic diversification of Meliponini.
For instance, vicariance events due to the Andean uplift likely promoted the subdivision of the species pair distributed on either side of the Andean mountain range, i.e., the Amazonian T. amalthea and the Central American/Western Ecuadorian T. silvestriana [74]. Although both taxa seemed to lack a genetic intraspecific coherence (possibly due to extrinsic errors affecting some BOLD specimens) (Figure 2), all Peruvian individuals were recognized as T. cf. amalthea, consistently with their localization (Table 1, Figure 2). A large-scale isolation-by-distance (IBD) process caused by barriers to dispersion along the Neotropical region (lasting until at least 50 Ma between South and North/Central America [75]) could explain the observed divergence between the Central American clade of T. cf. fulviventris (Belize, Costa Rica, Mexico, and Honduras) and the Peruvian lineage ( Figure 2). A regional IBD would instead account for the split of T. cf. guianae into two distinct but geographically coherent clades (Figure 2) (i.e., Clade A: POA+UT; Clade B: PY). These past events, that likely limited gene flow, further hamper the resolution of taxonomic uncertainties already challenging the systematics of T. fulviventris and T. guianae, as well as the boundaries of their actual distribution [51,76]. T. fulviventris mainly differs from T. guianae for the rust/yellowish metasome, but its abdomen can occasionally appear black-colored, as it is, yet, typical of T. guianae [46] (Figure 3). Then, although T. fulviventris is considered more widely distributed than T. guianae, their relative geographic (partially overlapping) ranges remain doubtful [51]. Concerning Peru, T. guianae has been previously recorded in humid forests of San Martin and Loreto [12,17,18]. In the Tumbesian dry forest, however, this same taxon was considered a subspecies of T. fulviventris (i.e., T. f. guianae) [16], although its nomenclature should be carefully re-examined. In fact, all Peruvian T. cf. fulviventris were collected in the seasonally dry forest of MA, suggesting a possible southward extension of the current geographic range limit of this species (placed at the southwestern sector of Ecuador [51]) to include the whole Tumbes-Piura dry forest ecoregion.
Being the second largest genus of Neotropical Meliponini [68], recognition of species of Plebeia through barcoding appeared particularly challenging. In fact, P. cf. frontalis was largely affected by polyphyly, yet apparently clustering in a geographically coherent manner ( Figure 4). The genetic peculiarity of MA specimens of P. cf. frontalis (Clade A) with respect to conspecifics from humid forests (Clade B) may reflect the outstanding biological value and elevated degree of endemism of Tumbes-Piura dry forests [77]. Members of P. cf. frontalis included in Clade C + D appeared with a different habitus ( Figure 5), and then they could belong to a different species. This, however, should be ascertained considering a thoroughly systematic revision of the genus. The barcode-based assignment of Peruvian specimens to Plebeia franki, not confirmed by morphology, is also puzzling, as this species has never been recorded in Peru [51].
Overall, our results emphasize the importance of conducting a taxonomic revision of the tribe by incorporating distinct methods to species identification of Meliponini, as recommended in integrative taxonomy. Consequently, morphology-based identification tools shall be revised, updated, and made available to non-expert taxonomists to attain reliable species diagnoses. Based on these prerequisites, the BOLD database shall be revised and expanded: in fact, genera of Neotropical stingless bees are unevenly represented, being some (e.g., Trigona, Melipona, and Partamona) more populated of public (downloadable) COI sequences than others (e.g., Nannotrigona, Trigonisca, Tetragona, Scaura, Geotrigona, and Ptilotrigona). This unequal representation of Neotropical Meliponini in BOLD originates from the application of the COI fragment to resolve, along with morphometric data, some taxonomic issues related to single genera/species, in particular to, e.g., highlight isolated reproductive units within Scaptotrigona hellwegeri [78], Melipona yucatanica [29], Melipona beecheii [30], and Tetragonula iridipennis [33]; ii) detect specific entities within (Mexican) Scaptotrigona [32] and Tetragonula [37]; and iii) resolve taxonomic ambiguities in Malagasy Liotrigona spp. [28] and Kenyan Hypotrigona spp. [36]. Thus far, only rarely were Meliponini "barcoded" to acquire knowledge on stingless bee biodiversity in some areas where faunistic monitoring was conducted [34][35][36][37].
Integrative taxonomy efforts would certainly allow scaling up the use of DNA barcoding to enlarge knowledge on biodiversity of Meliponini in more areas of the Neotropics. A combined use of the already available and tested mitochondrial and nuclear markers [74,75] would greatly improve taxonomy and identification of stingless bees. Due to their important role in pollination, this knowledge would be crucial to preserve Amazonia and avoid its degradation [79]. The Peruvian forests, especially the Amazonian ones, have suffered in the last decade a significant increase in deforestation, due to the expansion of cocoa plantations, coffee and, especially, oil palms, with a significant impact on the main rainforest regions of San Martin and Ucayali [80]. Much remains to be known about the communities of pollinators living in these forests and the consequences that their populations may suffer because of continuous habitat loss.
Beyond their importance in preserving forest ecosystems, stingless bees are also an integral part of the cultural knowledge of many indigenous peoples around the world [78]. In Peru, many native peoples use their honey for medicinal purposes [81], so the conservation of stingless bees has an additional value. Due to the lack of knowledge about the diversity and the ecology of Meliponini, any species is considered breedable in meliponiculture, so much so as to justify the uncontrolled extraction from the forest of all colonies. On the other hand, not all of them prove to be suitable, and meliponiculture could end up threatening the natural populations of bees instead of benefiting beekeepers. Melipona, Trigona, and Plebeia are among the most frequently reared genera for these purposes in Peru [81]. As already evidenced in Melipona [38], species within these genera present a large genetic diversity and an unsolved phylogeny (Figure 2 and Figure S3) possibly shaped by biogeographic events, which renders genotyping of breedable strains mandatory.
In the light of all these concerns, it is imperative that inventories of stingless bees should be extended to other areas of Peru and, at the same time, that revised species-level studies shall be promoted. Knowing the specific diversity of stingless bee communities in Peru is essential for any conservation action of these important pollinators. Hence, a reliable taxonomic identification of Meliponini is crucial to drive the rural development of New World communities and promote the maintenance of local biodiversity.