A molecular phylogeny of Glaphyridae (Coleoptera: Scarabaeoidea): evolution of pollination and association with ‘Poppy guild’ flowers

Among Scarabaeoidea, pollen feeding occurs in two major lineages, pleurostict Scarabaeidae and Glaphyridae. Here we infer for the first time the phylogeny of the scarabaeoid lineage Glaphyridae (Coleoptera) based on molecular data using partial gene sequences for 28S rRNA, cytochrome oxidase I (cox1) and 16S rRNA (rrnL) for 41 species. Based on the resulting tree topology, we inferred the timing of the origin of pollination and of their coevolution with different flower host taxa, with particular focus on the prominent red‐coloured ‘poppy guild’ flowers. All genera of Glaphyridae that were sampled with multiple species were recovered as monophyletic. According to this analysis, the origin of Glaphyridae was around 140 Ma, while crown group divergence was dated to have occurred c. 112 Ma. Pollen feeding originated in Glaphyridae only once and much later than in other important pollinator groups, between 97 and 67 Ma. According to the reconstruction of ancestral feeding traits, Asteraceae (Cicharioidae) were the first hosts of Glaphyridae. Presumably, a further adaptive radiation was triggered by feeding on and pollination of red flowers (poppy guild) which arose at a later stage. It occurred for the first time between 30 and 40 Ma, whereby the clades that use red Ranunculaceae (Pygopleurus spp.) are older than clades using exclusively red Papaveraceae (Eulasia spp.) (25–30 Ma). The rather young age of red Ranunculaceae would imply that Pygopleurus species only subsequently used red Ranunculus species as flower hosts, and that a broad parallel host shift probably from red Papaver spp. to red Ranunculus asiaticus has occurred rather recently.

In the eastern Mediterranean region, the poppy guild flowers blossom in early spring, as early as February in the Jordan valley, in a period when most efficient insect pollinators are not yet active. Thus, the 'early-spring phenology' of some species of Glaphyridae plays a very important role in pollination in the absence of other abundant insect pollinators (G. Sabatinelli, 2006-2019, personal observations, in Israel, Jordan, Syria, Lebanon, Turkey, Greece, Armenia, Iraq, Iran andTunisia). In March-May, the phenology of Glaphyridae and red flowers overlaps in a temporal succession of plant and beetle species that supports this close plant-pollinator interaction (Dafni et al., 1990).
The strong association of Pygopleurus species with poppy guild flowers is enabled by the presence of red-sensitive photoreceptors in addition to green-and UV-sensitive ones that are more common among anthophilous beetles (Martínez-Harms et al., 2012). Given the importance of floral colour as a cue for Glaphyridae, the presence of a red-sensitive photoreceptor in their visual system can be understood as a determinant character involved in the recursive processes that conserve this form of plant-pollinator interaction.
But flowers are also mating sites for insects (Thornhill & Alcock, 1983;Barth, 1985). This is the case for many species of Glaphyridae, whose adults feed on floral pollen and mate on flowers (Keasar et al., 2010). Several species of these beetles rely strongly on visual cues to find particular flowers, as receptive females on flowers are highly dispersed in space and to find them may be difficult for males, especially if the density of females is lower than that of flowers. Even if we assume pheromones as an essential mechanism in partner finding (Robbins et al., 2006), the chance of finding a mate may improve if both males and females orient to a subset of the flowers, and visit them preferentially. Thus, the evolutionary scenario of the pollinator-flower relationship could be two-fold.
The uniqueness of the glaphyrids has been recognized for a very long time. Nevertheless, the status of the group as well as the phylogenetic position of the Glaphyridae within the Scarabaeoidea has been the subject of debate (Machatschke, 1959;Zunino, 1988). The group is considered to represent 'intermediate' forms between ancient and modern lineages of Scarabaeoidea. Browne & Scholtz (1995) found the Glaphyridae to be the monophyletic sister group of the trogid subgroup [Trogidae, Bolboceratinae (Geotrupidae) and Pleocomidae] based on characters of wing articulation (d 'Hotman & Scholtz, 1990;Scholtz, 1990;Browne & Scholtz, 1995), while in the phylogenetic study of Ahrens et al. (2014) using mitochondrial and ribosomal DNA, results indicate Glaphyridae to be sister to Pleurosticti. Two Chilean genera (Lichnia Erichson and Arctodium Burmeister) previously included in Glaphyridae ended up among Melolonthinae lineages in a phylogenetic analysis based on 18S and 28S rDNA (Smith et al., 2006) and were consequently removed from Glaphyridae; they now comprise the melolonthine tribe Lichniini Burmeister, 1844 (Hawkins, 2006).
Except for a few species, the life history of most glaphyrids is poorly known. Adults are often brightly coloured, densely setose, active diurnally, and strong flyers. Adults of many species have coloured bands of setae on the abdomen and resemble bumble bees; they are feeding on pollen but some species do not feed at all. Pollen feeders either stay on the flower surface (surface feeder) or dig deeply into the bottom of the flower (deep feeder). Larvae feed on roots or on detritus (Medvedev, 1951(Medvedev, , 1952Ritcher, 1966;Westcott, 1976;Carlson, 1977Carlson, , 1980, with some species [Lichnanthe vulpina (Hentz, 1826)] considered to be pest species (Carlson, 1977(Carlson, , 1980Averill & Sylvia, 1998).
Here, we investigate the phylogeny of Glaphyridae using mitochondrial and ribosomal DNA markers. Based on this first phylogenetic hypothesis of the group, using the most comprehensive sampling so far available, we aim to reconstruct the evolution of feeding behaviour of these beetles and its time frame, with particular reference to host association. We are particularly interested in what drives host choice, and whether we can identify a closer evolutionary link to the red 'poppy guild' flowers.

Materials and methods
Taxon sampling, DNA extraction, and DNA sequencing Forty-five populations of 42 species of Glaphyridae were sampled from various biogeographic regions (Table 1). The analysis additionally included nine outgroup taxa. DNA was extracted from thoracic leg muscle tissue using Promega Wiz-ardSV extraction plates. Following DNA extraction, beetles were dry mounted. Vouchers are deposited at the Zoological Research Museum A. Koenig (ZFMK) in Bonn. Phylogenetic analysis was based on three markers: mitochondrial gene regions including cytochrome c oxidase subunit 1 (cox1) and 16S ribosomal RNA (rrnL) as well as nuclear 28S rRNA (domains D3-D6). PCR and sequencing was performed using primers Pat and Jerry for cox1 and 16Sar and 16S B2 for rrnL (Simon et al., 1994); 28S rRNA was amplified using primers FF and DD (Monaghan et al., 2007). Sequencing was performed on both strands using bigdye v. 2.1 and an ABI3730 automated sequencer. Sequences were edited manually using sequencher v.4.8 (Genecodes Corp., Ann Arbor, MI, USA). All sequences used in the present study have been submitted to GenBank (Table 1).
Host plants, in particular for Israel and neighbouring countries, were identified using Shmida (2005) and we consulted A. Shmida directly in case of doubts.

Phylogenetic analysis and divergence time dating
Multiple sequence alignment was conducted with the divide-and-conquer realignment technique as implemented in saté v.2.2.7 (Liu et al., 2012) using mafft v.5.8 (Katoh et al., 2002(Katoh et al., , 2005 as aligner and muscle v.3.7 (Edgar, 2004a,b) as merger. partitionfinder was used to identify best-fitting substitution models and partition scheme under linked and unlinked branch lengths (Table S1; Lanfear et al., 2012Lanfear et al., , 2014 for the data which was initially split into rrnL, 28S, and three codon positions of cox1. The optimal setting was chosen based on the Bayesian information criterion. Bayesian analysis of phylogenetic relationships and divergence times was conducted with beast 2.3.0 (Bouckaert et al., 2014) using the partitions inferred by partitionfinder (Nylander et al., 2004;Brandley et al., 2005). A lognormal relaxed clock model (Drummond et al., 2006) was applied, setting the number of discrete rates to the number of branches. Two fossil priors were available for tree calibration (Fig. 2): Cretoglaphyrus Nikolajev spp., being the oldest fossils of Glaphyridae (Krell, 2000(Krell, , 2006Nikolajev, 2005;Zhao et al., 2016) [node A, maximum age (range) for the origin of the Glaphyridae: 145.5-140.  (Zhao et al., 2016) cannot be associated with any of the extant genera due to their poor preservation and were thus not considered for calibration here. The North American fossil Lichnanthe defuncta whose genus is endemic to northern America, however, could be used as similar geographical distributions indicated the systematic position of the fossil.
Exponential priors were used to implement fossil calibration points for node A (offset 140.2 Ma) and node B (offset 33.9 Ma). The calibration distribution means were estimated by gamma distributed hyperpriors [node A: alpha = 10, beta = 0.5, i.e. 95% highest posterior density (HPD) for exponential prior mean varies about 5 Ma; node B: alpha = 10, beta = 2, i.e. 95% HPD for exponential prior mean varies about 20 Ma]. We ran four independent analyses with 1 × 10 8 generations; parameters and trees were logged every 10 000 generations to produce 10 000 samples. Convergence of parameters and stationarity of runs were assessed for each run separately in tracer v.1.6 (Rambaut & Drummond, 2013). Log files and species trees from independent runs were subsequently combined using logcombiner 2.3.0, after removing a burn-in of 10% and checked again in tracer. Trees were summarized using treeannotator 2.3.0. The .xml-file for a single Markov chain Monte Carlo run is provided as File S1. Model-based phylogenetic analyses under maximum likelihood (ML) were performed in iq-tree (v.1.6.7; Nguyen et al., 2015;Chernomor et al., 2016). Partition schemes and substitution models were set according to the partitionfinder analysis. Individual evolutionary rates for each partition were allowed. The best tree out of 20 full tree searches was chosen by its likelihood. Node support was estimated by 1000 replicates of ultrafast bootstrapping (Hoang et al., 2018). The tree was rooted according to the beast analysis (see later).

Phylogenetic comparative analysis
For each glaphyrid species feeding type, flower taxon host, and host colour was coded. Colour preferences are also based on results from plate trap experiments (Dafni et al., 1990). As no published data were available for many taxa, trait coding was principally based on our extensive field observation (G. Sabatinelli 2006(G. Sabatinelli -2019. Maximum parsimony reconstructions (Fitch algorithm) of ancestral feeding type and flower taxon and colour were done with phangorn v.2.4.0 (Schliep, 2011). Ambiguous character states were coded as such and inapplicable traits of non-feeders were treated as extra character state. As body size (as a proxy of body weight) might play a role in host choice in relation to the robustness of inflorescences and their stems, we analysed correlations of body size in relation to feeding behaviour, pollinated plant taxon, and the colour of the florescence with phylogenetic generalized least-squares regression (PGLS) in caper v.1.0.1 (Orme et al., 2018), taking the phylogenetic relationship of investigated species as inferred with beast into account. As no reliable mean values for body size were available, all analyses were done for minimum and maximum body sizes of all species. Because caper is not able to treat ambiguous character states, all ambiguities were treated as missing data. This seemed to be the cleanest way to handle such data with available phylogenetic comparative methods.

Phylogenetic analysis and divergence time dating
In the tree from Bayesian inference (BI) Glaphyridae were recovered as monophyletic and sister to Ochodaeus (i.e. Ochodaeidae). All genera of Glaphyridae were monophyletic.  The iq-tree analysis had a very similar generic topology, with exception of the position of Athypna meles (Fabricius), being nested within Pygopleurus (Fig. S1). Trichopleurus was positioned here as sister to the subgenus Eulasia + E. nitidicollis (Reiche), while in BI it was nested within two Rudeulasia species clades.

Phylogenetic comparative analysis
Non-feeding, being present in Amphicoma and Lichnathe (which also do not mate on flowers), was recovered as ancestral in the projection of the feeding mode (Fig. 3a) on the tree topology of BI. Superficial and deep feeding were inferred as derived characters with a unique origin. The first is encountered in the clade (Anthypna + Pygopleurus) + Eulasia, whereas the latter is restricted to the genus Glaphyrus. The type of feeding is also widely consistent with the morphology of mandibles of the species, showing an unarmed, blunt apex in non-feeding Amphicoma and Lichnathe, but sharp teeth at the distal part of mandible in feeding Glaphyrus, Pygopleurus and Eulasia. The only exception is represented by the superficial feeder Anthypna, showing a moderately blunt tip similar to that of Lichnathe (Fig. 3a).
By contrast, the systematic position of flower host and flower colour did not show a close link to the bumble beetle phylogeny (Fig. 3b). While several species of Eulasia and Glaphyrus feed on Cichoriodeae (Asteraceae), host use was unique for the latter genus only. Several species in Pygopleurus and Eulasia also showed multiple host groups with alternative colour properties, although the majority of included taxa are linked to red flowers of the poppy guild. Interestingly, we detected one well-supported monophyletic clade within Eulasia that mostly relates to red-coloured Papaveraceae. The reconstruction of Table 2. Trait correlation analysis between minimum/maximum body size and feeding traits/host characteristics from phylogenetic generalized least-squares regression analyses treating polymorphic characters as missing data showing coefficient of determination (r 2 ), F-statistics (F stat ), and P-value the ancestral flower colour visited by Pygopleurus was red, whereas the situation was ambiguous for Eulasia (Fig. 3c). Ancestral hosts were Ranunculaceae for the former and Asteraceae (Cicharioidae) for the latter. Asteraceae (Cicharioidae) were also inferred as ancestral hosts for all other deeper nodes (Glaphyrus + all other western Palearctic Glaphyridae; Athypna + Pygopleurus + Eulasia), while flower colour in these nodes was ambiguous.
The PGLS correlation analysis between body size and feeding traits, such as feeding behaviour, host taxonomy and colour of host flower, was only significant for plant taxonomy and flower colour, while feeding behaviour was not significant (Table 2). Analysing minimum or maximum body size did not alter conclusions. F-statistics and r 2 value for flower colour were highest for maximum body size, while plant taxonomy had slightly higher values for minimum body size.

Discussion
In this study, we present the first comprehensive phylogenetic analysis of Glaphyridae based on DNA sequences and extensive sampling of the entire group. While the sister-group relationship is so far still uncertain, due to contrasting signals and evidence (e.g. Hunt et al., 2007;Lawrence et al., 2011;Gunter et al., 2016), Ochodaeidae were recovered here as sister to Glaphyridae, similar to other previous studies (e.g. Hunt et al., 2007;Ahrens et al., 2014;Bocak et al., 2014). Although the taxonomic sampling was not designed to address this question, it seems apparent from most studies and also from most comprehensive and recent evidence (McKenna et al., 2019) that Glaphyridae belong to the early-diverging lineages of Scarabaeidae, together with Hybosoridae and Ochodaeidae.
According to our results, pollen feeding originated in Glaphyridae only once, between 97 and 67 Ma, and this is largely consistent with results from Ahrens et al. (2014). An exact timing cannot be given, as it is not clear whether the early members of the stem lineage of Glaphyrus + Pygopleurus/ Eulasia had already fed on pollen. According to current knowledge this was much later than in other important pollinators groups such as bees (c. 123 Ma;Cardinal & Danforth, 2013) or butterflies (119 Ma; Espeland et al., 2018). By contrast, in Pleurostict scarabs there were multiple (at least six) shifts from herbivory to pollen feeding . Similar to Glaphyridae, all occurred rather late in the Tertiary, from the Palaeocene to the Miocene.
The feeding type of first glaphyrid pollinators is still uncertain. Further evidence might come from a comprehensive comparative morphological analysis of all body traits between superficial and deep feeders. Deep feeding on flowers is likely to be linked to additional adaptations, such as the reinforcement of anterior legs for digging into the flowers, and the robustness of maxillary and labial palps (i.e. Glaphyrus), and thus distribution of these traits in ancestral and pollinating lineages might provide further insights. Such traits are also known in deep feeders of monkey beetles (Scarabaeidae: Hopliini; Midgley, 1992;Ahrens et al., 2011). In the latter case, males compete for females on the flowers, which might also be the case in Glaphyrus. Interestingly, body size was not correlated with feeding behaviour, possibly due to the fact that there was considerable variation among the superficial feeders.
According to the reconstruction of ancestral feeding traits, Asteraceae (Cicharioidae) were the first hosts of Glaphyridae. Presumably, a further adaptive radiation was triggered by feeding on and pollination of red flowers (e.g. poppy guild), which arose at a later stage. It occurred for the first time between 30 and 40 Ma, whereby the clades that use red Ranunculaceae (Pygopleurus spp.) are older than clades using exclusively red Papaveraceae (Eulasia spp.) (25-30 Ma). While the genus Ranunculus is known to have originated in the Oligocene , the red Ranunculus acting as host to the Pygopleurus species in the Near East has been dated to 3.9 Ma . This would mean that Pygopleurus species only subsequently used red Ranunculus species as flower hosts, and that a broad parallel host shift probably from red Papaver spp. to red R. asiaticus had occurred. Unfortunately, we had no access to dated phylogenies of Papaver (i.e. Papaveraceae) that might help in elucidating the origin of red flowers in this group in order to evaluate the incidence with Glaphyridae evolution. Nevertheless, as eastern Mediterranean Glaphyridae are among the most important pollinators of poppy guild flowers and also because considerable physiological and anatomical adaptation for vision of red colours has been found (Martínez-Harms et al., 2012), it has to be assumed that there always was a close evolutionary relationship between these flowers and Glaphyridae (i.e. species of Pygopleurus and Eulasia).
To explore the question of what determines the choice of flower host, we investigated the correlation between body size (as a proxy of body weight) and flower taxonomy, feeding behaviour and flower colour. Body weight (i.e. size) could play a role in host choice in relation to the robustness of flowers and their stem, in particular when considering that flowers are not only a food resource of Glaphyrids but also a mating site (Keasar et al., 2013). Sometimes many specimens are found on a single florescence. The correlation of body size with flower taxonomy and colour confirmed this expectation, as flowers of the same taxon and colour tend to have similar morphology, being either more or less robust. As no reliable mean values for body size were available, all analyses were done with known minimum and maximum body sizes of the species. Both had similar patterns of correlation, although the correlation for maximum size was clearer, probably indicating the determinant measure for this relationship. Species on red flowers, i.e. poppy guild flowers, are generally smaller. Keasar et al. (2013) argue that adaptation on poppy guild flowers would be a mode of aggregation behaviour, which is in turn promoted by plant-mediated chemical protection of beetles from predation (Keasar et al., 2013). Most poppy guild taxa include phenologically early inflorescences, and thus pollinators might not be too abundant. Therefore, it is likely that plants adapted to new pollinator resources that are available in early spring, and that bumble beetles switched to new food resources.
Although the evolution of pollination and the coevolution with red poppy guild flowers make up only one chapter of the evolution of Glaphyridae, it is the one with the most important impact on ecosystems by the group. The evolutionary plasticity evident in adopting new food resources in relatively young geological timescales is impressive and, at the same time, surprising for such a comparatively ancient lineage, and merits further investigation at a finer scale in order to better explore the mechanisms of this interesting coevolutionary relationship.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.  Table S1. Partition scheme and substitution models that were inferred by partitionfinder. The results of the analysis with linked branch lengths is shown, which was preferred based on the BIC. File S1. beast extensible markup language (xml) file for one out of four MCMC-run repeats of the dating analysis.