A molecular phylogeny of the long-horned bees in the genus Melissodes Latreille (Hymenoptera: Apidae: Eucerinae)

The genus Melissodes Latreille in the subfamily Eucerinae (Hymenoptera: Apidae) is a widespread and common group of bees. There are 129 described Melissodes species that range throughout the western hemisphere with the center of diversity in the warm deserts of southwestern North America. Despite its widespread nature and importance in agriculture, the evolutionary relationships among the species have never been investigated. Here, we present a molecular phylogeny using five loci for 89 species of Melissodes with representatives from all subgenera. We confirm all of the subgeneric delineations constructed by LaBerge, except for a paraphyletic M. ( Tachymelissodes ) and the placement of M. ( Heliomelissodes ), which renders M. ( Eumelissodes ) paraphyletic. We also discuss the unexpected placement of M. tristis Cockerell , M. paucipuncta LaBerge, M. dagosus Cockerell, and M. pexa LaBerge. Finally, we combine this analysis with previous data to support the placement of Melissodes within the subfamily Eucerinae.


Introduction
The genus Melissodes Latreille (Apidae: Eucerinae sensu Bossert et al. 2019) is a diverse genus of medium sized, setaceous bees. It is the second largest genus in the tribe Eucerini with 129 described species ranging from Canada to Argentina. Melissodes bees are either solitary or gregarious ground-nesters with most species emerging in mid-to late summer. They are widespread and important pollinators in both natural and agricultural settings and are mentioned as prominent pollinators of sunflower (Parker et al. 1981), canola (Morandin & Winston 2005), cantaloupe (Winfree et al. 2007), watermelon (Kremen et al. 2002), alfalfa (LaBerge 1956a), cotton (LaBerge 1956a), coffee (Jha & Vandermeer 2010), and anecdotally on many other crops. Many of the species are polylectic (able to pollinate plants in two or more families) while others In Part I, LaBerge (1956a)  While Melissodes ranges throughout North and South America, each subgenus has a more restricted distribution (Fig. 3). Melissodes (Ecplectica) has only ten described species with five in South America, three in the Antillean islands, one species ranging from Panama north to Mexico, and no locality information found for the tenth species. Melissodes s.s. has 24 species mostly from southwestern North America (including Mexico), with only five species that occur further north or east, and six species that occur in the Antillean islands (one of which also occurs on the mainland of South America crossing into northern Brazil). Melissodes (Tachymelissodes) has only four species. Three are in southwestern North America and one is more widespread throughout the western United States, north to Washington. Melissodes (Callimelissodes), with 14 species, has its center of diversity in the western United States with four species extending east of the Mississippi and two of these also extending south into Mexico. Melissodes (Apomelissodes) has four species that occur in the eastern United States, with two as far west as Texas. Melissodes (Heliomelissodes) consists of two species, both with large distributions in North America, one primarily east of the Mississippi and one west of it. Melissodes (Psilomelissodes) is a monotypic subgenus that is restricted to Texas, Kansas and Nebraska. Finally, the largest subgenus, M. (Eumelissodes), has 72 described species ranging from British Columbia to Maine and south to Cuba and Panama. Most of the diversity is in the western United States and Mexico with about one third of the species ranging further east, one fifth found further north into Canada, and only two species in Central America, and one in Cuba (Ascher & Pickering 2017;Michener 2007).
For the tribe Eucerini, Michener (2007) listed 33 genera with the acknowledgement that much work was needed to better circumscribe them. Two recent studies investigated the relationships among many of the Eucerinae genera (Cardinal & Danforth 2013;Praz & Packer 2014), including Melissodes. And most recently, the Eucera Scopoli complex has undergone a major reclassification based on molecular and morphological data, in which Dorchin et al. (2018a) have reassigned six previously recognized genera to subgeneric status in the genus Eucera and erected one new genus (Dorchin et al. 2018b). This leaves the tribe Eucerini with 28 genera, Eucera being the largest and only cosmopolitan genus. There are three genera that occur only in the eastern hemisphere and the rest are in the western hemisphere. Of the western genera, six span North and South America, fifteen are Neotropical, and only two genera are solely Nearctic. Even though relatively few genera occur in North America as compared to South America, the two largest genera, Eucera and Melissodes, contain most of the species and are primarily Nearctic.
The two primary objectives of this study are (i) to present the first comprehensive molecular phylogeny of the genus Melissodes and test monophyly of the genus and subgenera and (ii) to understand the evolutionary relationships among these groups and contrast them with the groupings proposed by LaBerge (1956aLaBerge ( , 1961. This work will provide a framework for the future revision of the genus and a reference for studying the evolution of host plant specialization within this genus. Our outgroup sampling included one species of Exomalopsis Spinola and one species of Anthophorula Cockerell in the tribe Exomalopsini. Within Eucerini, we included thirteen species of Eucera in five subgenera per Dorchin et al. (2018a); two Xenoglossa Smith, three Xenoglossodes Ashmead, one Peponapis (Say), one Syntrichalonia LaBerge, and six Synhalonia Patton, as well as one species each of Martinapis Cockerell and Florilegus (Cresson), and five species of the sister genus Svastra.

Taxon sampling
The  Supplementary Table S1. All identifications were made or confirmed by the first author using LaBerge's keys (LaBerge 1956a(LaBerge , b, 1961 and reference material from the USDA-ARS Pollinating Insects Research Unit in Logan, UT and the American Museum of Natural History in New York, NY. In addition, DNA sequences from 11 taxa that overlap the breadth of this analysis were acquired from GenBank. These data came from Cardinal & Danforth (2013) and included a complementary taxon from each of the outgroup taxa in Exomalopsini, as well as five taxa in Tapinotaspidini, Emphorini and Ancylaini, three overlapping taxa within Eucerini, and a single sequence from within Melissodes. These additional taxa along with the large number of outgroup taxa included by the authors were used to more strongly place Melissodes within the Eucerinae.
We included some specimens that could not be confidently identified. To confirm that the unknown specimens were molecularly distinct, we ran a preliminary phylogenetic analysis in a maximum likelihood framework. We compared the position of the unknown specimens to their sister species on the resulting phylogeny. For the pairs with extremely short or no branch length separating them, we assumed that the unknown specimen was conspecific with its sister and removed it from the analysis. If the branch length separating them was long, comparable to the branch length between known species pairs on the tree, we assumed the specimen represented a distinct taxon and was left in the analysis.

Character sampling
DNA was extracted from a single mesothoracic leg using Qiagen DNeasy® Blood and Tissue Kits (Valencia, CA, USA) following the manufacturer's protocol. Polymerase chain reactions (PCR) were performed using an Eppendorf Mastercycler ep gradient S Thermal Cycler® (Eppendorf, Hamburg, Germany) with TaKaRa Amplitaq TM (Applied Biosystems, Foster City, CA, USA). Various primers and temperature regimes were used to amplify five gene fragments including a fragment of mitochondrial DNA spanning cytochrome c oxidase I (791 aligned base pairs [bps]), tRNA-Leucine (132 bps), and cytochrome c oxidase II (251 bps) and four nuclear gene fragments. The nuclear loci were RNA polymerase II (839 bps), arginine kinase (547 bps) with one intron (445 bps), the F2 copy of elongation factor 1-alpha (730 bps) with one intron (281 bps), and opsin (421 bps) with two introns (908 bps). See Supplementary Table  S2 for a list of primers. In a few cases, introns were removed, usually from genera far removed from Melissodes that could not be aligned.
For the taxa from Cardinal & Danforth (2013) we used six loci from the Cardinal & Danforth (2013) study; three of which overlapped with the current analysis (EF1a, Opsin, and PolII) and an additional three; 18S (782 bps), NaK (1,460 bps), and Wingless (455 bps) that did not.
PCR amplicons were visualized using gel electrophoresis, cleaned with ExoSAP-ITTM (USB-Affymetrix, Cleveland, OH, USA), purified with Sephadex® G-50 (GE Healthcare, Uppsala, Sweden) to prepare for sequencing using ABI Prism Big DyeTM (v3.1; Invitrogen, Fairfax, VA, USA). Sequencing was conducted with an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) in the Biology Department of the University of New Mexico. We sequenced the amplicons in both directions (see Supplementary Table S1 for fragments that were only successful in a single direction) and the resulting data were edited in Sequencher® (Gene Codes 1999).

Data analysis
Individual gene sequences were aligned using MUSCLE (Edgar 2004) implemented in MEGA6 (Tamura et al. 2011) using default parameters. Introns were only problematic in a few cases and sections were removed only when alignment was ambiguous (Gen-Bank sequences were complete). Gaps were treated as missing data. All eight aligned loci were concatenated and organized in Mesquite (Maddison & Maddison 2018). In cases where multiple sequences for the same species were amplified and the combination of those sequences created a longer fragment than each alone, these fragments were made into contigs that were used in the final analysis, but only if there were no base pair differences in the sections that did overlap. This treatment resulted in a total of 5,345 bps from our own data and 8,059 total bps including Cardinal & Danforth data (6,293 coding, 1,634 intron, and 132 tRNA). All sequences were submitted to GenBank (Supplementary Table S1).
PartitionFinder (Lanfear et al. 2012) was used for model selection and finding the best-fit partitioning scheme. Both mitochondrial and nuclear protein-coding genes were partitioned by codon positions. Introns, tRNA, and 18S rRNA were treated as single data blocks, resulting in a total of 42 data blocks. The greedy search algorithm was used to find the best-fit scheme, which was determined by the Bayesian Information Criterion (BIC), implemented in PartitionFinder. PartitionFinder suggested 8 partitions, which were used for all subsequent analyses. Maximum Likelihood (ML) and Bayesian analyses (BA) were run in RAxML 7.2.8 (Stamatakis et al. 2008) andDOI 10.1163/1876312X-bja10015 MrBayes V3.2.6 (Ronquist et al. 2012), respectively on XSEDE (Extreme Science and Engineering Discovery Environment, https://www.xsede.org) through the CIPRES Science Gateway (Miller et al. 2011). For the ML analyses, we used the best-fit partitioning scheme recommended by PartitionFinder with the GTRCAT model applied to each partition and nodal support was evaluated using 1,000 replications of rapid bootstrapping implemented in RAxML. For the BA analyses, we also used the best-fit partitioning scheme and partition-specific models recommended by PartitionFinder and default priors, and two independent analyses were conducted each with four runs with four chains each for 100 million generations, sampling every 2,500 generations. We plotted the likelihood trace for each run to assess convergence in Tracer V1.6 (Rambaut et al. 2014), and discarded an average of 25% of each run as burn-in. The alignment and newick trees for both the ML and BA analyses were deposited in Mendeley (Wright 2019, doi: 10.17632/zjrrwkcbzd.1). Additional trees were inferred excluding the data from Cardinal & Danforth (2013) with no effect on the topology of the trees.

Phylogenetic relationships among outgroups
All tribal relationships were consistent with the phylogeny of Cardinal & Danforth (2013). Within Eucerini, Svastra was sister to Melissodes. Martinapis was sister to Melissodes + Svastra and the Eucera complex was sister to the larger clade. Florilegus + Svastrina were sister to the remainder of the Eucerini. The ML tree (Fig. 4) and the BA (Supplementary Material) consensus tree were in agreement except for the placement of Syntrichalonia. The ML tree placed Syntrichalonia as sister to the Synhalonia + Eucera clade whereas the BA tree placed it as sister to the entire Eucera complex.

Phylogeny of Melissodes
We recovered the genus Melissodes as monophyletic (Fig. 5)      have very sparse branching on their scopae, M. micheneri LaBerge and M. moorei Cockerell both have sharply acute triangular pygidial plates and tessellate galea in the females, and finally M. lutulentus LaBerge and M. pallidisignatus Cockerell both have irregular size and spacing of punctures in the T2 interband region. Because of these morphological similarities as well as the slightly better resolution, we present the ML analysis here (Fig. 5).

Discussion
This study represents the first comprehensive molecular phylogeny of the long-horned bee genus Melissodes. Few papers have been published on the Eucerinae in general. Cardinal & Danforth (2013) presented a Bayesian phylogeny of all bees that included four eucerine representatives (including one Melissodes) that placed Eucerini sister to Ancylaini. All relationships among taxa in the current study and in Cardinal & Danforth (2013) are equivalent, which is not suprising considering we used some of their data. In 2014, Praz & Packer presented a more in depth look at the Ancylaini and its relationship to Eucerini, which included 13 eucerine genera (again with one Melissodes). The relationships among taxa that are included in both this study and Praz & Packer are congruent except for the placement of Tetraloniella relative to Eucera. Finally, Dorchin et al. (2018a) provided a thorough revision of the Eucera complex, sinking many genera to subgeneric standing and erecting one new genus (Dorchin et al. 2018b). The Dorchin et al. study (2018a) included eight eucerine genera with four representatives of Melissodes. Both the Dorchin et al. (2018a) tree and the ML tree from this study place Syntrichalonia as sister to the Eucera + Synhalonia subgenera, whereas the BA tree from this study places it outside the entire Eucera complex.
The sister relationship among the subgenera M. (Ecplectica) and Melissodes s.s. is not surprising, however the original assessment of LaBerge (1956aLaBerge ( , 1961  Although the pre-Hennigian approach that LaBerge used was inadequate for inferring the relationships among the subgenera, his placement of most species in their respective subgenera is consistent with this phylogeny with a few exceptions. Melissodes (Apomelissodes) whereas the male antennal length is short and contradicts this placement. These two small clades are relatively well supported and, in a future revision, will require a more thorough study of their morphologies.
Melissodes (Eumelissodes) tristis and M. (Eumelissodes) pexa were unexpectedly placed as sister to M. (Eumelissodes) + M. (Callimelissodes). Melissodes (Eumelissodes) tristis is a widespread and morphologically distinctive species. It is a hyper-generalist, collecting pollen from at least 26 genera in 13 plant families (Wright 2018). The females have a broadly hyaline apical edge on T1, a single row of dark erect hairs posterior to the apical hair band on T3, few to no punctures in the interband zone of T2 and the disc of the scutellum, and a shiny boss on the clypeus. In the males, the antennal configuration is consistent with M. (Eumelissodes) but the clypeus is entirely dark and they do not share the synapomorphic M. (Callimelissodes) trait of having a broadly convex hyaline apical edge of S4. Nor does M. tristis share the ecological trait of oligolecty on Asteraceae as are many of the M. (Callimelissodes) and M. (Eumelissodes). Instead, its dietary patterns are more consistent with those of a broad generalist, as are most of the Melissodes s.s. (Wright 2018). Finally, Melissodes (E.) pexa is a morphologically typical but rare species with very little known about its biology. Retention of these two species in M. (Eumelissodes) would render it paraphyletic but the support for these two nodes is not strong, and without other compelling information, it would be premature to make any taxonomic changes based on molecular evidence alone.
Most likely 'nspwhite', 'mystery', 'affpersonatellus' and 'unk15' represent undescribed species because they are morphologically distinct from the known species. 'Ec_sp' and 'Unk11' may be described species, but we could not confirm the species identifications because of lack of reference material or lack of setae, respectively. Type specimens of M. (Ecplectica) and a group of small morphologically similar M. (Eumelissodes) would have to be examined to confirm this, but revisionary work was deemed outside the scope of this study.
Although individual small clades within M. (Eumelissodes) are well supported and consistent between the ML (Fig. 4) and the BA (Supplementary Figure S1) analyses, the backbone is poorly supported and the relationships among the smaller clades are unresolved.
Combining results from Cardinal & Danforth (2013), Praz & Packer (2014), Dorchin et al. (2018a) and this study, a comprehensive picture of the tribe Eucerini and its placement within the larger Eucerinae is coming to light. A meta-analysis combining all these data with exemplars from Gaesischia and Cubitalia and denser taxon sampling for Svastra, Melissoptila, Thygater, and Alloscirtetica would more fully round out our understanding of the tribe Eucerini.
According to the divergence time estimate analysis by Cardinal & Danforth (2013), the Melissodes + Svastra clade diverged from the Eucera clade 13.74 million years ago. If we accept this value, we can infer that Melissodes represents a very young radiation that occurred from the Middle to Late Miocene. This was a period of great change in North America. Although the timing of major events such as mountain, grassland and desert formation in North America are still debated, it is clear that during the Miocene, North America saw major changes including increased aridity and cooler temperatures (Wilson & Pitts 2010). The grasslands of North America may have started as patchy, isolated habitats as early as the Early-Middle Eocene, but the expansive grasslands that are notable today, may have formed as late as the Middle-Late Miocene (Strömberg 2011). Also, the warm deserts of southwestern North America may have formed as recently as 10,000 years ago (Wilson & Pitts 2010). Although the Asteraceae groups that thrive in these habitat types evolved much earlier, both the North American Clade of Astereae Noyes & Rieseberg 1999) and the Helianthus Alliance (Baldwin 2009) have similar patterns of a North American origin with disjunct distributions in South America. In fact, analyses of both groups point to a southwestern origin in the warm deserts of North America with a rapid Oligocene-Early Miocene diversification (Funk et al. 2009). The grasslands and arid regions, where Asteraceae thrive, were patchy and probably ephemeral as drier climates allowed grasslands and parklands to expand while the forests contracted (Strömberg 2011). These collectively suggest that the diversification of Melissodes followed the formation of the North American deserts and grasslands and the diversification of the two largest North American clades of Asteraceae. The information gained in this study can be used as a basis for future studies on the evolution of this group of important pollinators.