Our morphological data set contains 18 of the 22 Phlaeopterus species, and nine of these are represented in our molecular dataset (Table 2). The three Chinese species recently transferred from Lesteva to Phlaeopterus and the new species described by Shavrin (2020), Phlaeopterus rufimarginatus (Rougemont, 2000), Phlaeopterus ruzickai (Shavrin, 2015), Phlaeopterus septemmaculatus (Rougemont, 2000), and Phlaeopterus hromadkai Shavrin, 2020, were not included in our analyses, which were based on data obtained before publication of Shavrin (2020). However, our dataset includes one Palearctic species, Phlaeopterus czerskyi (Shavrin, 2001).

Wherever possible, multiple exemplars of each Phlaeopterus species were included in our molecular dataset, with specimens selected from the widest available geographic range for each species. Geographic ranges were estimated based on study of borrowed museum specimens. Of over 2,000 museum specimens databased (http://arctos.database.museum/saved/Phlaeopterus) each species was represented by an average of 142.2 specimens for these range estimates (range: 2–325 specimens/species). We prepared the map in Fig. 3 using these data with SimpleMappr (Shorthouse 2010). Repository details for these specimens can be found in Mullen et al. (2018). Our outgroup for the morphological dataset is composed of four Lesteva and one Unamis Casey, 1894 species, the latter corresponding to the same Unamis species for which we could obtain a DNA barcode sequence. Our outgroup for the molecular dataset is composed of all the publicly available Lesteva (six species) and Unamis (one sequence identified to species and two identified to genus) DNA barcode sequences that were available at the time, and those we generated, on BOLD. Lesteva and Unamis were selected as outgroup taxa based on their morphological similarity to Phlaeopterus (Moore and Legner 1979) and because all three genera belong to the tribe Anthophagini (Newton et al. 2000).

Our morphological dataset contains 35 external and 5 male genitalic characters, the majority of which derive from the unpublished work of J.M. Campbell. Character 35 was modified from Moore and Legner (1979). Some of Campbell’s original characters were modified after our study and measurements of specimens led us to disagree with his codings (e.g. character 9, interfacetal setae of eye, was changed from a three-state character to a two-state character). Codings reflect hypotheses of homology among species, although some characters are notable autapomorphies. When doubt existed about the homology of states in certain taxa those taxa were coded as ‘missing/inapplicable’ for those characters. Continuous characters such as ratios of lengths were evaluated to find gaps among species which were used to assign states. Within-species polymorphism was handled by assigning the majority state to that species when the state was clearly in the majority or assigning a separate state for the presence of polymorphism (e.g. character 30). When no state was in the clear majority, the species was coded as polymorphic with both states (for instance, P. lagrandeuri Hatch, 1957 in character 27). Polymorphism of length characters within a species was handled by using the average of multiple measurements. All characters of potential phylogenetic value were included even if some were missing data.

We observed characters with a Leica M165 C stereomicroscope (Leica Microsystems, Wetzlar, Germany), and coded morphological data in Mesquite 3.6 (Maddison and Maddison 2016). Our morphological data matrix and trees are archived at http://purl.org/phylo/treebase/phylows/study/TB2:S27393 and among Supplementary files S1, S2, and S3. Length data for character 16 were analyzed statistically using R version 4.0.0 (R Core Team 2020) in RStudio version 1.2.5042 (RStudio Team 2020) during the character coding process.

We coded characters from type specimens whenever possible, and specimens determined primarily by J.M. Campbell belonging to, or on loan to, the University of Alaska Museum Insect Collection from the California Academy of Sciences, San Francisco, California (David H. Kavanaugh, Jere Schweikert) and the Canadian National Collection of Insects, Ottawa, Ontario, Canada (Patrice Bouchard, Anthony Davies).

2.2.1. List of morphological characters

Character descriptions follow for our morphological data. All characters are unordered and equally weighted. Figures are primarily from Mullen et al. (2018), cited below as “fig. M#” or “figs M#,#” using a capital M to distinguish them from figures herein, cited as Fig. #.

1. Maximum size: (Figs 1, 2) (0) less than 5.0 mm; (1) greater than 6 mm. This character represents J.M. Campbell’s hypothesis that the large-bodied species form a monophyletic group. In Figs 2, 4–7, the large-bodied species, corresponding to state 1, are colored red, blue, or green.

Figure 4. 

Evidence that large body size evolved twice: a clade of large bodied species (red and blue colors) and a separate origination of large body size in P. elongatus (green). Concatenated Bayesian analysis of 40 morphological characters and 404–654 bp of COI with posterior probability and maximum likelihood bootstrap support values near each branch with 18 Phlaeopterus ingroup species, 1 Unamis and 4 Lesteva outgroup species. Topology is a Bayesian 50% majority rule consensus phylogram. Clades with ML bootstrap values less than 50% are indicated with a “–“. Taxon names in blue and red are species with maximum body sizes over 6.5 mm. Blue indicates the largest species, those with maximum sizes over 8 mm, illustrated by A) P. bakerensis; names in red are species with maximum body sizes 6.5 to 7.9 mm with wide bodies, illustrated by B) P. kavanaughi; the name in green is for C) P. elongatus with a maximum body size of 6.4 mm, with an elongate body; names in black are species with maximum body sizes below 5 mm, illustrated by D) P. obsoletus. Beetle inset illustrations are scaled to approximate relative body sizes.

2. Maximum size: (Figs 1, 2) (0) less than 7.8 mm; (1) greater than 8.1 mm. Character one and two reflect the hypothesis that the largest bodied species with maximum lengths over 8.1 mm are a clade within or sister to the large-bodied species that are over 6 mm in length. Because small body size is plesiomorphic, characters one and two do not represent double-weighting of small body size. In Figs 2, 4–7, the largest bodied species, corresponding to state 1, are colored blue.

3. Body length to width ratio: (0) wide, length to width ratio under 3 (figs M1, 2B–D, 3A, C, D, 4, 5); (1) elongate, length to width ratio greater than 3 (figs M2A, 3B).

4. Ocelli: (0) present (figs M34A–D, 35A–D); (1) absent (fig. M34E, F).

5. Color of elytra: (0) dark brown, light brown, or dark reddish brown, to black (figs M2, 4); (1) light reddish to yellowish brown (figs M1B, 3C); (2) bicolored, part dark to light brown, part yellowish to reddish yellow (fig. M3D).

6. Antocellar foveae (dorsal impressions between eyes) shape: (0) elongate narrow (fig. M27); (1) elongate wide (figs M34C–F); (2) oval to circular (figs M35A–D); (3) reduced (fig. M34A, B).

7. Nuchal constriction: (0) distinct, post-ocular region clearly divided into temple and neck with impression continuing across midline of neck (fig. M27); (1) vague (figs M34A, B, 35A–D; (2) absent (fig. M34C–F).

8. Interantennal impression: (0) present (fig. M35A–D); (1) obsolete. This is a transverse impression that runs across the frons between the bases of the antennae.

9. Interfacetal setae of eye: (0) present on most of eye (often glabrous in extreme dorsal portion) (fig. M35E, F); (1) absent from most of eye (glabrous in dorsal half, reduced in number or glabrous on ventral half) (fig. M36A–F).

10. Shape of labrum: (0) 0.3–0.5 times as long as wide (fig. M31A, B, E, F); (1) less than 0.3 times as long as wide (fig. M31C, D).

11. Dorsal surface of labrum: (0) without micropores (fig. M31A); (1) with micropores limited to setose region (fig. M31B–G); (2) with micropores across labrum, many posterior to setose region (fig. M31H).

12. Minimum width of gula: (0) narrow, less than 0.2 times as wide as mentum (fig. M41E, F); (1) wide, more than 0.2 times as wide as mentum (fig. M42).

13. Base of epipharynx: (0) with oblique, parallel rows of fine ridges (fig. M33A–D); (1) without rows of ridges (fig. M33E–H).

14. Apical portion of epipharynx: (0) with spines extending to apical margin (fig. M33A, B, E–H); (1) smooth (fig. M33C, D).

15. Center of epipharynx: (0) with many widely spaced large spines (fig. M33A, B); (1) with many closely spaced small spines (fig. M33C–H).

16. Maxillary palpi (length of fourth vs. third palpomere): (0) apical palpomere on average at least 3.5 times as long as third (figs M30E, 32A); (1) apical palpomere on average 3 to 2.55 times as long as third (fig. M32B, C); (2) apical palpomere on average less than 2.47 times as long as third (fig. M32D). Five specimens of each species (except P. czerskyi – we only had two) were measured which revealed an inverse relationship between a species’ average 4^th vs 3^rd maxillary palpomere ratio and increasing body length (R² = 0.5742, p = 0.00027). States were chosen to correspond to breaks in the data. Larger bodied Phlaeopterus have 3^rd palpomeres considerably longer than smaller bodied species, relative to the length of their 4^th palpomeres.

17. Fine cilia at base of hypopharynx: (0) in oblique rows (fig. M40A, B); (1) not in rows (figs M40C–F, 41A–C).

18. Micropits of antennae: (0) present and with numerous non-protruding papilliform structures (fig. M37B–E); (1) absent or reduced to tiny pits (fig. M37A); (2) present and with pore-like openings (fig. M37H, I); (3) present and with numerous protruding papilliform structures (fig. M37F, G, J, K).

19. Molar area of mandible: (0) with fine parallel rows of short setae or setiform projections (fig. M28A, B); (1) with no setae or setiform projections (fig. M28C–F, 29A–C). It is unknown if these structures are true setae or setiform projections.

20. Prosthecal fringe of mandible: (0) does not extend apicad of the subapical mandibular tooth (fig. M28A, B); (1) extends beyond subapical mandibular tooth (figs M28C–F, 29A–C).

21. Mesal margin of mandible: (0) not curved, angulate, or excavate apicad of the molar area (fig. M28A–D); (1) curved, angulate, or excavate apicad of the molar area (figs M28E, F, 29A–C).

22. Lateral area of pronotum: (0) with no, or only trace of impression (fig. M15A, B); (1) moderately impressed (figs M15C–H, 17E); (2) deeply impressed, foveiform (fig. M16).

23. Base of pronotal disk: (0) without pair of transverse shallow impressions (figs M15, 16, 17A, B, E); (1) with pair of transverse shallow, sometimes confluent, impressions (fig. M17C, D).

24. Pronotal width to head width: (0) less than 1.5 (figs M2A, B, 3A, B, D, 4B, C); (1) more than 1.55 (figs M1A–C, 2C, D, 3C, 4A, D 5A, B).

25. Pronotal lateral margins posterior to lateral impressions: (0) not deflexed (figs M15C–E, 16, 17A, B); (1) deflexed posterior to deeply impressed lateral impressions (fig. M17C, D); (2) deflexed posterior to obsolete or vaguely impressed lateral impressions (figs M15F–H, 17E).

26. Anterior ridge of mesoventrite: (0) with no projection or only a small median posterior lobe (figs M19A–D, 20F); (1) prolonged posteriorly into a long projecting median tooth (figs M19E, F, 20A–E).

27. Second tooth posterior to first on midline of mesoventrite: (0) absent (fig. M19E); (1) present (fig. M20E). Those without a first tooth also lacked a second tooth so were coded inapplicable for this character. This second tooth is easily seen in a lateral view of most specimens that have been card mounted.

28. Mesoventral carina: (0) present along midline (figs M19A–D, 20E); (1) vague or absent (figs M19E, F, 20A, B, F).

29. Apical margin of elytra: (0) convex or subtruncate (figs M1, 2A, C, D, 3, 4B–D, 29F); (1) prolonged at suture only in females (figs M2B, 29G); (2) prolonged at suture in both sexes (figs M4A, 29E).

30. Metathoracic wings: (0) fully developed; (1) reduced or absent; (2) fully developed in most individuals but reduced in some.

31. Glabrous portion of mesotibia: (0) absent; (1) present (fig. M18).

32. Apical tooth on metatrochanter: (0) absent; (1) present (fig. M41D).

33. Humeral angles of elytra: (0) convex, epipleural carina not projecting (figs M1B–D, 2–5); (1) more rectangular, epipleural carina projecting and explanate (fig. M1A).

34. Shape of wing-folding spicule patches on tergite 5: (0) round and widely separated; (1) broadly oval and narrowly separated or so close as to appear combined into a single transverse band (fig. M17F–H); (2) absent. The presence of these spicules is logically correlated with the presence of wings, and indeed, P. czerskyi is brachypterous and lacks these spicules, however, P. houkae is brachypterous and has these spicules.

35. First metatarsomere: (0) shorter than ultimate tarsomere; (1) subequal to or longer than ultimate tarsomere (figs M1–5).

36. Paramere length to median lobe length ratio: (0) 1.16 or less (figs M21B–F, 22, 23, 24, 26); (1) 1.19 or greater (figs M21A, B, 25B, D). Measurements were taken from the figures in Mullen et al. (2018) with averages used for species with multiple figures except the subspecies of P. castaneus, which were not averaged. Measurements were taken from the median junction of the parameres to the paramere and median lobe apices.

37. Internal sac shape (inverted): (0) not rectangular (figs M21, 22, 23, 24A, B, 25B–D); (1) rectangular wide (figs M24D, 25A); (2) rectangular elongate (fig. M24C).

38. Carina of median lobe apex: (0) absent (figs M21A–D, 22B–D, 23, 24A, C, D, 25, 26); (1) present (figs M21E, F, 22A, 24 B).

39. Internal sac spinature: (0) few to no microspinules; (1) sclerites and many sparse microspinules (fig. M21A–D); (2) covered with microspinules (figs M21E, F, 22–26).

40. Internal sac microspinule pattern: (0) not denser basally (figs M21, 22, 23B–D, 24A, C, D, 25, 26); (1) denser basally (figs M23A, 24B). Those with few to no microspinules were coded inapplicable for this character.

Our molecular dataset contains 164 partial sequences of the mitochondrial gene cytochrome c oxidase subunit 1 (COI), representing nine Phlaeopterus species, six Lesteva species, and one Unamis species (Table 2). We sequenced COI fragments from previously unsequenced specimens and added 107 public sequences from BOLD after filtering out sequences less than 400 bp in length, those flagged as possible misidentifications or contaminations, those with 5 or more ambiguous character states, and those with stop codons. These included data from Hendrich et al. (2015), Pentinsaari et al. (2014), and Sikes et al. (2017). We also included one non-public sequence for Phlaeopterus loganensis provided to us with permission by the Canadian National Collection of Insects, Arachnids and Nematodes (CNC), Ottawa, Ontario, Canada (Table 2).

Sequences ranged from 403 to 654 bp in length. We extracted DNA from whole hind legs using Qiagen DNeasy extraction kits following the “Purification of Total DNA from Animal Tissues” protocol in the DNeasy Blood and Tissue Handbook that came with the extraction kit. We amplified a 658-bp region of COI using standard COI barcoding forward and reverse primers, LCO-1490 and HCO-2198, respectively (Folmer et al. 1994). We ran all PCR reactions at 25-μl volume. Typical PCR reaction solution included: 12.5 μl GoTaq Green Master Mix, 1.0 μl forward primer, 1.0 μl reverse primer, 8.5 μl H₂O, and 2.0 μl template DNA. We used the thermocycler protocols of Folmer et al. (1994) except where noted. Typical thermocycler protocol included a single cycle of 95°C followed by 35 cycles of 95°C for 30 seconds of denaturation, 45°C for one minute of annealing, and 72°C for two minutes of elongation, followed by a single extension of 72°C for 10 minutes.

We viewed sequence data with 4Peaks (Griekspoor and Groothuis 2005) and aligned sequences by eye in MacClade 4.08 and Mesquite 3.6 (Maddison and Maddison 2005, 2016). We created consensus sequences from bidirectional reads, aligned by eye, then checked our alignment by translation to amino acids and alignment by amino acid codon position using the “minimize stop codons” option in Mesquite to find the reading frame. The final alignment was free of stop codons and aligned with a published Lesteva longoelytrata (Goeze, 1777) COI sequence (GenBank accession: KM442270.1, Hendrich et al. 2015) in both nucleotide and amino acid alignment. Voucher specimens used for DNA extractions can be found in the following insect collections: University of Alaska Museum (UAM), Canadian National Collection of Insects (CNC), California Academy of Sciences (CAS), Santa Barbara Museum of Natural History (SBMNH), Field Museum of Natural History (FMNH), Smithsonian Museum of Natural History (NMNH), University of Idaho William F. Barr Entomological Museum (WFBM), and Brigham Young University Monte L. Bean Life Science Museum (BYUC). Digital records of specimen collection data we used to generate new sequences are available through the UAM Arctos database (https://arctos.database.museum/saved/PhlaeopterusBold). Our molecular data matrix and trees are archived at http://purl.org/phylo/treebase/phylows/study/TB2:S27393 and as Supplementary files S1, S2, and S3. GenBank and BOLD accession / process ID numbers for our sequence data are provided in Table 2.

We tested the monophyly of prior taxonomic hypotheses, as well as unpublished taxonomic groupings (Table 1) using estimated posterior probability (PP) support values obtained from our Bayesian analyses. Taxon relationships that did not occur in any of the sampled trees in an MCMC run were assumed to have a PP of < 1/number of trees sampled (Miller et al. 2002; Holder and Lewis 2003). In our analyses, such clades have a PP < 0.0001. Of particular interest was the unpublished hypothesis of Campbell ‘Phlaeopterus sensu stricto’ which states that all the large-bodied species are descended from a large-bodied common ancestor, i.e. large body size evolved once in the genus.

To provide an additional test for this hypothesis, we used Bayes factors (Kass and Raftery 1995; Ronquist et al. 2020) to compare the means of the marginalized likelihood scores of two stepping stone MrBayes (Xie et al. 2011) analyses, one with a topological partial constraint enforced that required all large-bodied species except P. bakerensis, which was free to join anywhere, to form a clade, and one with a partial constraint that required all large-bodied species except P. elongatus Mullen and Campbell, 2018, which was free to join anywhere, to form a clade (the constrained clade had 12 taxa in both cases). This approach includes one analysis forcing P. elongatus to belong to the clade of large-bodied species and one allowing P. elongatus to join anywhere but keeps the size of the clade the same in both analyses. For the COI data, we balanced the topological constraint so that there were 84 constrained sequences in both analyses. Comparing two analyses with the same number of constrained taxa eliminates the differences in the topology priors that would result from comparing a constrained to an unconstrained analysis, or to a constrained analysis with a different number of constrained taxa. Bergsten et al. (2013) documented the important influence of topology priors on Bayes factors tests of analyses with topological constraints, and ways to improve the reliability of these tests. Stepping stone runs with the morphology-only data resulted in high variance among estimated marginal likelihood scores that inexplicably fell into three bins (scores 109–112, 235–238, and 380–383). Although the results from the morphology-only Bayes factor comparisons, within bins, agreed with the conclusions found from the concatenated and COI analyses, this high variance led us to doubt the reliability of the morphology-only Bayes factors results so they are not shown.

To test an alternative hypothesis that the most recent common ancestor of the genus Phlaeopterus was large-bodied with subsequent reversals to small-bodied adults in some species, we performed ancestral character state reconstruction with Mesquite 3.6 (Maddison and Maddison 2016) to estimate the body size of the most recent common ancestor of the genus using the Mk1 model, the morphology only matrix, and the concatenated analysis tree (Fig. 4). We also explored an asymmetric two-parameter likelihood model and a step matrix parsimony model to see what ratio of gains to losses would be needed to estimate the common ancestor as large-bodied (Cunningham et al. 1998). A likelihood ratio test performed in Mesquite rejected the asymmetric two parameter model as overfitting the data (Chi-square, p = 0.9233) so we present only the results under the one parameter model. We tested the monophyly of Phlaeopterus species using the full 164-sequence COI dataset.

If a monotypic genus concept failed to match our trees due to paraphyly resulting from the missing species being unknown at the time the genus was erected, such as Motschulsky’s (1853) original concept for Phlaeopterus, we considered the concept untestable. Casey’s (1893) concept of Phlaeopterus being a junior synonym of Tilea and Scheerpeltz’s (1933) reversal to make Tilea a junior synonym of Phlaeopterus (in accordance with the principle of priority) we considered untestable using phylogenetics because neither of these concepts makes phylogenetic predictions. However, Fauvel’s (1878) erection of Tilea as a monotypic genus for P. cavicollis Fauvel, 1878 could be tested because, following current classification practices, it predicts this species would fall outside the clade containing the rest of the Phlaeopterus species. The same is true for Casey’s (1885) erection of the monotypic genus Vellica for P. longipennis Casey, 1885. However, all hypotheses relating to the monophyly and generic status of Phlaeopterus are only weakly testable given our taxon sampling. To strongly test these would require much denser sampling within Lesteva, Unamis and other anthophagine genera to test the monophyly of genera relative to all of the genera within the tribe.

All the large-bodied species are reported as snowfield-associated, and two of the small-bodied species are also (P. lagrandeuri and P. houkae Hatch, 1957). To determine if there is a relationship between body size and use of snowfields as habitats, one additional character representing habitat was added to conduct a character correlation analysis, also known as a test of independent evolution, using Pagel’s (1994) method as implemented in Mesquite 3.6 (Maddison and Maddison 2016). The habitat character was binary: (0) not snowfield associated; (1) snowfield associated, and was not used in the phylogenetic inference. The topology from the concatenated analysis was used with the option to seek any relationship (i.e. any effect of X on Y or Y on X) between character 1 (maximum body size) and this binary habitat character, with 1000 simulations to estimate the p-value. Habitat data were taken from Mullen et al. (2018) and specimen label data. This test compares the likelihoods of a four-parameter model to an eight-parameter model (we used ten iterations per simulation).

Our COI alignment is 654 bp long and comprises 164 sequences including the outgroup taxa (113 ingroup sequences). Base composition is as follows: A = 29.5%, C = 17.63%, G = 16.44%, T = 36.42%. These values are within the range typically reported for insect mitochondrial DNA (Dowton and Austin 1997). Our full COI dataset contains 216 parsimony-informative sites and 425 constant sites, with 81% (n = 175) of the parsimony informative sites found in the 3rd codon position. The Chi-square test of homogeneity of base frequencies across taxa as implemented in PAUP* 4.0 resulted in a Chi-square value = 182.0 (df = 489) with a p-value = 1.0, suggesting that nucleotide frequencies across taxa were stationary.

When P. loganensis is removed from our data, uncorrected p-distances within and among Phlaeopterus species’ COI sequences are above 4% for among-species comparisons (min 4.23%, mean 10.52%, max 14.96%) and below 4% for within-species comparisons (min 0.00%, mean 0.53%, max 3.08%). Phlaeopterus castaneus has a minimum among-species distance of 0% with P. loganensis, and is also the only species with within-species distances above 2.08% (8 of 15 within-species comparisons for this species had distances ranging from 2.29 to 3.08%).

Our most inclusive estimate of the phylogeny of Phlaeopterus resulted from the concatenated COI+morphology data (Fig. 4). The genus Phlaeopterus is well supported as monophyletic (PP = 0.98, MLBS = 96). A clade with 12 large-bodied species is well supported (PP = 0.99, MLBS = 92) but emerges from a four-branch polytomy of smaller-bodied species, including the large-bodied P. elongatus – indicating that large-bodied adults evolved twice in the genus. Within the clade of large-bodied species is nested a clade of the largest-bodied Phlaeopterus (maximum body size over 8 mm). This clade was weakly supported in the Bayesian analysis (PP = 0.72), but not supported in the maximum likelihood bootstrap analysis (MLBS < 50).

A 50% majority rule consensus tree resulting from Bayesian analysis of the COI dataset is shown in Fig. 5. This tree is better resolved than the combined COI+morphology phylogeny (Fig. 4), although the backbone has some weak branch support with posterior probabilities in the 0.70s and maximum likelihood bootstrap values < 50. This phylogeny differs from the combined COI+morphology and morphology-only phylogenies in the placement of the small-bodied P. lagrandeuri and the large-bodied P. elongatus as sister species with high Bayesian support (PP = 0.96) but weak ML bootstrap support (MLBS = 57).

This phylogeny supports the monophyly of seven of the nine Phlaeopterus species in the analysis with strong posterior probabilities (all PP > 0.98) but mixed strong to weaker maximum likelihood bootstrap support (MLBS 67–100). However, the two P. loganensis sequences nest within a clade of P. castaneus sequences. Phlaeopterus castaneus is one of two ingroup species with multiple well-supported clades separated by relatively long branches, the other being P. fusconiger.

Figure 5. 

Bayesian analysis of 164 sequences, 404–654 bp in length, of the mitochondrial gene COI with posterior probabilities and Garli maximum likelihood bootstrap values above each branch. Included are 9 Phlaeopterus ingroup species and 1 Unamis and 6 Lesteva outgroup species Topology is a Bayesian 50% majority rule consensus phylogram. One intermixed clade of the species P. castaneus + P. loganensis is indicated at (a), with (a’) indicating the subspecies P. castaneus castaneus and (a’’) indicating the subspecies P. castaneus cascadiensis. Taxon names in blue are species with maximum body sizes over 8 mm; names in red are species with maximum body sizes 6.5 to 7.9 mm with wide bodies; names in black are species with maximum body sizes below 5 mm; the name in green is for P. elongatus with a maximum body size of 6.4 mm, with an elongate body.

The 50% majority rule consensus tree of the morphology dataset (Fig. 6) contains one well supported large clade of large-bodied species that is nested within a three-branch polytomy of small-bodied species. The basal polytomy includes all the small-bodied Phlaeopterus species (maximum body size less than 5 mm, 5 species) and also includes the large-bodied P. elongatus. The larger, well supported clade (PP = 0.99, MLBS = 88) includes all the large-bodied Phlaeopterus species (maximum body size more than 5 mm, 12 species) except P. elongatus. Body lengths for all Phlaeopterus species are shown in Fig. 2.

A summary of Bayesian PP support for taxonomic hypotheses of Phlaeopterus in the COI+morphology (Fig. 4), COI-only (Fig. 5), and morphology-only (Fig. 6) phylogenies is given in Table 3. All three analyses strongly rejected Fauvel’s (1878) and Casey’s (1885) hypotheses suggesting P. cavicollis and P. longipennis were lineages separate from Phlaeopterus deserving of generic rank status (Tables 1, 3). Of the unpublished species group hypotheses, all but two that could be tested were strongly rejected as polyphyletic or paraphyletic, including the ‘Phlaeopterus sensu stricto’ hypothesis that all the large-bodied species descended from a common large-bodied ancestor (Figs 4–6).

Figure 6. 

Bayesian analysis of 40 morphological characters with posterior probabilities and Garli maximum likelihood bootstrap values near each branch with 18 Phlaeopterus species and 1 Unamis and 4 Lesteva outgroup species. Topology is a Bayesian 50% majority rule consensus phylogram. Taxon names in blue are species with maximum body sizes over 8 mm; names in red are species with maximum body sizes 6.5 to 7.9 mm with wide bodies; names in black are species with maximum body sizes below 5 mm; the name in green is for P. elongatus with a maximum body size of 6.4 mm, with an elongate body.

The stepping stone Bayes factor analyses using partial constraints, one forcing P. elongatus to belong to the clade of large-bodied species and one allowing P. elongatus to join anywhere, found strong to very strong evidence for large body size evolving twice in the genus (Table 4). Analyses that forced P. elongatus to join the clade of large-bodied species were significantly worse than those that did not. This result agrees with the strong posterior probabilities for the two branches separating P. elongatus from the clade of large-bodied species in the COI analysis (PP = 0.91 and 1.0, Fig. 5).

Ancestral character state reconstruction of character one, using the morphology only matrix and the tree in Fig. 4, under the Mk1 model, estimated the most recent common ancestor of Phlaeopterus as small-bodied with a proportional likelihood of 0.9978. Under a step matrix parsimony model for character one, the cost of evolving a large body size would have to be 2.6 times or more than the cost of evolving a small body size to estimate the common ancestor of the genus as large-bodied.

Campbell’s unpublished hypotheses for the longipennis species group, which suggested a sister species relationship between the only two species in the genus without ocelli, was strongly supported in the COI+morphology and morphology analyses (Table 3) but was untestable in the COI-only analysis in the absence of DNA sequence for P. obsoletus Mullen and Campbell, 2018. The unpublished hypothesis for the castaneus species group (P. castaneus and P. kavanaughi Mullen and Campbell, 2018, as sister species) was strongly supported in the morphology-only analysis but rejected in the other two analyses because P. loganensis was intermixed with the clade of P. castaneus sequences in the majority of sampled trees (Table 3, Fig. 5). However, the concatenated analysis did not reject this hypothesis strongly – fully one third of the post-burnin MCMC trees included P. castaneus and P. kavanaughi as sister taxa (Table 3).

The character correlation test to determine if character 1 (maximum body size) is correlated with the use of snowfields as a habitat failed to reject the null hypothesis of independence (p-value = 0.06), albeit weakly.

Differences in topology, resolution, and branch support were observed between the concatenated COI+morphology (Fig. 4), COI-only (Fig. 5), and morphology-only (Fig. 6), analyses. A strict consensus tree of all three final topologies, limited to the taxa for which we have DNA data, shares only two clades within Phlaeopterus – a polytomous clade of the wide, large-bodied species nested within the small-bodied species and P. elongatus (i.e., a clade of all Phlaeopterus species). Using the full taxon set, comparing the topologies of the concatenated and morphology-only data, three shared sister-species clades exist: P. hatchi + P. filicornis, P. frosti + P. fusconiger, and P. longipennis + P. obsoletus. These analyses also share a clade containing all Phlaeopterus species. Additionally, both datasets and all three analyses place P. castaneus and P. kavanaughi as close relatives with moderate to strong support – differing only in their relationship to P. loganensis as described in section 4.2. below. No other groupings are shared between the COI and morphology results.

The mitochondrial genome has an effective population size that is ¼ that of the somatic portions of the nuclear genome in diploid organisms, resulting in much shorter coalescence times, which should help detect closely timed speciation events but is more likely to suffer from saturation than nuclear DNA (Wiens and Penkrot 2002). Low PP and ML bootstrap values in the backbone of the COI phylogeny (Fig. 5) suggests that saturation may be confounding the phylogenetic signal of this dataset at deeper ingroup nodes. Morphological data provides a useful contrast to mtDNA, as each morphological character is coded by potentially multiple nuclear markers, and therefore may approximate their phylogenetic signal. Each type of data has strengths and limitations, and the strengths of each can complement the limitations of others. Hillis and Wiens (2000) and Wiens (2004) provide excellent summaries of the arguments for the use of both morphological and molecular data in phylogenetic analysis. In cases of conflict between morphology (Fig. 6) and COI (Fig. 5), such as the relationship of P. elongatus to P. lagrandeuri, the concatenated phylogeny (Fig. 4) recovered a topology more similar to the morphology-only phylogeny (Fig. 6), suggesting the signal in the morphological data is stronger than that in the COI data in that instance.

One Phlaeopterus species pair forms an intermixed clade in the full COI phylogeny (Fig. 5) viz. P. loganensis + P. castaneus. In an analysis of uncorrected within- and among-species p-distances these two species have uncorrected among-species p-distances of less than 2%. Specimens of P. castaneus (Table 2, UAMObs:Ento:232749) and P. loganensis (Table 2, UAMObs:Ento:232748) have 100% identical COI haplotypes. Further investigation suggests that human error is unlikely to be the explanation for these identical sequences in these two species. The two specimens were collected from different localities and sequenced by different laboratories, and the identifications of these two specimens has been confirmed by Anthony Davies of the Canadian National Collection of Insects and the second author. Morphologically, these two species are easily distinguished with the naked eye: P. loganensis has the autapomorphy of the elytra being prolonged at the suture in both males and females (fig. 29E in Mullen et al. 2018) whereas P. castaneus has broadly rounded elytral margins typical of Phlaeopterus and related taxa (fig. 29F, G in Mullen et al. 2018). Furthermore, these results are corroborated by a second pair of P. castaneus and P. loganensis sequences with less than 1% divergence. The known distribution of P. loganensis is entirely within that of P. castaneus.

Interestingly, P. castaneus forms two distinct clades in the full COI phylogeny (marked a’ and a’’ in Fig. 5) each corresponding with the subspecies recognized by Mullen et al. (2018), which corresponds with the high within-species COI distances we document in this species (2.29–3.08%). Only the inland Rocky Mountain P. castaneus castaneus clade has intermixed P. loganensis COI sequences. The other clade, P. castaneus cascadiensis Hatch, 1957, contains only two sequences, both from coastal mountain ranges, one from Haines, Alaska, and one from Mt. Garibaldi near Vancouver, British Columbia (see fig. 6B in Mullen et al. 2018 for a distribution map of P. castaneus). This pattern suggests that P. castaneus may have hybridized with P. loganensis in the Rocky Mountains (P. loganensis is only known from the Rocky Mountains and Selkirk Mountains, fig. 12A in Mullen et al. 2018) and that the Rocky Mountain P. castaneus castaneus haplotype replaced the P. loganensis haplotype. Or, alternatively, the P. loganensis haplotype replaced the P. castaneus castaneus haplotype. However, there is no morphological evidence suggesting these two species are closely related (Fig. 6), which suggests the two subspecies would not both join with P. loganensis in the COI analysis, thereby making the former hypothesis more likely. The two P. castaneus sequences from the coastal mountains correspond to the subspecies P. castaneus cascadiensis (marked a’’ in Fig. 5) these two subspecies overlap in Garibaldi and Manning Provincial Parks, British Columbia. These subspecies fall into two BINs on BOLD with BOLD:ACH1347 for P. c. cascadiensis and BOLD:AAP7088 for P. c. castaneus. Further study of these subspecies, particularly in their zone of sympatry, is warranted to rigorously test their subspecies status.

These results may also be due to incomplete lineage sorting (Maddison 1997) or infections of the bacterium Wolbachia (Whitworth et al. 2007). Studies of sequence divergence within various insect groups have reported large variation of within- and among-species distances (Cognato 2006; Trewick 2008) and high rates of non-monophyly of animal species (Funk and Omland 2003). These issues might disappear with a larger dataset including additional mitochondrial and nuclear loci, which are less likely to be introgressed. Therefore, in this case we have based our conclusions regarding species status on the morphological phylogenetic signal.

Although our survey of character states in the outgroup taxa was rather cursory and relied in part on J.M. Campbell’s assessments, these analyses suggest the following character states might be of value for further study of inter-generic relationships. Phlaeopterus and its presumed sister group, Unamis, share the following: the labrum has micropores (chr. 11), the minimum width of the gula is wide, more than 0.2 times as wide as the mentum (chr. 12), the center of the epipharynx has many closely spaced small spines (chr. 15), the fine cilia at the base of the hypopharynx are not in rows (chr. 17), the molar area of the mandible is without setae (chr. 19), the prosthecal fringe of the mandible extends beyond the subapical tooth (chr. 20), the lateral area of the pronotum is moderately to deeply impressed, foveiform (chr. 22), and the first metatarsomere is subequal to or longer than the ultimate metatarsomere (chr. 35). Unamis is autapomorphic in having a narrow labrum less than 0.3 times as long as wide (chr. 10) and the apical portion of the epipharynx is smooth and without spines (chr. 14). Additional characters diagnostic of Unamis are provided in Newton et al. (2000). Characters diagnostic of Phlaeopterus are listed in Mullen et al. (2018) and Newton et al. (2000).

The phylogeny of the tribe Anthophagini, to which Phlaeopterus belongs, is largely unknown. The tribe contains 27 North American genera, of which only six have been revised taxonomically (Campbell 1978, 1979, 1982, 1983, 1984; Mullen et al. 2018) and 43 genera globally (Kim et al. 2019). Furthermore, the tribe Anthophagini has been referred to as a taxonomic dumping ground, as it lacks synapomorphies and is likely not monophyletic (Newton et al. 2000; Kim et al. 2019). Kim et al. (2019) performed a phylogenetic analysis of the subfamily Omaliinae using three mitochondrial and three nuclear markers. Their analysis did not support the monophyly of the subfamily Omaliinae or the tribe Anthophagini, which formed two widely separated clades in their Bayesian analysis – one entirely of species of Lesteva, and the other with 11 other anthophagine genera and four coryphiine genera. No species of Phlaeopterus or Unamis were included in the taxon sampling of Kim et al. (2019). Given that a number of the deep branches in the phylogeny of Kim et al. (2019) were not strongly supported, despite their large dataset of six genetic markers, it will probably require much denser taxon sampling combined with next-generation sequencing and phylogenomic methods to confidently resolve the phylogeny of the Omaliinae.

Testing the monophyly of Phlaeopterus properly would also require much denser taxon sampling, including more Lesteva species, the addition of all or most of the known Unamis species, and representatives of most other anthophagine genera, in combination with a much larger genetic dataset. Although the large- and wide-bodied species form a well delimited clade, the delimitation between the remaining Phlaeopterus, especially the newly added Asian species, and the outgroup taxa is relatively subtle. Greater focus on inter-generic characters and relationships is warranted in future studies.

Here, we present the first phylogeny of the genus Phlaeopterus, and the first modern phylogenetic reconstruction of species-level relationships within the rove beetle subfamily Omaliinae using both morphology and molecular data. It is our hope that the coming years will see the production of much needed phylogenetic and taxonomic revisions of other anthophagine genera, resolving tribal, generic, and sub-generic relationships within the Omaliinae.