Integrative approach resolves the taxonomy of Eulaema cingulata (Hymenoptera, Apidae), an important pollinator in the Neotropics

Species delimitation is a rich scientific field that often uses different sources of data to identify independently evolving lineages that might be recognized as species. Here, we use an integrative approach based on morphometrics, COI-barcoding


Introduction
Accurate species diagnosis remains a major challenge for many groups of organisms and is often described as the "taxonomic bottleneck" (Kim and Byrne 2006). On one hand, species identification can be challenging for evolutionary units that diverged recently because they may exhibit very similar morphology making them difficult to differentiate based on morphological traits alone. Difficulty in properly identifying cryptic species leads to underestimation of species diversity. On the other hand, conspecific individuals can display discrete polymorphism in morphological characters (e.g., Quezada-Euán et al. 2015;Lepeco and Gonçalves 2018) that, if used for taxonomic identification, may lead to the erroneous description of different species.
In order to circumvent some of the limitations of taxonomic classifications exclusively based on discrete qualitatively morphological characters, integrative taxonomy incorporates evidence from multiple independent datasets such as geometric morphometrics and molecular markers to inform species recognition and species boundaries (Goldstein and DeSalle 2010;Padial et al. 2010;Schlick-Steiner et al. 2010). Basing a species description on a variety of characters from different and independent datasets is generally regarded as the best practice (DeSalle et al. 2005). When species are considered as independently evolving lineages different lines of evidence are additive to each other (de Queiroz 2007).
Orchid bees (Hymenoptera: Apidae: Euglossini) are endemic to the Neotropical region. They receive their common name because males of these bees actively collect volatile compounds from orchids, and in the process pollinate the flowers (Dressler 1982). A variety of social behaviors are found among the euglossine genera, from solitary behavior to relatively complex social interactions in some species (Cameron 2004;Faria and Melo 2020). Five genera have been described in the tribe Euglossini: Euglossa Latreille, 1802, Eufriesea Cockerell, 1908, Eulaema Lepeletier, 1841, Exaerete Hoffmannsegg, 1817, and Aglae Lepeletier and Serville, 1825, encompassing approximately 240 species (Michener 2007;Moure et al. 2012). The genus Eulaema includes the largest orchid bees, with total body size varying from 18 to 30 mm in length (Oliveira 2000;Melo 2014), and is divided into two subgenera, Apeulaema and Eulaema, which are easily recognizable. Males of Apeulaema differ from Eulaema in having whitish-yellow spots on the clypeous, parocular areas and terga integument black, the latter sometimes turning to blue-violet (Moure 2000).
Eulaema (Apeulaema) cingulata (Fabricius, 1804) has a widespread distribution (from southern Brazil to southern Mexico), where it is associated mainly with wet and dry forested areas ( Fig. 1A-C). This species is frequently sampled in faunistic inventories (e.g. Pires et al. 2013;Silveira et al. 2015;Machado et al. 2018) and reported in a wide range of ecological studies related to pollination in agroecosystems (Marques et al. 2017;Gutiérrez-Chacón et al. 2018), pollination of plants in Amazon rainforest (Martel et al. 2019;Watteyn et al. 2021), and as vectors of cleptoparasitic beetles (Rocha-Filho and Garófalo 2015).
Eulaema cingulata was initially described as Centris cingulata with no information on the number and sex of the specimens studied. The only information about these specimens is their location as "America meridionalis" (according to Moure (1960) probably Guyana) and the institution in which they were deposited, the Zoologisk Museum in Copenhagen. In the species description, the first tergum and the base of the second are described as a single tergum, which led some authors to subsequent identification mistakes related to the color pattern of these terga. Lepeletier (1841), having specimens of C. cingulata, did not recognize them as such and described them as E. fasciata, based on a female, and E. cajennensis, based on a male (Moure 1960(Moure , 2000. Oliveira (2006) described E. (Apeulaema) pseudocingulata for specimens from the Amazon Forest ( Fig. 1D-F) that show a morphologically distinct feature in the velvety area of the mid-leg of males. Males of E. cingulata have a wider posterior edge in the velvety area of the middle tibia than in E. pseudocingulata (Fig. 1B, E). In addition, the coloration of the male abdomen is darker in E. pseudocingulata than in E. cingulata. Figure 1. A, B male specimen of Eulaema cingulata from Sergipe, Brazil A body in dorsal view B velvety area in the mid tibia C geographical distribution of E. cingulata in green D, E male specimen of E. pseudocingulata from Pará, Brazil D body in dorsal view E velvety area in the mid tibia F geographical distribution of E. pseudocingulata in green.
Morphological differences between females of these two forms have not yet been found ( Moure et al. 2012) and the female of E. pseudocingulata remains undescribed. Nemésio (2009) synonymized E. pseudocingulata under E. cingulata and proposed a new species (E. marcii) for the Atlantic Forest population, based on the misunderstood conclusion that this latter form had not been found in the type locality of E. cingulata (Melo in Moure et al. 2012). Moure et al. (2012) recently placed Eulaema marcii in synonymy with E. cingulata and adopted E. pseudocingulata as a taxonomically valid species as proposed by Oliveira (2006).
Despite the distinct differences in coloration and the velvety area of the mid-leg between males of E. cingulata and E. pseudocingulata, the taxonomic status of these nominal species has not been investigated with additional datasets. A comparative phylogeographic study of E. cingulata based on mitochondrial and nuclear markers revealed a lack of structure for E. cingulata throughout the whole range of the species (López-Uribe et al. 2014). That study included a small set of specimens from both E. cingulata and E. pseudocingulata and did not indicate genetic differentiation between these species. In the present study, we incorporate morphometric data (wings and heads) and molecular data from mitochondrial and genome-wide markers to test the status of both species. We used the character of the velvety area of the mid-leg of males for species identification of the specimens. Genetic and/or morphological clustering of specimens taxonomically identified based on this morphological character would support the presence of two evolutionary units.

Morphometric data
We studied 107 specimens of E. cingulata and E. pseudocingulata from across their geographic ranges. Most of the specimens were collected in the Amazon forest where both species are sympatrically distributed ( Fig. 2A, Suppl. material 1: table S1). Using a Leica DFC 295 camera, we photographed specimens attached to a stereomicroscope Leica M205C. For the analysis of wing morphometrics, we separated the right forewing from the body at the base of the radial vein using forceps and fixed it on glass microscopy slides. The heads were photographed in frontal view. For the image analysis, we saved photographs into TPS files using the software TpsUtil 1.60 (Rohlf 2013) and identified landmarks using the software tpsDig version 2.26 (Rohlf 2006). We selected 10 landmarks for the head (Fig. 3A), and 18 landmarks on the vein intersections for the wing (Fig. 3B). We choose the landmarks based on previous studies of bees (Quezada-Euán et al. 2015;Souza et al. 2015;Costa et al. 2020).
For comparisons of the overall wing and head sizes between species, we extracted the centroid size and used a generalized Procrustes analysis to capture shape variables without the effects of orientation and position of the images. A regression analysis between size and shape was performed to quantify the effect of allometry. Afterwards, we removed this allometric effect to independently quantify shape variation. The resulting landmark configurations retained only shape information (Klingenberg 2015).
To visually compare all individuals in multivariate trait space, we used a Principal Component Analysis (PCA) using the relative Cartesian coordinates of each landmark after alignment. The shape difference between species were tested using a Discriminant Analysis followed by a leave-one-out cross-validation test (Lachenbruch 1967). All analyses were performed using the software MorphoJ (Klingenberg 2011). Additionally, we used the percentages of correct classification to evaluate the discriminatory power of wing and head shapes.

DNA sampling, extraction, amplification, and sequencing
We extracted total genomic DNA from the hind legs or thoracic muscles of specimens using Qiagen DNeasy Kit (Qiagen) with modifications to maximize DNA yield for dry specimens as incorporated in Evangelista et al. (2012). Before extraction, we placed pinned specimens in a humid chamber for 24 hours to facilitate the removal of the leg or abdomen muscle without damaging the whole specimen. For tissue digestion, we added ATL buffer and 20μL of proteinase K (20 mg/ml) for the first six hours, and 10μL additionally after six hours of incubation at 55 °C .
We amplified sequences of the cytochrome oxidase I (COI) region using universal barcoding primers LCO (5'-GGTCAACAAATCATAAAGATATTGG-3') and HCO (5'-TAAACTTCAGGGTGACCAAAAAATCA-3') (Folmer et al. 1994). Polymerase Chain Reaction (PCR) amplifications were performed in a final volume of 17μl including 0.13 μL of Taq Polymerase (Qiagen) , 2μL of genomic DNA, 1.3μL of each primer (10μM), 3.4μL of dNTPs (10mM), 0.85μL of MgCl, 1.7μL of 10× Qiagen Buffer and 6.32μL of purified water. We performed amplification with an initial step of three minutes at 94 °C, followed by 40 cycles of 30 seconds at 94 °C, 30 seconds for annealing at 50 °C and 1 min at 72 °C. After 40 cycles, we performed a final step at 72 °C for 10 minutes. We performed all steps related to DNA extraction and COI amplification at the Molecular Biology Laboratory of the Museu de Zoologia da Universidade de São Paulo, Brazil (MZSP). We sent the amplified COI fragments to Macrogen (Seoul, South Korea) for post-PCR purification and Sanger sequencing. We corroborated sequence quality by the quality scores provided by Macrogen and by visually examining the chromatograms using the software Ugene (Okonechnikov et al. 2012).
To include genome-wide markers in our dataset, we sequenced 2,180 Ultra Conserved Element (UCE) loci using a recently published bait set specific to bees, ants, and other apoid wasps ("hym-v2-bee-ant-specific"; Grab et al. 2019). This bait set is a subset of the principal Hymenoptera bait set first reported in Branstetter et al. (2017). For each sample, we sheared the DNA using a Qsonica Q800R2 acoustic sonicator, with the target fragment size range being 400-600 bp (60-120 secs shear time, 25% amplitude, 10-10 sec pulse). For older samples with more degraded DNA, we adjusted the shearing times to between 30-60 seconds. We cleaned fragmented DNA at 3× volume using a homemade SPRI-bead substitute ("speedbeads"; Rohland and Reich 2012). We generated Illumina sequencing libraries for each sample using Kapa Hyper Prep Kits (Roche Sequencing and Life Science, Wilmington, MA) and custom, dualindexing adapters (Glenn et al. 2019).
We amplified libraries for 12 cycles, purified them using speedbeads (Rohland and Reich 2012), and quantified them using a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA). To enrich UCE loci, we pooled 10 libraries at equimolar concentrations. Then, up to 500 ng of each pool was enriched following the manufacturer's protocol for day 1 (MYcroarray enrichment protocol v3.02) and the standard UCE protocol for day 2 (enrichment protocol v1.5 available at ultraconserved.org). The custom bait set was diluted 1:4 (1μL bait, 4μL H2O) with enrichment incubation at 65 °C for 24 hours using strip tubes and a PCR thermal cycler. For the second day of enrichment, we used 50uL of streptavidin beads per sample and performed on-bead PCR following the three heated (65 °C) wash steps. We amplified the enriched pools for 18 cycles and the resulting products were purified with SPRI beads at 1× volume. We sent the sequencing pools to the University of Utah Genomics Core for sequencing on an Illumina HiSeq 2500 (2×125, v4 chemistry).

Molecular approaches for species phylogeny and delimitation
For the molecular analyses, we generated DNA sequences from 19 males (less than 10 years old after collection) for the amplification of COI and UCEs. We preserved remaining body parts and associated DNA extractions for future studies (Suppl. material 1: For the mitochondrial data, we aligned COI sequences using the multiple sequence alignment online tool MAFFT (Katoh et al. 2019) and manually edited them using the software Ugene (Okonechnikov et al. 2012). We inferred a phylogenetic tree using Bayesian Inference (BI) method. Genetic distances within and between species were calculated in MEGA-X (Kumar et al. 2018) using 10,000 bootstraps. We performed Bayesian phylogenetic analyses in MrBayes v.3.2 (Ronquist et al. 2012) with COI sequences using the best nucleotide model estimated by PartitionFinder2 (Lanfear et al. 2017). The Markov Chain Monte Carlo (MCMC) was run for 20 million generations sampled every 1000 th generation. We discarded 25 percent of the first trees as burnin. We visualized and edited the Bayesian trees in FigTree v. 1.4.4 (Rambaut 2018).
For the UCE dataset, we performed most data processing steps using the software package Phyluce (Faircloth 2016). We cleaned the reads for adapter contamination and low-quality bases using Illumiprocessor (Faircloth 2013), which functions as a wrapper around the software Trimmomatic (Bolger et al. 2014). We assembled cleaned reads de novo for each individual using SPAdes v. 3.12.0 (Bankevich et al. 2012). To identify UCE regions from the bulk of assembled contigs and to remove paralogs, we used the Phyluce script match_contigs_to_pobes and the HymV2-bee-ant UCE bait files from Grab et al. (2019). We aligned all the loci individually using MAFFT (Katoh and Standley 2013) as implemented in Phyluce package, and trimmed resulting alignments using Gblocks with reduced stringency parameters (Castresana 2000 ). We removed loci that had data for fewer than 75% of taxa and generated a concatenated matrix from the resulting alignment set.
For phylogenetic reconstruction using UCE data, we inferred phylogenetic trees using Maximum Likelihood (ML), and multi-species coalescence reconstruction . We performed maximum likelihood analyses using two different strategies: single concatenated alignment and partitioned based on the best-fitting partitioning scheme. For the concatenated alignment, we obtained the substitution model (TVM+F+R2) using ModelFinder (Kalyaanamoorthy et al. 2017) which is part of the IQ-TREE v2.1.1 software (Nguyen et al. 2015). The best-fitting partitioning scheme was obtained using Sliding-Window Site Characteristics (SWSC), which divides each UCE into three data blocks corresponding to the right flank, core, and left flank (Tagliacollo and Lanfear 2018). We analyzed the resulting data subsets using PartitionFinder2 (Lanfear et al. 2017) using the rclusterf algorithm with AICc model selection criterion and GTR+G model of sequence evolution obtained by ModelFinder. We used the likelihood-based program IQ-Tree v2.1.1 for phylogenetic reconstruction of both partitioning schemes. To assess branch support, we performed 1000 replicates of the ultrafast bootstrap approximation (UFBoot; Hoang et al. 2018) and 1000 replicates of the branch-based Shimodaira-Hasegawa approximate likelihood-rate test (SH-aLRT; Guindon et al. 2010) using the command '-alrt'. Only clades with support values of UFBoot ≥ 0.95 and SH-aLRT ≥ 0.80 were considered robust. To account for heterogeneous gene histories that may influence phylogenetic accurate resolution, we inferred a species tree under the multi-species coalescent model using the program ASTRAL-III v.5.7.3 (Zhang et al. 2018), using gene trees provided by IQ-Tree as input. Support was assessed as local posterior probability, with ≥ 0.95 considered robust.
We also estimated a coalescent-based species tree using *Beast (Heled and Drummond 2010) and the BEAST2 package (Bouckert et al. 2019). Due to computational constraints, we selected a subset of loci for the *BEAST analysis. We ran the command 'phyluce_align_get_informative_sites' from Phyluce to identify the 100 most informative genes. From these most informative genes, we then selected those that were present in all samples, resulting in 88 UCEs. The analysis was run for 100 million generations sampling every 10,000 generations under a strict clock model with a constant population model, and a Yule model as a tree prior. We used a GTR model (unlinked across loci) for the nucleotide substitution model that was provided by PartitionFinder2 (Lanfear et al. 2017). To examine the convergence across the four runs performed and the ESS values of sampled parameters, we used Tracer v1.7 (Rambaut et al. 2018). A maximum clade credibility was constructed in TreeAnnotator and visualized the tree using Densitree, both included in the BEAST package (Bouckert et al. 2019).

Data availability
The raw UCE sequence reads have been uploaded to the NCBI Sequence Read Archive under BioProject accession PRJNA875942.

Morphometric geometrics
Individuals identified as E. cingulata and E. pseudocingulata formed one cluster based on the shape of the head and wings (Fig. 4). The PCA captured the variation of head and wing shape in 8 and 32 PCs, respectively. The first three components of the head variation explained 45.28%, 19.16%, and 12.37% of the covariance, respectively, total ing 76.81% (Suppl. material 1: table S3). The cross-validation test correctly located 62.07% and 66.67% of specimens of E. cingulata and E. pseudocingulata, respectively. The first three components of the wing shape explained 13.32%, 10.84%, and 10.02% of the variance, respectively, explaining a total of 34.18% of the variation (Suppl. material 1: table S4). Similar to the head results, the cross-validation test correctly located 66% and 68% of specimens of E. cingulata and E. pseudocingulata, respectively (Table 1). The multivariate regression analysis showed that, after 10,000 permutation rounds, the influence of the allometric effect was statically significant (P < 0.0001) with 32.28% and 3.07% predicted shape variation of the head and wing shape, respectively. Nonetheless, even after removing this allometric effect, the groups remained undifferentiated.

Genetic distance
The 655 bp fragment of the mitochondrial COI region resulted in an average pairwise genetic p-distance within E. cingulata of 1.3% and within E. pseudocingulata of 0.7%, while between the two species was 0.9%. The greater amount of genetic differentiation within the group of E. cingulata than between specimens from the two nominal species indicated no support for the presence of two evolutionarily independent units. The average pairwise genetic distances among E. cingulata and the outgroups (Eulaema mocsaryi, E. nigrita, and E. meriana) are greater than 5.6% (Table 2).

Phylogenetic relationships
Both the COI and UCE phylogenetic reconstructions provided no support for the species differentiation of E. cingulata and E. pseudocingulata (Figs 5, 6). For COI, PartitionFinder2 identified GTR+I as the best nucleotide substitution mod-  el. The consensus tree obtained from MrBayes showed that the COI fragments grouped E. cingulata and E. pseudocingulata into one highly supported clade sister to E. mocsaryi (Fig. 5). For the UCE data, we recovered a total of 21,470,135 reads    229 -1,831). The maximum likelihood inference of the UCEs markers recovered identical topologies from both partitioned and concatenated schemes, with most nodes showing high support (Fig. 6A, B). Individuals identified as E. cingulata and E. pseudocingulata form a monophyletic group with maximum support value. The difference between these topologies is limited to the phylogenetic position of E. pseudocingulata TA03 and TA04 within a small pseudocingulata clade. The only individual from the Atlantic Forest (E. cingulata TA09) is recovered with maximum support as sister to the individual from Mato Grosso (E. cingulata TA07), a state in the Amazon Forest. The ASTRAL and * Beast species tree showed maximum support values and both recovered the same topology (illustrated by the ASTRAL tree on Fig. 6C). Even though the COI and UCE data included different samples, both methods present very similar topologies with no significant differences between the individuals grouped in the two nominal species. Both trees showed one clade with individuals from the two species with high statistical support.

Discussion
Our results indicate that specimens identified as E. cingulata and E. pseudocingulata do not show morphological or genetic differentiation. Geometric morphometrics of the forewings has been previously used as a powerful technique to discriminate bee species (Francoy et al. 2012;Combey et al. 2013), subspecies (Oleksa and Tofilski 2014;Silva et al. 2015), cryptic species (Francisco et al. 2008;Hurtado-Burillo et al. 2016) and geographical ecotypes (Francoy et al. 2011;Grassi-Sella et al. 2018;Carneiro et al. 2019). Using both landmarks and outlined-based methodologies, Francoy et al. (2012) showed that the use of this approach to discriminate Euglossa species was more effective than studies using allozymes and restriction patterns of mitochondrial genes. Quezada-Euán et al. (2015) identified differences in wing shape of the two different morphs of Euglossa viridissima Friese, 1899 that had been only otherwise identified by the number of mandibular teeth. We applied the same methodology to this study to assess the taxonomic status of the two focal species in the genus Eulaema. Similarly, the results obtained with the head measurements were congruent with the results obtained with the forewings. Despite being used less frequently than the wings, the use of landmarks on the head has also been informative to recognize intercastes in honey bees (Souza et al. 2015), and to discriminate the morphologically indistinguishable females of the Psychodopygus complex (Diptera) (Godoy et al. 2018).
Because the level of morphological differentiation among recently divergent lineages is sometimes insufficient to recognize species, we used mitochondrial and UCE data to test the presence of genetically distinct groups of individuals among the specimens studied. DNA barcodes are increasingly becoming a standard tool used by taxonomists and its association with morphological characters has proven useful at discriminating species in several groups of bees (Gibbs 2009). Although this method has been criticized by some authors (Rubinoff et al. 2006;Wheeler 2008), it gives additional support to the recognition of species when considered along with other data sources (Padial and De La Riva 2007;Packer et al. 2009). Based on our data, the pair wise sequence divergences within and between the two species were below 3%, which according to Hebert et al. (2003) is a result compatible with the expected variation within a single species. More importantly, specimens of E. cingulata and E. pseudocingulata did not form a distinct monophyletic clades and the genetic distance among individuals of E. cingulata was greater than the genetic distance among individuals of E. cingulata and E. pseudocingulata. Dick et al. (2004) found that mtDNA divergences within Euglossini species were consistently low, with divergences among populations separated by the Andes averaging 1.1% (collection sites cover 3,000km). These findings and other studies have indicated the presence of high levels of long-distance gene flow between orchid bee populations. Rocha-Filho et al. (2013) observed a comparatively high dispersal ability in E. cingulata through genetic analyses comparing mainland and island populations. López-Uribe et al. (2014) also found low values of mitochondrial nucleotide divergence between populations of three species of Eulaema, showing a minimum value of 0.39% within species divergence in E. cingulata. According to the authors, the low sequence divergence between populations of E. cingulata was partially explained by the recent origin of this species.
The topology obtained with genomic data also supports the monophyly and recognition of one clade that comprises E. cingulata and E. pseudocingulata. UCEs have been successfully used as a tool for species discrimination as they provide sufficient variation at shallow time scales (Smith et al. 2013;Gueuning et al. 2020). Combined phylogenetic and population genetic approaches have been effectively used to investigate boundaries between complexes of wild European bees suspected to harbor cryptic diversity, mitochondrial introgression, or mitochondrial paraphyly (Gueuning et al. 2020). Using COI and UCEs with the multispecies coalescent method (BPP), Gueuning et al. (2020) also concluded that UCEs can provide robust species hypotheses and outperform COI in species delimitation. The adoption of delimitation methods based on the multispecies coalescent model has been criticized by Sukumaran and Knowles (2017), who argued that these methods tend to delimit population structure instead of species. The authors' concern was raised by the possibility of taxonomic inflation if species are described based only on molecular data. However, most studies using genetic data in species delimitation also incorporate addiotal sources of data such as morphology or morphometrics.
Herein, we conclude that E. pseudocingulata is not an independent evolutionary lineage from E. cingulata suggesting that the morphological differences observed in the velvety area of the mid-leg of males are a species polymorphism in E. cingulata. The difference in the shape of the mid tibia velvety area, proposed as the diagnostic character between the species, can be interpreted as a variable condition: it can be narrower and farther from the rear edge in some morphs occurring in the Amazon forest, or wider and closer to the rear edge in morphs occurring throughout the species distribution. Initially, the color of the abdomen was also proposed as a diagnostic character to differentiate E. cingulata and E. pseudocingulata, but a gradient can be observed in the two morphs, varying from a yellowish tone to orange (Fig. 7A-D). Additionally, color can be a highly variable trait within bees making it a difficult character for taxonomic identification. Variable color patterns on the abdomen have been described for several species of bumble bees (e.g., Carolan et al. 2012;Huang et al. 2015 ). Color variation in Eucerini bees was reported by Grando et al. (2018), in which they found two distinct color patterns in sympatric populations of Melissodes nigroaenea (Smith, 1854) in Brazil. Variation in color and shape was observed in Augochlora amphitrite (Schtottky, 1909) by Lepeco and Gonçalves (2018) using morphometric analyses and studying the male genital capsules. The authors did not find any character that support the recognition of distinct color morphs and macrocephalic females as different species. Similar to our findings, a genetic study of different color morphs of Euglossa species from the Atlantic Forest did not support the recognition of different species (Ferrari and Melo 2014).
The mechanisms responsible for maintaining the variation of the velvety area of the mid tibia in individuals of E. cingulata in the Amazon Forest remain unknown. However, a plausible explanation for the lack of differentiation in morphological and genetic markers is that there is an ongoing speciation process driven by sexual selection in the Amazonian population of E. cingulata. Such a rapid speciation process has been described in two sympatric Euglossa species from southern Mexico: E. viridissima and E. dilemma (Eltz et al. 2011). These sister species can be partially differentiated by the number of mandibular teeth: E. dilemma males possess three teeth on the mandibles while E. viridissima mostly show two teeth with some individuals expressing three teeth. However, these species are unequivocally distinguished by chemical characters (cuticular hydrocarbons found in the hind tibia) as well as by highly variable DNA markers (microsatellites and SNPs) (  A similar process could be occurring in E. cingulata and E. pseudocingulata but chemical information and characterization of genetic polymorphism across the entire genome of these bees would be necessary to properly investigate these questions. Further investigations are necessary to understand potential processes of incipient speciation in characters that have been characterized thus far.

Diagnosis and comments.
We recognized the taxon present in the Amazon forest under the name E. pseudocingulata as a junior synonym of E. cingulata, in light of the morphometric and phylogenetic data presented here. E. cingulata includes individuals with velvety areas of the midtibia that vary in width of the smooth area near the posterior edge (Fig. 1B, E). The basal tuft can be narrower or wider and slightly sloping. Other characters in this species include the labrum with lateral carina smoothly curved at the apex and ending distantly the edge of the clypeus, and the middle carina much shorter than lateral carina. The color of the metasomal hairs varies from orange to a light yellowish color ( Fig. 7A-D). The material examined is the same as referenced in the Suppl. material 1: tables S1, S2. Lectotype of Eulaema cingulata was examined by the images available in Nemésio (2009) and the holotype of E. pseudocingulata in Almeida et al. (2020).

Conclusion
We tested the species hypothesis of Eulaema cingulata and E. pseudocingulata by integrating multiple independent datasets: geometric morphometrics, phylogenetics using mitochondrial DNA, and phylogenomics using ultraconserved elements. All results across methods were congruent, showing no separation between morphs previously recognized as different species. Our results also suggest that the morphology of the mid tibia of E. pseudocingulata, proposed as the diagnostic character between the morphs, appears as a variable condition in some individuals of E. cingulata from the Amazon basin. Besides the variation in the mid leg, there is also color variation across samples. The evolutionary drivers of this variability are currently unknown. In the interests of nomenclatural stability, we have designated E. pseudocingulata as a junior synonym of E. cingulata. Orchid bees are important pollinators in Neotropical forests, being widely used in environmental quality studies, and they have become a good model for evolutionary genetics studies. However, as shown here, there is still a need to improve knowledge of the alpha taxonomy of these important group of Neotropical pollinators. Integration of DNA sequence analyses with geometric morphometrics of heads and wings can more rigorously test species boundaries than traditional morphological assessments alone and can ultimately improve species descriptions and identification tools.