Combining Morphology and DNA Barcoding Resolves the Taxonomy of Western Malagasy Liotrigona Moure, 1961 (Hymenoptera: Apidae: Meliponini)

ABSTRACT The taxonomy of the western Malagasy stingless bees (Liotrigona spp.) has been complicated by the high degree of morphological similarity between the described species. This group includes one of the smallest bee species of the world, Liotrigona bitika (Brooks & Michener, 1988). However, the description of this species is solely based on size differences, and observation of individuals intermediate in size between Liotrigona bitika and Liotrigona madecassa (Saussure, 1890) have cast doubt on the distinctness of the former taxon. I examined specimens, collected in the dry deciduous forest of Kirindy, representing all three of the described Liotrigona species of western Madagascar using a combined approach of size measures and DNA sequence analysis of the cytochrome oxidase I barcoding region. The analysis revealed all described taxa to be valid and furthermore uncovered the existence of a fourth taxon, described as Liotrigona kinzelbachi sp. n. The sympatric species exhibit non-overlapping body size ranges, indicating a strong influence of body size in niche differentiation of this otherwise morphological highly similar clade.


INTRODUCTION
Stingless bees (Meliponini) are eusocial bees with a pantropical distribution (Roubik 1989). Madagascar is home to seven described species in the genus Liotrigona (Moure 1961;Pauly et al. 2001). They are most abundant in the dry western part of the island, where three species are known, including one of the smallest bee in the world, Liotrigona bitika (Brooks & Michener 1988) with a body length of less than 2 mm. However, the de scription of these morphologically similar species is almost exclusively based on size differences. Because of an observed continuous variation in size between collected specimens of L. bitika and L. madecassa (Saussure 1890), the taxonomic distinctness of L. bitika has been questioned (Pauly et al. 2001). Furthermore, Brooks and Michener (1988) reported a group of specimens collected in the Tulear Province that is intermediate in size between L. madecassa and L. bitika, but did not describe them because of overlap with both species. To resolve the taxonomy of this group, specimens corresponding in size to all three described species were collected in the dry deciduous forest of Kirindy, western Madagascar. All specimens were measured for the metrics reported in Brooks & Michener (1988) and the cytochrome oxidase I barcoding region was sequenced for a selection from each morphological cluster. DNA barcoding has recently been applied to bees, and was shown to reliably distinguish between different species ). Integrating DNA barcoding with morphology has also previously helped in resolving taxonomically difficult groups of bees and other organisms (Gibbs 2009;Packer et al. 2009).

Collection
Bees attracted to sweat were collected at several localities in the dry deciduous forest of Kirindy, western Madagascar during November 2009. Bee specimens were stored in a DMSO/EDTA/NaCl buffer (Seutin et al. 1991). The specimens from this study were deposited in the collections of the Natal Museum, Pietermaritzburg, South Africa, and the Biosystematics Division of the Plant Protection Research Institute, Agricultural Research Council, Pretoria, South Africa (ARC-PPRI).

Measurements
Bee specimens (n = 76) were photographed under a Leica Wild M8 stereo light microscope together with a micrometer scale. The following measurements on the digital pic tures were taken using ImageJ (Rasband 1997(Rasband -2009: body length (in lateral view from frontal head margin to the tip of the abdomen), forewing length (from posterior tegular margin to wing tip), pterostigma length, scutal width (greatest width in front of the tegulae) and head width (greatest width in facial view). For body length, specimens were measured while immersed in ethanol without manipulation of the extension of the abdominal terga; furthermore, attempts were made to measure specimens in a straight position of the abdomen relative to the rest of the body, resembling a natural resting position.

Data analyses
For all morphological measurements, means and standard deviations were calculated. To identify discrete morphological groups, a principal component analysis (PCA) was carried out in JMP 7.0 (2007) on a covariance matrix of all measured variables. Prior to the PCA, all variables were log transformed to linearize the relationship between them and standardize the variances (Bookstein et al. 1985). The first two principal components were plotted.

Molecular protocols
DNA was extracted with a Wizard SV Genomic DNA extraction kit (Promega) from the whole thorax and abdomen of 34 of the 76 specimens used for the morphometric measurements, representing the range of variation observed for the four morphological clusters found in the PCA.
The cytochrome oxidase I (COI) barcoding region was PCR amplified using primers LepF (5'-ATT CAA CCA ATC ATA AAG ATA T-3') and LepR (5'-TAA ACT TCT GGA TGT CCA AAA A-3') ). For specimens showing poor PCR success, LioR (5'-CCA AAA AAT CAA AAT AAA TGT TG-3') was used as reverse primer. The 25 µl PCR-mix consisted of 0.5 µl of each primer (10 µM), 5 µl 5x PCR buffer, 0.125 µl of each dNTP (50 µM), 0.75U of Taq polymerase, 15.35 µl of ultrapure water and 3 µl of DNA template. The following cycling regime was applied: 1 min at 94 °C followed by 35 cycles of 1 min at 94 °C, 1.5 min at 51 °C, 1.5 min at 72 °C and a final extension of 5 min at 72 °C. Amplification products were checked on a 1.5 % agarose gel, and PCRs generating a single band were used for cycle sequencing. 10 µl of PCR product was digested with Exonuclease I (10U) and Shrimp Alkaline Phosphatase (1U) for 15 min at 37 °C prior to cycle sequencing to remove leftover primers and dNTPs. Products were sequenced from both ends with BigDye version 3.1 and primers LepF or LepR/LioR, and were analysed on an ABI 3730xl sequencer. Sequences were deposited in GenBank (accession numbers HQ012801-HQ012834).

Analyses of molecular data
Sequences were aligned with ClustalW (Thompson et al. 1994, http://align.genome. jp). An appropriate model of sequence evolution was determined in jModelTest 0.1 (Posada 2008) using the model with the lowest Akaike information criterion. A maximum likelihood phylogenetic tree was produced using a GTR+γ model in PhyML (Guindon & Gascuel 2003). Branch support was assessed with 500 bootstrap replicates. The tree was midpoint rooted because of the lack of a cytochrome oxidase I sequence of a suitable outgroup. A maximum parsimony tree was generated in TNT version 1.1 (Goloboff et al. 2008) using a search with TBR branch swapping and a random addition sequence with 100 replicates, the same settings were used for 500 bootstrap replicates. Furthermore, a neighbour joining tree was generated in MEGA version 4.0 (Tamura et al. 2007) with a Kimura 2-parameter (K2P) model and pairwise deletion of missing data as well as 500 bootstrap replicates. Mean sequence divergence within and between species was assessed using a K2P model in MEGA.

Morphology
The results of a principal component analysis on the morphological measurements are summarized in Figure 1. The collected specimens cluster in four distinct groups, separated by PCA axis 1. Detailed measurements for these four groups are listed in Table 1. The measurements for the individuals of three of the four clusters agree with those of previously described species of Liotrigona (Brooks & Michener 1988;   largest specimens correspond to L. mahafalya (Brooks & Michener 1988), the second largest to L. madecassa, and the smallest to L. bitika. The cluster of a fourth and novel species (described in this paper as L. kinzelbachi sp. n.) falls in between L. madecassa and L. bitika. While overlapping in body length and scutal width with both L. madecassa and L. bitika, this species can be distinguished by its distinct forewing length, pterostigma length and head width (Table 1).

Molecular Analysis
A midpoint-rooted maximum likelihood tree for all sequenced specimens is presented in Fig. 2 with branch support values from the maximum likelihood, neighbour joining and maximum parsimony bootstrap analyses. Individuals group into four well supported clades corresponding without exception to the four clusters observed in the PCA on the size measurements (Fig. 1).
The average genetic distance between species is more than one order of magnitude higher than within one species. Within-group distances average between 0 % (L. bitika) and 0.14 % (L. madecassa). Between-group distances range from 2.1 % (L. madecassa to L. mahafalya) to 3.3 % (L. kinzelbachi sp. n. to all three other species). Both L. bitika and L. kinzelbachi sp. n. are thus clearly distinguishable in terms of morphology and COI sequence from the other two species. Queen and male. Unknown. Comparison: Very similar to L. bitika and L. madecassa, but intermediate in size between the two. Distinguishable from both species by a combination of forewing length, pterostigma length, head width and scutal width (see Table 1). Also distinguishable from other species of Malagasy Liotrigona by its COI sequence (GenBank accession numbers HQ012826-HQ012834).

DISCUSSION
This study provides another example of the successful integration of morphology and DNA barcoding to resolve the species boundaries of taxonomically difficult taxa. Since the different species of Liotrigona are extremely similar in their morphology, it is difficult to define species in this group using morphological characters alone. Three of the four clusters of the morphometric analysis (Fig.1), however, correspond well with the measurements given for the three species described in Brooks and Michener (1988 ;  Table 1), and can therefore be identified as L. bitika, L. madecassa and L. mahafalya, even though type material could not be examined in this study. The fourth cluster appears to belong to a novel species, described as L. kinzelbachi sp. n. in this paper. Still, only L. mahafalya is clearly separated from the other species by distinct gaps in size measurements related to wing, head and thorax. The other species rather show a gradient in sizes from L. bitika over L. kinzelbachi sp. n. to L. madecassa. This has probably led to confusion over the validity of L. bitika as a distinct species (Pauly et al. 2001) and has caused L. kinzelbachi sp. n. to be overlooked in the past, even though the possible existence of an additional species similar to L. madecassa was already hinted at by Brooks and Michener (1988). However, adding COI sequence data as independent additional evidence shows these four clusters to be genetically distinct taxa. Further molecular and morphological studies should investigate the relationship of these species to the only relatively recently discovered Liotrigona species in the eastern rainforests of Madagascar (Michener 1990;Pauly et al. 2001). The high genetic similarity between the western Malagasy species indicates recent speciation events. Similarly, Rasmussen and Cameron (2010) date the split of L. ma decassa and L. mahafalya at less than 10 Mya. However, comparison with the two Lio trigona species on the African mainland (Eardley 2004) is needed to reveal if Madagascar was colonized multiple times or a radiation happened on the island.
The sympatric western Malagasy stingless bees are separated by the first PCA axis, which is generally seen to represent size (Jolicoeur & Mosimann 1960). This is also apparent by non-overlapping ranges of the measures for head width, forewing length and pterostigma length. The measured body length is overlapping between the different species, but differences in body posture and extension of the abdomen render this a less reliable character to distinguish size (Brooks & Michener 1988). Size has been implicated before in playing a role in niche differentiation in stingless bees. Whereas larger bee species can win in direct competition for nesting sites, smaller bee species are able to evade this competition by colonizing smaller nesting cavities (Hubbell & Johnson 1977;Michener 2001). While the nesting habits for only two of the western Malagasy stingless bees are known, this hypothesis is supported by the largest species (L. mahafalya) nesting in tree cavities and the smaller L. madecassa nesting in the smaller cavities of hollow bamboo stems (Brooks & Michener 1988). Alternatively, body size has also been shown to influence thermal properties in stingless bees, with smaller bees losing and gaining heat more quickly (Pereboom & Biesmeijer 2003). In an environment like the dry forests of Madagascar with large daily temperature fluctuations this could lead to different temporal activity patterns. Larger bees are able to retain heat more efficiently, and can thus forage earlier during the day, while smaller bees can more efficiently regulate their body temperature during hotter periods of the day (Biesmeijer 1997;Pereboom & Biesmeijer 2003). Both mechanisms could play a role for niche differentiation in the Malagasy stingless bees, but further research is needed to compare the daily activity patterns of the different species and to find the nesting sites of L. bitika and L. kinzelbachi sp. n. Due to the abundance and ease of collection, the western Malagasy Liotrigona could represent an ideal study subject for the role of body size in niche differentiation of closely related insect species.