Phylogenetic and Morphological Characteristics Reveal Cryptic Speciation in Stingless Bee, Tetragonula laeviceps s.l. Smith 1857 (Hymenoptera; Meliponinae)

Simple Summary The status of DNA barcoding in the cryptic species of stingless bees from Borneo, Tetragonula laeviceps sensu lato (s.l.) (Smith 1857), is poorly known. The T. laeviceps s.l. samples used in this study, which contain worker bee individuals grouped according to morphological characteristics and morphometric variations, could potentially have a similar grouping of COI haplotypes, but this has not yet been investigated. In this study, we investigate whether individuals of T. laeviceps s.l. worker bees grouped according to the same or nearly the same morphological traits show similar COI haplotype cluster groupings. The specimens were first classified according to the most obvious morphological characteristics, i.e., hind tibia color, hind basitarsus color and body size and grouping was based on morphological characteristics important for distinguishing the four groups within T. laeviceps s.l. The most distinctive features of the morphological and morphometric characteristics measured by PCA and LDA biplots that distinguished Group 1 from the other groups were the blackish-brown antennae scape (ASC) and the black (TC). Group 2 had a yellowish-brown ASC and a dark TC, while Group 3 had a dark brown ASC; a black TC; and the largest TL, FWW, and FWL. As for phylogenetic relationships, 12 out of 36 haplotypes showed clear separation with good bootstrap values (97–100%). The remaining haplotypes did not show clear differentiation between subclades that belonged together, regardless of their morphology and morphometric characteristics. Thus, the combination of DNA barcoding for species identification and phylogenetic analysis, as well as traditional methods based on morphological grouping by body size and body color, can be reliably used to determine intraspecific variations, such as the possible occurrence of subspecies within T. laeviceps s.l. Abstract Tetragonula laeviceps sensu lato (s.l.) Smith 1857 has the most complicated nomenclatural history among the Tetragonula genera. The objective of this study was to investigate whether T. laeviceps s.l. individuals with worker bees are grouped in the same or nearly the same morphological characteristics and have similar COI haplotype cluster groups. A total of 147 worker bees of T. laeviceps s.l. were collected from six sampling sites in Sabah (RDC, Tuaran, Kota Marudu, Putatan, Kinarut and Faculty of Sustainable Agriculture (FSA)), but only 36 were selected for further studies. These specimens were first classified according to the most obvious morphological characteristics, i.e., hind tibia color, hind basitarsus color and body size. Group identification was based on morphological characteristics important for distinguishing the four groups within T. laeviceps s.l. The four groups of T. laeviceps s.l. had significantly different body trait measurements for the TL (total length), HW (head width), HL (head length), CEL (compound eye length), CEW (compound eye width), FWLT (forewing length, including tegula), FWW (forewing width), FWL (forewing length), ML (mesoscutum length), MW (mesoscutum width), SW (mesoscutellum width), SL (mesoscutellum length), HTL = (hind tibia length), HTW (hind tibia width), HBL (hind basitarsus length) and HBW (hind basitarsus width) (p < 0.001). Body color included HC (head color), CC (clypeus color), ASC (antennae scape color), CFPP (Clypeus and frons plumose pubescence), HTC (hind tibia color), BSC (basitarsus color), SP (leg setae pubescence), SP (Thorax mesoscutellum pubescence), SPL (thorax mesoscutellum pubescence length) and TC (thorax color) (p < 0.05). The most distinctive features of the morphological and morphometric characteristics measured by PCA and LDA biplot that distinguish Group 1 (TL6-1, TL6-2 and TL6-3) from the other groups were the yellowish-brown ASC and the dark brown TC. Group 2 (haplotypes TL2-1, TL2-2 and TL2-3 and TL4-1, TL4-2 and TL4-3) had a dark brown ASC and a black TC, while Group 3 (haplotypes TL11-1, TL11-2 and TL11-3) had a blackish-brown ASC, a black TC and the largest TL, FWW and FWL. As for phylogenetic relationships, 12 out of 36 haplotypes showed clear separation with good bootstrap values (97–100%). The rest of the haplotypes did not show clear differentiation between subclades that belonged together, regardless of their morphology and morphometric characteristics. This suggests that the combination of DNA barcoding for species identification and phylogenetic analysis, as well as traditional methods based on morphological grouping by body size and body color, can be reliably used to determine intraspecific variations within T. laeviceps s.l.


Introduction
Tetragonula laeviceps sensu lato (s.l.) (Smith, 1857) is widely distributed in the Indo-Malayan region [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. This species differs from other species of the Tetragonula genus by having a predominantly shiny black body, longer wing length, a black mesoscutum, a metasoma that is completely covered with yellowish setae [1] and a hind wing with five hamuli [8][9][10]. According to Rasmussen and Michener [11], Tetragonula laeviceps s.l., formerly known as Trigona laeviceps s.l. (Smith, 1857), has the most complicated nomenclatural history among the genera of Tetragonula. Previously, Schwarz [12] considered specimens of this species to be Tetragonula iridipennis (Smith 1854), but later studies showed that T. iridipennis specimens are restricted to India and Sri Lanka, while other specimens of T. iridipennis identified by Schwarz in other areas were renamed T. laeviceps s.l. by Over the last few decades, DNA barcoding [20] has emerged as an important genetic tool for the taxonomic identification of bee species [21][22][23][24][25] that were difficult to morphologically identify and/or did not have dichotomous keys [26,27]. To date, most sequences have been obtained from the Meliponini cytochrome oxidase I (COI) barcode fragment to support morphometric data and address specific taxonomic questions [18,[28][29][30][31][32][33][34][35][36][37] or conduct phylogenetic and phylogeographic research [38][39][40]. DNA barcoding with the COI gene of mitochondrial DNA (mtDNA) is a protein-coding gene that contains more differences than the ribosomal gene [21]; therefore, it can be used to distinguish closely related species in addition to the morphological characteristics and morphometric variations of insects. To date, the status of DNA barcoding for the cryptic species of stingless bees from Borneo, T. laeviceps s.l., is poorly known. There is a possibility that T. laeviceps s.l. samples in this study, which contain worker bee individuals grouped in a particular set of morphological characteristics and morphometric variations, have a similar grouping of COI haplotypes, but this has not yet been investigated.

Stingless Bee Samples
A total of 147 T. laeviceps s.l. workers were collected directly from the nest entrance of the colony at six sampling sites in Sabah, including RDC, Tuaran, Kota Marudu, Putatan, Kinarut and the Faculty of Sustainable Agriculture (FSA) at Universiti Malaysia Sabah. The stingless bees were kept in collection bottles, stored in 100% ethanol and taken to FSA Sandakan for molecular analysis. The voucher specimens were deposited in the insect collection in the FSA herbarium. However, only 36 T. laeviceps s.l. individuals were selected for further systematic studies. These specimens were first classified according to the most obvious morphological characteristics, namely body size, hind tibia color and hind basitarsus color, which can be detected by eye as an important morphological feature, to create four distinct groups within T. laeviceps s.l. (Table 1). Additional T. laeviceps s.l. individuals with similar characteristics were excluded from this study. ences than the ribosomal gene [21]; therefore, it can be used to distinguish closely related species in addition to the morphological characteristics and morphometric variations of insects. To date, the status of DNA barcoding for the cryptic species of stingless bees from Borneo, T. laeviceps s.l., is poorly known. There is a possibility that T. laeviceps s.l. samples in this study, which contain worker bee individuals grouped in a particular set of morphological characteristics and morphometric variations, have a similar grouping of COI haplotypes, but this has not yet been investigated.

Stingless Bee Samples
A total of 147 T. laeviceps s.l. workers were collected directly from the nest entrance of the colony at six sampling sites in Sabah, including RDC, Tuaran, Kota Marudu, Putatan, Kinarut and the Faculty of Sustainable Agriculture (FSA) at Universiti Malaysia Sabah. The stingless bees were kept in collection bottles, stored in 100% ethanol and taken to FSA Sandakan for molecular analysis. The voucher specimens were deposited in the insect collection in the FSA herbarium. However, only 36 T. laeviceps s.l. individuals were selected for further systematic studies. These specimens were first classified according to the most obvious morphological characteristics, namely body size, hind tibia color and hind basitarsus color, which can be detected by eye as an important morphological feature, to create four distinct groups within T. laeviceps s.l. (Table 1). Additional T. laeviceps s.l. individuals with similar characteristics were excluded from this study.
To ensure that the species obtained in this study most closely resembled the morphological characters of T. laeviceps s.l., we also collected four species of the genus Tetragonula as comparative species, namely Tetragonula melanocephala s.l. (i. GenBank Accession Nos. ii. ON458757 (TL = 3.57 mm, SW = 1.05 mm)), but they were excluded from the present analysis. Tetragonula fuscobalteata s.l., which was used as an outgroup species for the phylogenetic analysis, was described by the smallest size, six distinct longitudinal pubescence and five glabrous areas on the mesoscutum and reduced submarginal of the forewing cells (see Figure A1) [4].  In this study, the overall measurements of the morphological characteristics of T. laeviceps s.l. were carried out as follows:

Body Color and Pubescence
Eleven (11) qualitative measurements of body pubescence and color of head, antennae and thorax were used in this study: head color (HC), clypeus color (CC), antennae scape color (ASC), clypeus and frons plumose pubescence (CFPP), hind tibia color (HTC), hind basitarsus color (BSC), hind leg setae pubescence (LSP), hind leg setae pubescence length (LSPL), thorax mesoscutellum pubescence (SP), thorax mesoscutellum pubescence length (SPL) and thorax color (TC) ( Table 2). The measurements of the body color and pubescence of T. laeviceps s.l. were based on the quantitative color code described by Abu In this study, the overall measurements of the morphological characteristics of T. laeviceps s.l. were carried out as follows:

Morphological Characteristics and Morphometric Parameters of T. laeviceps s.l.
In this study, the overall measurements of the morphological characteristics of T. laeviceps s.l. were carried out as follows:

Body Color and Pubescence
Eleven (11) qualitative measurements of body pubescence and color of head, antennae and thorax were used in this study: head color (HC), clypeus color (CC), antennae scape color (ASC), clypeus and frons plumose pubescence (CFPP), hind tibia color (HTC), hind basitarsus color (BSC), hind leg setae pubescence (LSP), hind leg setae pubescence length (LSPL), thorax mesoscutellum pubescence (SP), thorax mesoscutellum pubescence length (SPL) and thorax color (TC) ( Table 2). The measurements of the body color and pubescence of T. laeviceps s.l. were based on the quantitative color code described by Abu To ensure that the species obtained in this study most closely resembled the morphological characters of T. laeviceps s.l., we also collected four species of the genus Tetragonula as comparative species, namely Tetragonula melanocephala s.l. ii. ON458757 (TL = 3.57 mm, SW = 1.05 mm)), but they were excluded from the present analysis. Tetragonula fuscobalteata s.l., which was used as an outgroup species for the phylogenetic analysis, was described by the smallest size, six distinct longitudinal pubescence and five glabrous areas on the mesoscutum and reduced submarginal of the forewing cells (see Figure A1) [4].

Morphological Characteristics and Morphometric Parameters of T. laeviceps s.l.
In this study, the overall measurements of the morphological characteristics of T. laeviceps s.l. were carried out as follows:

Body color and pubescence
Head color (HC) 18.

DNA Extraction and PCR Amplification
To detect a possible differentiation of the mtDNA haplotypes, a total of 37 stingless bee individuals, comprising 36 T. laeviceps s.l. workers and 1 worker of Tetragonula fuscobalteata (outgroup species), were examined for molecular analyses ( Table 1). As indicated by Koch [33], DNA was extracted from the entire body of the stingless bee using the Wizard Genomic DNA Extraction Kit (Promega, Madison, WI, USA), as the thorax alone did not yield enough DNA for a subsequent analysis. The entire body of the stingless bee was determined to be the best source of mitochondrial DNA [41]. The targeted cytochrome c oxidase I (COI) primers LepF1 (forward) 5 -ATTCAACCAATCATAAAGATATTGG-3 and LepR1 (reverse) 5 TAAACTTCTG GATGTCCAAAAAATCA-3 were used for PCR amplification and DNA sequencing [33].

Sequence Editing and Alignment
All analyses were performed on the 680-nucleotide portion of the alignment of 37 sequences. The quality of the reads was checked using Chromas 2.6.6 (Technelysium DNA Sequencing Software). Raw sequence reads were processed using Mega X version 10.0.5 [42] and aligned using Clustal X [43]. Additionally, Mega X version 10.0.5 was used to assemble and edit the contigs (the consensus sequence from the region of overlap between the 5 and 3 sequences). The sequences were then imported into DnaSP v6 [44] to generate haplotypes. Clustal X and MEGA X version 10.0.5 were utilized to align the sequences with the T. fuscobalteata outgroup.

Correlations, PCA and LDA Biplot Analysis
The PCA biplot analysis used variables of the correlation matrix of 21 morphological and morphometric traits and the calculation matrix of the four T. laeviceps s.l. groups. All correlations between the main axes, the vectors of the variables (loadings) and the T. laeviceps s.l. groups (scores) are presented in the PCA biplot. To obtain a correct predictive classification between the groups of T. laeviceps s.l., the sample in the LDA was split from a dataset of the variables identifying a validation set of 3 groups of traits based on the selection algorithm: (i) A (TL, HW, HL, CEL, CEW, FWLT, FWW, FWL, ML, MW, SW, SL, HTL, HTW, HBL and HBW); (ii) B (CC, ASC and TC) and (iii) C (ML and CFPP). In this study, LDA used a PCA-based multivariate evaluation where the resulting canonical axes (Axis 1 and Axis 2) showed the largest separation between all groups.

Sequence and Phylogenetic Analysis
The maximum likelihood method (ML) was used to reconstruct the phylogenetic relationships between the T. laeviceps s.l. haplotypes. For phylogenetic reconstruction, a ML optimality criterion (as implemented in Paup 4.0a169 [45]) was used with the Kimura-2 parameter (K2P) model and T. fuscobalteata as the outgroup. An optimal model of nucleotide evolution of the ML analysis was determined using ModelTest 2 [46]. Branch swapping was performed through Tree Bisection Reconnection (TBR). An estimate of support for each branch was determined with the bootstrap test using 1000 pseudoreplications.

Identification of T. laeviceps s.l.
In Figure 1, the identification of T. laeviceps s.l. workers were confirmed; the main features include four longitudinal hairs that are not well banded, three glabrous areas of the mesoscutum, reduced submarginal cells of the forewing and the elliptical disc of the basitarsus, as described by Abu Hassan [4], Purwanto et al. [14], Suriawanto et al. [47], Siti Fatimah et al. [13] and Sakagami [7]. In the morphometric analyses, four groups with different morphological characteristics were found among the 36 T. laeviceps s.l. specimens,

Morphological and Morphometric Parameters
Analysis of T. laeviceps s.l.

Qualitative Analysis of Body Color and Pubescence
The Kruskal-Wallis test for body color and pubescence between the four groups revealed that the HC, CC, ASC, CFPP, HTC, BSC and TC were significantly different at p < 0.001, while LSP, SP and SPL were significantly different at p < 0.05. However, a nonsignificant difference was detected for the LSPL (Table 4; p > 0.05).

PCA and LDA Biplot Analysis
The PCA biplot in Figure 2 displays the correlation between the principal axes (PC1 and PC2), 20 morphological and morphometric parameters (loadings) and T. laeviceps s.l. groups (PCA scores). The first two principal components (PC1 and PC2) describe 91.56% of the total variance and have eigenvalues greater than the unit. Small angles indicate a strong correlation between the morphological and morphometric parameters. Thus, the group of parameters based on PCA are as follows: (i) A (body sizes; TL = total length, HW = head width, HL = head length, CEL = compound eye length, FWLT = forewing length (+tegula), FWW = forewing width, FWL = forewing length, MW = mesoscutum width, SW = mesoscutellum width, SL = mesoscutellum length, HTL = hind tibia length, HTW = hind tibia width, HBW = hind basitarsus width, and HBL = hind basitarsus length), (ii) B (compound eye width (CEW) and body color (CC = clypeus color, ASC = antennae scape color and TC = thorax color) and (iii) C (mesoscutum length (ML) and clypeus and frons plumose pubescence (CFPP)). The linear discriminant analysis (LDA) determines the group means and calculates the probability of each specimen belonging to different groups ( Figure 3). Based on the PCA multivariate analysis, the LDA also showed a clear separation of the four groups between canonical axes 1 and 2, and Groups 3 and 4 are also more closely related. The PCA biplot in Figure 2 indicates the correlation between the principal axes (PC1 and PC2), 21 morphological and morphometric parameters (loadings) and the T. laeviceps s.l. groups (PCA scores). The first two principal components (PC1 and PC2) describe 91.14% of the total variance and have eigenvalues greater than the unit. Small angles indicate a strong correlation between the morphological and morphometric parameters. In the PCA biplot, the four groups of T. laeviceps s.l. specimens are represented by convex hulls (Groups 1, 2, 3 and 4) that are completely separated from each other (Figure 2). Compared to the other groups, Groups 3 and 4 are the most closely related groups with the least differences in morphological and morphometric parameters in terms of body size, body color and pubescence; in these two groups, most endpoints of the respective parameters coincide. Groups 3 and 4 also have a larger body size than groups 1 and 2, and both groups have dark brown ASC and black TC. However, Groups 1 and 2 have the most distinctive features of the overall group parameters, with Group 1 having the most distinctive features of large CEW, as well as blackish brown ASC and black TC, and Group 2 having the most distinctive features of relatively larger sizes of FWLT and FWL although it is the third smallest after groups 3 and 4, as well as yellowish-brown ASC and dark brown TC.
The linear discriminant analysis (LDA) determines the group means and calculates the probability of each specimen belonging to different groups (Figure 3). Based on the PCA multivariate analysis, the LDA also showed a clear separation of the four groups between canonical axes 1 and 2, and Groups 3 and 4 are also more closely related. The linear discriminant analysis (LDA) determines the group means and calculat the probability of each specimen belonging to different groups (Figure 3). Based on t PCA multivariate analysis, the LDA also showed a clear separation of the four grou between canonical axes 1 and 2, and Groups 3 and 4 are also more closely related.

Group, Haplotypes
Morphological and Morphometrical Characteristics In this study, the indeterminacy of morphological differences within workers of T. laeviceps s.l. was observed, with differences in body size (Table 3; 16 traits, p < 0.001) and in body color and pubescence ( Table 4; 10 out of 11 traits, p < 0.05). The morphological traits of the four groups of T. laeviceps s.l. workers differed in size (see four groups of T. laeviceps s.l. in Table 1). A similar finding was also reported by Purwanto and Trianto [15] in Indonesia and Patel and Pastagia [48] in India. Table 3  Previous studies by May-Itzá et al. [18], Novita et al. [49], Quezada-Euán et al. [50] and Raffiudin et al. [51] found that the availability of food and nesting resources in the environment affects the size of the morphological features in stingless bees, such as body segments. As shown in Table 3, the largest mean morphological traits of T. laeviceps s.l. worker bees in relation to size for Group 3 are from the agricultural area where many oil palms were planted, followed by Group 4 from the natural and secondary forest, Group 2 from the urban vegetation and Group 1 from the secondary forest and teak plantation (Table 1). Changes in temperature or environmental conditions cause organisms to morphologically adapt to the environment, including flight and foraging activities [51]. This was further confirmed by Roubik [52] and Ruttner [53], who found that body size variations in worker bees are a form of adaptation in the foraging and exploitation of floral resources. Foraging distances also depend on environmental conditions, such as the density and distribution of food resources and the general physical carrying capacity of different habitats [54,55]. Larger bees with bigger thoraces tend to have greater flight distances, allowing them to access more extensive areas for food and nesting material sources from the environment [51][52][53][54][55][56]. Therefore, variations in some traits, such as the length of wings and legs, may be an important factor in the morphometry of stingless bees and can be used to estimate size differences between worker bees [57,58]. Furthermore, Sun and Chaplin-Kramer [55] found that the habitat suitability for pollinators is altered by changes in forest and pasture covers, which affect the morphological characteristics of stingless bees.
For the traits of body color and pubescence in T. laeviceps s.l. workers, 10 out of 11 morphometric traits differed significantly between the groups (Table 4; p < 0.05). In this study, the four groups can be clearly distinguished from each other on the basis of the following characteristics: head color (HC), clypeus color (CC), antennal color (ASC), clypeus and frons plumose pubescence (CFPP), hind tibia color (HTC), basitarsus color (BSC), leg setae pubescence (LSP), mesoscutellum pubescence (SP), mesoscutellum pubescence length (SPL) and thorax color (TC) (see four groups of T. laeviceps s.l. in Tables 1 and 4). Group 1 had a black HC, dark brown CC, blackish-brown ASC, black HT and BSC, dense LSP and SP, long LSPL and SPL and black TC. Group 2 can be distinguished by a blackish-brown HC, yellowish-brown CC and ASC, dark brown HT and BSC, sparse LSP and SP, short LSPL and SPL and dark brown TC. Group 3 showed a black HC, dark brown CC and ASC, blackish-brown HT and BSC, sparse LSP and SP, short LSPL and SPL and black TC. Group 4 had a black HC, dark brown CC and ASC, blackish-brown HT and BSC, sparse LSP and SP, short LSPL and SPL and black TC. Table 4 shows that the most obvious criterion for distinguishing the four groups by eye is the color of the hind legs (HTC) and the basitarsus (BSC). The worker bees have corbiculates, a specialized structure on their hind legs that is used to collect pollen and resin [59]. Therefore, the variations in the morphological characteristics of the hind legs of stingless bees, such as the color and shape of the corbiculae in this study, could be due to their coevolution with angiosperms and gymnosperms as pollinators [60,61]. Bomtorin et al. [61] found that the development of hind leg diphenism, characteristic of corbiculae bee species, is driven by several castepreferentially expressed genes, such as those encoding cuticular protein genes, P450 and Hox proteins, in response to the naturally diverse diet offered to honeybees during the larval period. As for the other characteristics of body color and hair of the T. laeviceps s.l. worker group, it was difficult to distinguish between them. Even though two of the three groups generally had the same color and hair patterns, these groups differed statistically, and therefore, further research is needed.

PCA and LDA Biplot and Phylogenetics Relationships of T. laeviceps s.l.
The principal component analysis (PCA) is a commonly used analytical technique in taxonomic research, because it can identify the role of each trait in each group formed [62,63]. In this study, both PCA and LDA biplots showed a complete separation of morphological and morphometric traits (variable Group A: HW, CEL, FLT, FW, FWL, MW, SW, HTL, HTW, HBL, HBW and ASC; variable Group B: TL, HL, CEW, SL, SPL and TC and variable Group C: CFPH) between groups 1, 2, 3 and 4 of T. laeviceps s.l. (Figures 2 and 3). However, Groups 3 and 4 had the most distinct group parameters, followed by Groups 1 and 2, indicating that the cryptic Bornean species T. laeviceps s.l. has morphological sizes and morphometric characteristics that are difficult to distinguish by eye and vary between habitats. This species occurs in a wide range of habitats in the Malaysian state of Sabah, from urban areas to forests and agricultural fields, and it can generally disperse over short distances. The results of this study suggest that the geological structure and diversity of tropical vegetation in Sabah, which may have been an important barrier to dispersal between populations (Table 1), is the most likely explanation for the observed patterns of variation in the morphological and morphometric characteristics of T. laeviceps s.l.
Of the 36 haplotypes, only 12 haplotypes from Group 3 (TL11-1, TL11-2 and TL11-3);  Group 2 (TL2-1, TL2-2 and TL2-3 and TL4-1, TL4-2 and TL4-3) and Group 1 (TL6-1, TL6-2 and TL6-3) of T. laeviceps s.l. showed clear separation with good bootstrap values (97-100%) (Figure 4). The other haplotypes did not show clear differentiation between the subclades, which belonged together regardless of their morphology and morphometric characteristics. In this study, the most important characteristics to distinguish between the groups of T. laeviceps s.l. were the color of the antennal scape (ASC) and thorax (TC), forewing length (FWL) and forewing width (FWW) ( Table 5). The result is further supported by the PCA and LDA biplots, which show that Groups 3 and 4 also have larger body size than groups 1 and 2, and both groups have dark brown ASC and black TC. Group 1 has a relatively large CEW, although it has the smallest body size, as well as blackish brown ASC and black TC, and Group 2 has a relatively larger FWLT and FWL, although it is the third smallest after groups 3 and 4, as well as yellowish-brown ASC and dark brown TC. Invertebrate species with low dispersal ability, such as bees and wasps, tend to be more affected by physical barriers that act as obstacles to movement between habitats, thereby affecting the community composition, leading to low or limited genetic exchanges [64,65]. Genetic differentiation can occur during evolution even if the morphological traits do not change [66,67]. However, certain external morphological traits may be retained due to selection pressure on important adaptive traits [68]. This study found that cryptic and polymorphic characteristics within species can complicate traditional taxonomic research. Both issues are common in communities of stingless bees such as T. laeviceps s.l. and deserve more attention during systematic identification.
The color of the hind tibia (HT), hind basitarsus (BSC), antennal scape (ASC) and the thorax (TC) are the most significant morphological features that distinguish the three groups of T. laeviceps s.l. (Tables 1 and 5, Figure 4). According to Castanheira and Contel [69], the color of the mesepisternum distinguishes the two subspecies of Tetragonisca angustula angustula and Tetragonisca angustula fiebrigi. In some cases, interbreeding between the species could also result in hybrids with different mesepisterna colors [70][71][72]. This could be related to differences in the spatial regulation of the expression of other genes during development, such as the Hox genes, which, in turn, lead to morphological differences within and between taxa [73]. Jeong et al. [74] found that the decline in male abdominal pigmentation in the lineage of Drosophila santomea was related to the selective loss in expression of both light brown and yellow pigmentation genes in the posterior abdominal segments of males. Previous studies have also demonstrated that changes in the expression and function of Hox genes are key to the evolution of novel body structures in insects, such as segmental differences within the body and thoracic legs in some Orthoptera and Hemiptera species, and possibly abdominal appendages in Ephemeroptera [75,76]. Averof [76] found that Hox gene mutations can convert characteristic structures of one segment into the corresponding structures of another segment. Furthermore, the gene responsible for bee coloration in this study could be determined by the F1 generation, which is expressed in both queens and worker bees and confers a similar phenotype [77][78][79].
Another important feature of T. laeviceps s.l. identified by PCA and biplot analyses is the size of the FWL and FWW (Table 3, Figures 2-4). In this study, the bees with the largest mean total body length also had higher mean values for other morphological characteristics by body size. Group 3 is the largest of the T. laeviceps s.l. groups, followed by Group 2 and Group 1 ( Table 5). Group 2 is characterized by relatively large FWLT and FWL values, although it is the third smallest after Groups 3 and 4 (Table 3, Figure 2). Rasmussen and Michener [11] also pointed out that the earliest studies of T. laeviceps s.l. by Sakagami [7] also showed that different body sizes can be considered as different species, but further systematics studies are needed. In Guatemala and Mexico, morphometric analyses based on body size in combination with molecular tools (DNA barcoding) have been used to identify complex species of the stingless bee Melipona yucatanica [34,36]. However, Figure 4 does not show a clear distinction between the subgroups of some haplotypes of T. laeviceps s.l., which belong together regardless of their morphology and morphometric characters. Based on this study, the morphologically cryptic species T. laeviceps s.l. needs to be further characterized, as the insects found on the islands may have different speciation processes [80].
In addition, morphometric measurements of the hymenopteran forewing (length and width of the wing, length of marginal and basal veins) can be used in the high-throughput protocol [81], as is the case with DNA barcoding. The high-throughput is an advantage over traditional identification methods, which are slow and require a high level of expertise [32]. The ability of DNA barcoding to discriminate between different species is reflected in the low intraspecific variation and relatively high interspecific variation [33]. In this study, the COI sequence variation analysis, together with the traditional methods of a PCA and biplot analysis, was useful in identifying the morphological cryptic complexes within T. laeviceps s.l. In addition, further analyses such as the multivariate ratio analysis (MRA) based on the insect body ratio [81,82] and geometric morphometrics (GM) based on the insect body outline [83] could be conducted in the future to confirm the results of the present study.

Conclusions
This study concludes that the combination of DNA barcoding for species identification and phylogenetic analysis, as well as traditional methods based on morphological grouping by body size and body coloration, can be reliably used to evaluate intraspecific variations such as the possible occurrence of subspecies within T. laeviceps s.l. The most significant morphological characteristics of T. laeviceps s.l. for distinguishing between groups or possibly subspecies were body size (FWL and FWW) and body coloration (TC and ASC). Furthermore, there were relatively high bootstrap values and genetic distances within the phylogenetic relationships of T. laeviceps s.l. These led to the separation of three distinct haplotype groups (12 out of 36 haplotypes) and suggest the existence of intraspecific cryptic speciation that warrants further investigation.

Institutional Review Board Statement:
The authors declare that all research complied with ethical standards in accordance with the "Universiti Malaysia Sabah: researcher's guidelines on code of practise for the care and use of animals for scientific purposes" (https://www.ums.edu.my/v5 /files/2017/Guidelines-for-UMS-Animal-Ethics.pdf, accessed on 11 February 2023). In this study, no vertebrates were used as test animals and the test material consisted only of lower order animals (i.e., insect, including the trades colonies of stingless bees) that permit testing.

Informed Consent Statement:
Not applicable for studies not involving humans.

Data Availability Statement:
Data is contained within the article.

Acknowledgments:
The authors wish to thank the Sabah Forestry Department and Forest Department Sarawak for the fieldwork permission, the Faculty of Sustainable Agriculture, and the Centre for Research and Innovation, Universiti Malaysia Sabah. The authors also wish to thank the editors and reviewers.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of the data; in the writing of the manuscript or in the decision to publish the results. Figure A1. The identification of the worker Tetragonula fuscobalteata s.l. (a-c) described according to Abu Hassan [4] and Siti Fatimah et al. [13], which has the smallest size within the genus Tetragonula and also has a brown gastral tergite, a black mesoscutum with pale yellow pubescence, antenna color below the testaceus to ferruginous, submarginal cells with reduced forewing venation (d) and six distinct longitudinal pubescence (e).