Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction?

Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction? | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction? View ORCID Profile E.K. López-Estrada , View ORCID Profile I. Sanmartín , J.E. Uribe , View ORCID Profile S. Abalde , M. García-París doi: https://doi.org/10.1101/2021.01.04.425192 E.K. López-Estrada 1 Museo Nacional de Ciencias Naturales (MNCN-CSIC) . José Gutiérrez Abascal, 2, 28006 Madrid, España 2 Real Jardín Botánico (RJB-CSIC). Plaza de Murillo , 2, 28014. Madrid, España Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for E.K. López-Estrada For correspondence: lokaren21{at}gmail.com I. Sanmartín 2 Real Jardín Botánico (RJB-CSIC). Plaza de Murillo , 2, 28014. Madrid, España Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for I. Sanmartín J.E. Uribe 1 Museo Nacional de Ciencias Naturales (MNCN-CSIC) . José Gutiérrez Abascal, 2, 28006 Madrid, España Find this author on Google Scholar Find this author on PubMed Search for this author on this site S. Abalde 1 Museo Nacional de Ciencias Naturales (MNCN-CSIC) . José Gutiérrez Abascal, 2, 28006 Madrid, España 3 Centro de Estudios Avanzados de Blanes (CEAB-CSIC). Accéss a Cala Sant Francesc , 14, 17300 Blanes, España Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for S. Abalde M. García-París 1 Museo Nacional de Ciencias Naturales (MNCN-CSIC) . José Gutiérrez Abascal, 2, 28006 Madrid, España Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Changes in life history traits, including reproductive strategies or host shifts, are often considered triggers of speciation, affecting diversification rates. Subsequently, these shifts can have dramatic effects on the evolutionary history of a lineage. In this study, we examine the consequences of changes in life history traits, in particular host-type and phoresy, within the hypermetamorphic clade of blister beetles (Meloidae). This clade exhibits a complex life cycle involving multiple metamorphoses and parasitoidism. Most tribes within the clade are bee-parasitoids, phoretic or non-phoretic, while two tribes feed on grasshopper eggs. Species richness differs greatly between bee and grasshopper specialist clades, and between phoretic and non-phoretic genera. We generated a mitogenomic phylogeny of the hypermetamorphic clade of Meloidae, including 21 newly generated complete mitogenomes. The phylogeny and estimated lineage divergence times were used to explore the association between diversification rates and changes in host specificity and phoresy, using State-Dependent Speciation and Extinction (SSE) models, while accounting for hidden factors and phylogenetic uncertainty within a Bayesian framework. The ancestor of the hypermetamorphic Meloidae was a non-phoretic bee-parasitoid, and independent transitions towards phoretic bee-parasitoidism or grasshopper specialization occurred multiple times. Bee-parasitoid lineages that are non-phoretic have significantly higher relative extinction rates and lower diversification rates than grasshopper specialists or phoretic bee-parasitoids, while no significant differences were found between the latter two strategies. This suggests that these two life strategies contributed independently to the evolutionary success of Nemognathinae and Meloinae, allowing them to escape from the evolutionary constraints imposed by their hypermetamorphic life-cycle, and that the “bee-by-crawling” strategy may be an evolutionary “dead end”. We show how SSE models can be used not only for testing diversification dependence in relation to the focal character but to identify hidden traits contributing to the diversification dynamics. The ability of blister beetles to explore new evolutionary scenarios including the development of homoplastic life strategies, are extraordinary outcomes along the evolution of a single lineage: the hypermetamorphic Meloidae. Diversification rates (speciation minus extinction) are rarely constant along the evolutionary history of a lineage ( Mooers & Heard 1997 ). Speciation and extinction rates tend to vary over time and among clades as a response to changing abiotic and biotic conditions, such as rapid climate change, geographic range fragmentation, the appearance of a new trait (“key innovation”) or an ecological opportunity resulting from the invasion of a novel niche ( Miller 1949 ; Stanley 1975 ; Sanderson & Donoghue 1996 ; Barraclough et al. 1998 ; Losos & Miles 2002 ; Ricklefs 2007 ; Maddison et al. 2007 ; FitzJohn et al. 2009 ; Stadler 2011 ; Rabosky et al. 2013 ; Donoghue & Edwards 2014 ; Donoghue & Sanderson 2015 ; Freyman & Höhna 2019 ). Shifts in the rate of diversification can have dramatic effects on the evolutionary history of a lineage, and their impact will depend on the magnitude and speed of the change, but also on the interaction with changes in other traits and with the physical template ( Vermeij 2001 ; De Queiroz 2002 ; Donoghue 2005 ; Condamine et al. 2018 ). The last decade has witnessed the advent of sophisticated methods to infer speciation and extinction rates from phylogenies of extant taxa ( Morlon 2014 ; Sanmartín & Meseguer 2016 ). Though it is mathematically possible to estimate the extinction rate from a reconstructed time tree including no fossil extinct lineages ( Stadler 2011 ), temporal and clade-specific deviations from constancy in diversification rates leads to inaccurate and often underestimated extinction rates ( Rabosky 2010 ; Morlon 2014 ). In recent years, different time-dependent and clade-dependent diversification models have been developed to overcome these issues, some assuming continuous time variation while others model change as discrete time steps ( Stadler 2011 ; Morlon et al. 2011 ; Rabosky et al. 2014 ; May et al. 2016 ; Culshaw et al. 2019 ; Höhna et al. 2019 ). Nevertheless, concerns that these models lack statistical power and can generate an infinite array of undistinguishable diversification histories remain ( Louca & Pennell 2020 , but see Morlon et al. 2020 ). Recently developed State-Dependent Speciation-Extinction models (SSE) allow researchers to test for a statistical association between the heterogeneity in diversification rates observed within a clade and the rates of evolution of a focal character that is thought to be driving diversification, e.g., a key innovation or ecological opportunity ( Maddison et al. 2007 ; FitzJohn 2010 ; Beaulieu & O’Meara 2016 ; Herrera-Alsina et al. 2019 ; May & Moore 2020 ). Models such as BiSSE (Binary State Speciation and Extinction, Maddison et al. 2007 ;) jointly estimate the transition rates among a character states and changes in speciation and extinction rates, and can therefore account for potential interactions between these processes. A perceived advantage of SSE models is that, by explicitly modeling all outcomes in the evolution of a trait with speciation and extinction, they might be able to account for extinct and unobserved lineages and therefore alleviate the problem of undistinguishable diversification histories ( Höhna et al. 2019 , but see Louca & Pennell 2020 ). There are, however, warnings about the risk of overconfidence in SEE models, especially the problem of “pseudoreplication” ( Maddison & FitzJohn 2015 ) and inflated Type I error, i.e. an association is detected when there is none. If the phylogeny exhibits high heterogeneity in diversification rates among clades, SSE models may lead to mistaken inference of state-dependent diversification ( Rabosky & Goldberg 2015 ). This issue has spurred the development of Hidden State-dependent Speciation and Extinction models (HiSSE, Beaulieu & O’Meara 2016 ), allowing for the detected variation in diversification rates to be unrelated to the focal character but instead explained by an unobserved character trait. Typically, the HiSSE model is used to corroborate or reject the association between the rate of evolution of the focal trait and changes in the rate of diversification detected by “standard” SSE models ( Condamine et al. 2018 ; Fernandez et al. 2018 ; Gajdzik et al. 2019 ; Nakov et al. 2019 ). However, these models could also be used to identify the “hidden trait” or unknown causal force whose interaction with the observed trait is behind the heterogeneity in diversification rates; to our knowledge, this other use has never been explored. A prerequisite for estimating extinction and speciation rates from reconstructed time trees is a sound, well-supported phylogeny. Evolutionary radiations are difficult to tackle with statistical phylogenetic methods because they often exhibit short internal branches, leading to stochastic error, and heterogeneity in gene-tree topologies (due, among others, to incomplete lineage sorting and gene flow), leading to systematic error ( Philippe et al. 2011 ; Cai et al. 2020 ). In the last decade, the use of genomic approaches has revolutionized phylogenetics, with the possibility to solve radiations in non-model organisms (Barret et al. 2016; Villaverde et al. 2018 ; Allio et al. 2020 ; Young & Gillung 2020 ; but see Cai et al. 2020 ). In bilaterian animals, including insects, shotgun sequencing of whole mitochondrial genomes has become a favorite tool for reconstructing robust phylogenetic hypotheses at an affordable cost; this is partly due to the small size of the mitochondrial genome, the presence of multiple copies in the cell, and its haploid non-recombinant nature ( Finstermeier et al. 2013 ; Gibb et al. 2016 ; Yuan et al. 2016 ; Liu et al. 2019 ; Yan et al. 2019 ; Irisarri et al., 2020 ). Changes in life history traits, including development of new reproductive strategies and host shifts, are often considered powerful triggers of speciation bursts ( Bonett & Chippindale 2004 ; Hardy & Otto 2014 ). Host shifts can produce a major turnaround in the evolutionary fate of a parasitic lineage, and might result in a “wild” increase in species number or in new levels of biological complexity ( Erwin 1992 ; Ricklefs & Fallon 2002 ; Silva et al. 2012 ). Though host specialization is common in parasites (i.e. the “one-parasite-one host” rule), there is an increasing body of literature showing that changes in host specificity are not infrequent ( Hardy & Otto 2014 ; Nylin et al. 2018 ; Braga et al. 2020 ). By jumping hosts, parasites can escape extinction and increase their probability to persist over long evolutionary times ( Thines 2019 ) or over ecological time scales ( Brooks et al. 2006 ; Calatayud et al. 2016 ). Host shifts have been widely documented among organisms, especially in humans, where the majority of pathogens originate through host changes, including HIV, malaria, and the most recent SARS-CoV-2 ( Wolfe et al. 2007 ; Zhang et al. 2020 ). Many factors intervene in the evolutionary success of a host shift, including physiological similarity between parasite and host ( Runge & Thines 2012 ; Thines 2019 ); phylogenetic, ecological and geographical distance between current and potential hosts ( Göker et al. 2004 ; Engelstädter & Fortuna 2019 ); or differences in parasite and host mutation rates – i.e. higher mutation rates may allow the parasite to overcome host defensive responses – ( Gandon & Michalakis 2002 ). In other words, the likelihood of a new host-parasite interaction and the specific evolutionary outcome of the new relationship are often not determined by a change in a single character (a “key innovation”) but by the interaction of multiple causal agents or changes in different character traits. Changes in different character traits can either appear “simultaneously” during a single speciation event, or as additive factors acting synergistically across several nested speciation events (a “synnovation”, sensu Donoghue & Sanderson 2015 ). Trait evolution such as host specialization or the development of intricate life-history strategies is deeply affected by the effect of “historical contingency”, defined as “the series of historical details that lead to a specific evolutionary outcome” ( Gould 1989 ). The link between contingency and evolutionary radiation has attracted considerable attention ( Losos et al. 1998 ; Vermeij 2001 ; De Queiroz 2002 ; Donoghue 2005 ; Beatty 2006 ; Blount et al. 2008 ; Givnish et al. 2014 ; Kriebel et al. 2020 ). According to Beatty (2006) , Gould’s contingency” offers two interpretations, which are not mutually exclusive but complementary: a) causal-dependence: a change on a trait is “contingent upon” an earlier change in a different trait, either promoting ( De Queiroz 2002 ; Losos et al. 1998 ), or restraining it ( Vermeij 2001 ); b) unpredictability: previous states of the character trait are necessary but insufficient to lead to the resulting outcome. So far, contingency in host specialization has never been explored using SSE models. Yet, because host shifts often imply a major change in ecological interactions ( Calatayud et al. 2016 ), they are probably driven by multiple interacting causal factors, which may act in confluence (“synnovation”) or contingent upon one another. Another interesting consequence of causal-dependence in the context of host-parasite associations is the possibility of becoming an “evolutionary dead-end” ( Vamosi et al. 2003 ). This process is generally associated to the acquisition of a character state with high extinction or low speciation rates, or with irreversibility in transition rates (Goldberg & Igic 2008; Goldberg et al. 2017 ). Host specialization has sometimes been considered an evolutionary dead-end because the acquisition of a narrow set of food resources (hosts), and concomitant trait adaptations, may limit future diversification, in contrast with phenotypic plasticity in generalist species ( Hardy & Otto 2014 ). Host-jumps that require changes in multiple levels of biological complexity (e.g., morphological, anatomical, physiological, and ecological) are especially difficult to reverse. Commonly known as “blister beetles”, Meloidae is a family of Coleoptera (Tenebrionoidea), which includes circa 3000 described species ( Bologna et al. 2008 ). Its vernacular name is related with the capacity to synthetize cantharidin, a sesquiterpenoid toxin ( Percino-Daniel et al. 2013 ; Bravo et al 2017 ). This substance is a powerful systemic poison, acting mainly in tissue degradation (lethal dose in humans: 0.5 mg/kg), and a high deterrent for invertebrate and vertebrate predation ( Kaiser & Michl 1958 ; Carrel & Eisner 1974 ; Bertaux et al. 1988 ). Cantharidin has also been associated to parasitic regulation in birds ( Bravo et al. 2014 ). Because of the pharmacological usage of cantharidin, blister beetles have drawn scientific attention since the origins of Zoology ( Dioscorides 1636 ; Fischer 1827 ; Amor Mayor 1860 ). Meloidae includes three lineages: Eleticinae, with approximately 100 species and a Gondwana-like distribution (absent from Australia); Nemognathinae, comprising ∼520 species and distributed in all continents except New Zealand and Antarctica; and Meloinae, the richest subfamily with almost 2500 species, sharing its geographic distribution with Nemognathinae ( Pinto & Bologna 1999 ; Bologna 2009 ; Bologna & Pinto 2002 ). In the most recent phylogeny of Meloidae ( Bologna et al. 2008 ), Eleticinae was reconstructed as sister to the clade formed by Meloinae and Nemognathinae. Eleticinae exhibits the non-parasitic life cycle typical of most Tenebrionoidea, which is characterized by a single event of metamorphosis ( Pinto et al. 1996 ; Bologna & Di Giulio 2011 ). In contrast, Meloinae and Nemognathinae exhibit a unique, intricate “hypermetamorphic” life cycle. This type of unusual development usually involves three metamorphoses before the imago, with at least four larval phases, each of them so different from each other in morphology and behavior that they could be easily assigned to other families different from Meloidae ( Fig. 1 ; Bologna & Pinto 2001 ; Bologna et al. 2008 ). Download figure Open in new tab Figure 1. Life cycle of hypermetamorphic Meloidae parasitizing sub-social bees. Eggs are laid in the ground (Meloinae) or on the phyllaries of flowers (Nemognathinae) mainly Asteraceae ( Enns 1956 ). When the eggs hatch, a highly mobile larva emerges from the hatching site and searches for bee nests. This first larva (“triungulin”) can either, depending on the lineage, climb a flower and wait for a bee to visit the flower and then attach to it and be transported to the nest (phoresy), or wander around the ground until it finds an entrance to the bees’ nest (active searching). Once the first larva reaches the bees’ nest, it starts eating the provisions, eggs or larvae of usually a single cell. A first metamorphosis then occurs, and the second larva known as “first grub larva” emerges; it presents reduced motility, but feeds nearly continuously until its metamorphosis. In the second metamorphosis, the first grub larva changes into a “coarctate larva”; the larva loses its appendages and enters into diapause. A third metamorphosis occurs, and a second grub larva morphologically similar to the second larva develops; this larva recovers the motility, although it does not feed. Finally, the larva pupates, and in a few days the adult emerges and the cycle begins again. Though the hypermetamorphic life cycle is present in all species of Nemognathinae and Meloinae, two traits exhibit variation at the tribal and generic levels: the mode of locomotion used by the first instar larva to reach the food source ( Fig. 1 ), and the host itself ( Fig. 2 ). The four tribes within Nemognathinae and six of the eight tribes included in Meloinae are parasitoids of many different species of solitary or subsocial bees (superfamily Apoidea) feeding on all resources available at the nest: eggs, bee larvae, and provisions ( Fig. 1 ). The two exceptions are Epicautini and Mylabrini, which first instar larvae feed on eggs within the egg-pods of grasshoppers of the family Acrididae ( Fig. 2 ; Bologna 1991 ; but see Bologna & Di Giulio 2011 , for a possible case of parasitizing Sphecidae). Bologna et al. (2008) proposed that Epicautini and Mylabrini were not sister-taxa and as a consequence each event of host-jump to Acrididae was independent from each other (homoplastic). The mode of locomotion of first instar larvae to reach the food source is another trait that varies across taxa; this includes phoresy, i.e. passively latching onto a bee to reach the nest, and active crawling. In phoretic taxa, first instar larvae climb to flowers and attach to passing bees ( Fig. 1 ; Hafernik & Saul-Gershenz 2000 ); in non-phoretic species, larvae wander on the ground, actively searching for bee nests ( Fig. 1 ) or grasshoppers’ egg-pods ( Fig. 2 ). All species of Nemonagthinae are phoretic parasitoids of bees (with the probable exception of Stenodera ; Bologna et al. 2002 ; Bologna & Di Giulio 2011 ), whereas some tribes of Meloinae that are bee-parasitoids (e.g. Meloini) include phoretic ( Meloe ) and non-phoretic ( Physomeloe ) genera. Phoresy is hypothesized to have evolved at least two times independently in bee-parasitoid Meloinae ( Bologna & Pinto 2001 ; Bologna et al. 2008 ). Download figure Open in new tab Figure 2. Life cycle of hypermetamorphic Meloidae specialized in grasshopper eggs. Meloid eggs are laid in the ground. A highly mobile larva emerges from the hatching site and actively searches for grasshoppers’ egg-pods. Once the first larva reaches the pod, it starts eating the eggs. A first metamorphosis then occurs, and the second larva known as “first grub fase” emerges. The following larval phases are as those of the lineages that parasitize bees’ nests. Jumps and reversions between larval stages have been observed mainly in Epicautini and Mylabrini ( Selander & Mathieu 1964 ; Selander & Weddle 1969 ). Bologna et al. (2008) argued that host specificity could explain the remarkable difference in species richness observed between Nemognathinae (∼500 species), and Meloinae, with ∼2500 species. It may also explain differences among tribes of Meloinae: the two tribes feeding on grasshopper eggs, Mylabrini and Epicautini, include circa 600-700 species each, while the largest tribes that parasitize bees (Meloini and Pyrotini) do not exceed 300 species together ( Table 1 ). There is no obvious pattern of species richness across phoretic and non-phoretic genera ( Table 1 ). For example, in the tribe Lyttini, the non-phoretic genus Lytta comprises 109 species, while Lagorina only includes two. Some phoretic genera of Meloini like Meloe are species-rich (153 species), while others ( Spastonyx, Lyttomeloe … ) comprise but a few species ( Table 1 ). Bologna & Pinto (2001) suggested that phoresy could be evolutionarily advantageous in Meloinae, because it allows long-distance dispersal of the first instar larvae, thus potentially increasing the geographic ranges of species; and, second, because transporting the larvae to their food source ensures successful larval rearing. However, these authors also questioned whether this mode of locomotion is more effective in finding the host than the non-phoretic, active crawling behavior present in other genera ( Bologna & Pinto 2001 ). For example, the grasshopper specialist tribes, Mylabrini and Epicautini, are both species-rich and non-phoretic. View this table: View inline View popup Download powerpoint Table 1. Taxa included in this study, with associated taxonomic diversity, presence of a phoretic behavior, and host-type specialization. Species richness were obtained from Pinto & Bologna (2009) , Bologna & Pinto (2002) and Campos-Soldini et al. (2018) . In this study, we reconstruct phylogenetic relationships and estimate lineage divergence times in the subfamilies Nemognathinae and Meloinae (with special reference to the latter), and use the resulting time tree and macroevolutionary statistical likelihood methods within a Bayesian framework to examine the role played by host-type and phoresy as drivers of diversification in the hypermetamorphic clade of Meloidae. To date, hypotheses of relationships among major clades of Meloinae and Nemognathinae have been proposed based on morphological characters ( Denier 1935 ; MacSwain 1956 ; Selander 1964 , Kaszab 1969 ; Bologna & Pinto 2001 ), or with combined morphological and molecular ( 16S and ITS2 ) datasets ( Bologna et al. 2008 ). These phylogenies, however, exhibit low statistical clade support, and were constructed using traits directly related to the problems we are interested in, i.e. morphological traits associated to phoresy. To provide an independent robust phylogenetic hypothesis, we sequenced whole mitochondrial genomes for 15 genera of Meloinae, covering 75% representation of the tribal diversity and 20% generic diversity ( Table 1 ). Additionally, we sequenced five genera of Nemognathinae, representing 50% of tribal diversity and 10% generic diversity ( Table 1 ). In total, we generated 21 new mitogenomes by shotgun sequencing. To account for incomplete taxon sampling we developed an integrative approach incorporating the mitogenomic phylogeny, clade species richness, divergence times, and birth-death simulations, to mimic a larger, clade-representative, taxon sampling. We then used state-dependent speciation and extinction (SSE) models that tie changes in diversification rates to transitions between states in a focal character, while accounting for potential interactions with unobserved (hidden) traits ( Maddison et al. 2007 ; FitzJohn 2012 ; Beaulieu & O’Meara 2016 ) in a Bayesian framework. Specifically, with this approach we wanted to answer the following questions: (i) are shifts in host-type responsible for the observed difference in species richness among tribes and between subfamilies? (ii) What is the role of phoresy in diversification dynamics? We demonstrate that the heterogeneity in diversification rates observed in the hypermetamorphic clade of Meloidae was the product of a life strategy change, involving either a move towards phoresy or a change of host-type. We further discuss how contingency may have shaped these life strategy changes and how changes in the diversification processes offered new pathways to escape from possible evolutionary dead-ends. Additionally, we show how SSE models can be used to identify an unobserved character, whose interaction with the focal trait is shaping diversification dynamics. MATERIALS AND METHODS Taxon Sampling, DNA Extractions, Sequencing and Assembly Our ingroup dataset was composed of 29 mitogenomes of Meloidae (Tenebrionoidea). It included 23 species from 15 genera of Meloinae, representing six out of the eight tribes of the subfamily (Epicautini, Eupomphini, Lyttini, Meloini, Mylabrini, and Pyrotini; representatives of Cerocomini and Tetraonycini could not be included, Table 1 ), and six species from five genera of the sister subfamily Nemognathinae, representing two of the four recognized tribes (Nemognathini and Horiini; representatives of Palaestrini and Stenoderini were missing; Table 1 ) ( Fig. 3 ). From this dataset, mitogenomes were newly generated for 21 specimens (Table S1). These specimens were collected in the field, stored in 100% ethanol and deposited at the Museo Nacional de Ciencias Naturales (MNCN-CSIC, Madrid, Spain). Mitogenomes of the remaining eight species were obtained from GenBank. To root the phylogeny, we selected from GenBank ten additional mitogenomes from six families. According to the most recent phylogenetic hypotheses in Coleoptera ( Timmermans et al. 2015 ; Yuan et al. 2016 ), these ten taxa represent a wide spectrum of the superfamily Tenebrionoidea. A complete list of the specimens included in this study and GenBank accession numbers is shown in Table S1. Download figure Open in new tab Figure 3. Habitus of representative species of Meloidae used for this study. Live adult specimens of: A) Apalus guerini (Mulsant 1858), a phoretic bee specialist (from Perales de Tajuña, Spain) (Nemognathinae: Nemognathini); B) Pyrota palpalis Champion 1893, a non-phoretic bee specialist (from Lordsburg, New Mexico) (Meloinae: Pyrotini); C) Oenas fusicornis Abeille de Perrin 1880, a non-phoretic bee specialist (from Tielmes, Spain) (Meloinae: Lyttini); D) Epicauta tenella (LeConte 1858), a grasshopper specialist (from Needles, California) (Meloinae: Epicautini); E) Physomeloe corallifer (Germar 1818), a non-phoretic bee specialist (from Serranillos, Spain) (Meloinae: Meloini); F) Croscherichia paykulli (Billberg 1813), a grasshopper specialist (from Moulay Bousselham, Morocco) (Meloinae: Mylabrini). Photographs by MGP. Tissue samples were obtained from the thoracic muscle of the hind coxae. Total genomic DNA was extracted using the “BioSprint 15 DNA Boold Kit” (Quiagen ®) following the protocol described by the manufacturer. A pair of primers (forward and reverse) were designed in the ribosomal gene 16S [Mel-16SF and Mel-16SR] (Appendix S1). These two primers, that virtually match in all the 16S fragments studied in Meloidae, were combined with a set of cox1 pairs of primers (forward and reverse) designed for each tribe. The 16S - cox1 primers were used to amplify by long polymerase chain reactions (PCR) the mitochondrial genomes in two overlapped amplicons (with ∼13 kb and ∼5 kb long; the amplification strategy and the sequence of each primer is described in Appendix S1). Amplifications were conducted following the protocol of Uribe et al. (2017) . For each mitogenome, two fragments were amplified using the following PCR conditions: denaturing step at 94°C for 60 s; 45 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and extension at 68°C for 60 s per kb; final extension step at 68°C for 12 min. PCR products were cleaned following an ethanol precipitation protocol. Sequencing and assembly schemes were carried out following the protocol described by Abalde et al. (2017) . Briefly, an indexed library per sample was constructed using the NEXTERA XT DNA library prep kit (Illumina, San Diego, CA, USA) and sequenced in an Illumina MiSeq platform at Sistemas Genómicos (Valencia, Spain). Mitochondrial genomes were assembled mapping the raw reads back to the 16S sequences previously Sanger-sequenced, using Geneious ® 11.0.5. We set for the first two iterations a minimum overlap of 60% and a minimum match identity of 99% in order to account for potential sequencing errors. These parameters were increased in the following iterations to 75% and 100%, respectively, until completion of the mitochondrial genomes. Number of reads and mean coverage for each genome are presented in Table S1. Annotation of tRNA genes was carried out through the MITOS2 Web Server ( Bernt et al. 2013 ) and tRNAscan-SE Web Server ( Chan & Lowe 2019 ), which infers cloverleaf secondary structures. The protein-coding genes were checked by aligning them against related mitogenomes to confirm the correct annotation of the start and stop codons. Ribosomal genes were identified and annotated by similarity with the homologous genes of already published mitochondrial genomes of blister beetles and assumed to extend to the boundaries of adjacent genes ( Boore et al. 2005 ). Phylogenetic Inference and Divergence Time Estimation Phylogenetic reconstruction was performed using only protein-coding and ribosomal RNA genes extracted from the complete mitogenomes ( Abalde et al. 2017 ). Sequences of protein-coding genes were extracted into single-gene matrices, subsequently aligned based on the corresponding amino acid translations using TranslatorX Web Server ( Abascal et al. 2010 ), with the MAFFT algorithm ( Katoh et al. 2005 ). Amino acid and nucleotide raw alignments were trimmed with Gblocks ( Castresana 2000 ) using the following specifications: excluding many contiguous non-conserved positions and allowing gap positions within the final blocks. Ribosomal genes were aligned and cleaned through MAFFT and Gblocks online services ( Katoh et al. 2017 ; Talavera & Castresana 2007 ). We constructed two datasets using the single-gene matrices: (a) NT-matrix, including all nucleotide sequences (12095 bp) and (b) AA+rNT-matrix, including DNA nucleotide sequences for the non-coding ribosomal genes and the translated amino acid sequences for the coding genes (5170 sites). Maximum Likelihood (ML) phylogenetic analyses were run on the NT-matrix, using the software RaxML v7.3.1 ( Stamatakis 2006 ) and the “GTRGAMMA” option, which implements the GTR model ( Tavaré 1986 ), with four gamma-distributed rate categories ( Yang 1994 ) for all partitions. ML analyses were conducted with default parameters using the rapid hill-climbing algorithm and 1000 bootstrap pseudo-replicates. Bayesian Inference (BI) analyses were performed using MrBayes v.3.2.6 ( Ronquist et al. 2012 ). PartitionFinder v2 ( Lanfear et al. 2016 ) was used to select the best partition scheme and molecular evolutionary models for the NT and AA+rNT matrices, under the Bayesian Information Criterion (BIC; Schwarz 1978 ). MrBayes analyses consisted on two simultaneous chains of 100 million generations each, sampling trees every 10000th generation. Convergence and mixing among chains were evaluated by checking the average standard deviation of split frequencies ( 200). A majority consensus tree for each analysis was reconstructed after discarding the first 25% of the sampled trees as burn-in. Differences in biochemical profiles across sites have been shown to be an important source of systematic error in deep-time phylogenomics ( Philippe et al. 2011 ). We used the Bayesian software PhyloBayes ( Lartillot et al. 2009 ) to analyze the NT and the AA (only) data matrices under the site-heterogeneous CAT molecular model: this model implements an infinite mixture model, in which sites are allowed to have different stationary frequencies or biochemical profiles ( Lartillot & Philippe 2004 ). We conducted analyses under the CAT-GTR model, which allows rate variation across sites under the GTR model ( Tavaré 1986 ), and the simpler CAT-Poisson model, assuming a single transition rate for all sequence positions. Analyses were run including both variable and constant sites, since excluding the latter has been shown to mislead phylogenomic inference (i.e. Thode et al. 2020 ). For every analysis, we ran 4 independent chains of one million generations, sampling every 10th state. After 20000 states, we checked every 5000th generation the convergence of chains. We used the tracecomp command to evaluate differences in effective samples sizes (EES) estimated by each chain, with a minimum value of 50 and using a burn-in of 10000 tree samples (100000 generations). Differences in bipartition frequencies between topologies were evaluated with the bpcomp and the readpb commands; results were summarized using a 10% burn-in, and a sub-sampling of every 10th trees from the posterior distribution. A value equal or lower than 0.1 was used as the maximum accepted difference ( maxdiff ) between chains as a measure of successful convergence, following the PhyloBayes manual ( http://www.phylo.org/tools/pb_mpiManual1.4.pdf ). We tested the hypothesis that host shift (bees to grasshoppers) occurred independently in Mylabrini and Epicautini ( Bologna et al. 2008 ). We specified two hypotheses of relationships between Mylabrini and Epicautini to be tested against each other; one involving monophyly of the grasshopper-specialists (e.g. Mylabrini and Epicautini are sister taxa, and consequently host shift occurred once in the history of hypermetamorphic Meloidae) ( H 0 ), and the alternative, in which Mylabrini and Epicautini do not for a monophyletic group ( H 1 ). We used Bayes Factor to compare the marginal likelihood of two models representing each one of the hypotheses obtained applying the topological constraint function in MrBayes to the Bayesian topology ( Fig. 4 ). One model was constructed by constraining the topology reconstructing Epicautini and Mylabrini as a monophyletic group ( H 0 ). For the alternative model ( H 1 ) we used the unconstrained topology ( H 1a ) ( Fig. 4 ), and also a constrained topology in which Mylabrini was sister to all other Meloinae (including Epicautini) ( H 1b ). Download figure Open in new tab Figure 4. Mitogenomic phylogeny of hypermetamorphic Meloidae. Ultrametric tree obtained using BEAST based on the concatenated mitochondrial dataset. The translated amino acid sequences for the coding genes + the non-coding ribosomal genes matrix (AA+rNT-matrix) was used for this analysis. The chronogram shows mean ages for lineage divergences times using Bayesian relaxed clocks; gray circles near nodes indicate a posterior probability (PP)>0.95 and bootstrap values above 70; purple horizontal bars show 95% HPD values; color shades represent different tribes; red color branches represent phoretic; the characteristic host of each tribe is represented next to its clade. To estimate the marginal likelihood of each model, we computed power posteriors using the stepping-stone sampling method of Xie et al. (2011) implemented in MrBayes, with the default values for the α–shape parameter of the beta distribution (0.4) and a burn-in of 250000. A MCMC chain of 250000 MCMC steps was sampled every 2500 generations for each of 50 power posteriors, with ß-values ranging between 1 (posterior) and 0 (prior). We compared the two marginal likelihood values using the likelihood ratio test, 2Ln ( H 0 – H 1 ), with values larger than 2 indicating positive support for one model over the other, and larger than 6, indicating strong positive support ( Kass & Raftery 1995 ). Since the hypothesis of monophyly ( H 0 ) was rejected, we did not consider necessary adding further constraints to the alternative topology ( H 1a, b ), because according to Bergsten et al, (2013) the standard Bayes Factor test is biased toward acceptance of the hypothesis of monophyly. Divergence times were estimated using Bayesian relaxed clocks implemented in BEAST v1.8.2 ( Drummond et al. 2012 ). Analyses were performed using the concatenated AA+rNT dataset. We used for the coding genes the WAG model of amino acid substitution ( Whelan & Goldman 2001 ), while for ribosomal genes we used the HKY model of molecular substitution ( Hasegawa et al. 1985 ). A fossil specimen from the Dominican amber (Cordillera Septentrional, Dominican Republic) was identified and described by Poinar (2009) as a larva of Meloe dominicanus . This fossil represents the oldest record for the genus Meloe and was used to calibrate the node that clusters the two samples of this genus included in our data set ( M. cavensis and M. mediterraneus ). Since the age of Dominican amber is controversial, ranging from 15 to 45 Ma (Cêpek in Schlee 1990 ; Iturralde-Vincent & MacPhee 1996), we used a log-normal prior distribution with offset=15, mean=10, and standard deviation=1.25 to span the proposed time window for the Dominican amber. Molecular clocks were unlinked across genes, and modeled through an uninformative prior (gamma distribution, initial value=0.01, shape=0.01, offset=0). We used the Birth-Death model with incomplete taxon sampling ( Stadler 2009 ) as the tree growth prior, to account for missing taxa in our dataset. The analysis was run for 100 million generations, sampling every 1000th. We inspected the trace plots and effective sample sizes in Tracer 1.8.0 ( Drummond & Rambaut 2007 ), after discarding the first 20 million generations as burn-in, to assess appropriate convergence and mixing of the chain. All analyses were run in the web public resource CIPRES Science Gateway version 3.3 ( Miller et al. 2010 ). State-Dependent Diversification and Incomplete Taxon Sampling To test for a causal correlation between observed differences in diversity and species variation in traits such as host specificity or mode of locomotion in Meloinae, we used State-dependent Speciation and Extinction (SSE) models ( Maddison et al. 2007 ; FitzJohn 2012 ; Beaulieu & O’Meara 2016 ). Incomplete taxon sampling has been shown to severely bias the estimation of diversification rates, especially for the relative extinction parameter (extinction / speciation), under both rate-constant and rate-variable birth-death models ( Stadler 2009 ; Höhna 2014 ; Morlon 2014 ; Sanmartín & Meseguer 2016 ; Louca & Pennell 2020 ). Our phylogeny represents less than 1% of the extant species in the family Meloidae. SSE diversification models can potentially account for incomplete taxon sampling by including a parameter ρ that incorporates the global sampling fraction, i.e. the proportion of taxa included in the phylogeny relative to the extant diversity ( Höhna 2013 ; Beaulieu & O’Meara 2016 ). However, this parameter assumes that missing species are uniformly distributed across clades, which is often not the case ( Höhna 2014 ). Ignoring the fact that some clades are better represented than others in the phylogeny, can potentially bias estimates of extinction rates, since the probability that a lineage in state i at time t would go completely extinct by the present is 1 minus the sampling fraction, 1– ρ ( Maddison et al. 2007 ; Louca & Pennell 2020 ). One option to deal with this problem is to distribute incomplete taxon sampling values across clades in the phylogeny, i.e. including clade-specific sampling fractions ( Rabosky et al. 2014 ). However, such a procedure may lead to incorrect estimation of extinction rates in SSE models: artifactual rate shifts in extinction will be inferred wherever any two clades that differ in their sampling fractions coalesce, rendering the likelihood estimation incorrect ( Moore et al. 2016 ; Beaulieu 2020 ). To avoid this bias and mimic a larger, clade-representative taxon sampling, which allow us to use SSE diversification models, we developed an integrative approach that incorporates taxonomic information on species richness, the backbone phylogeny, the divergence times previously estimated, and the use of birth-death simulations. First, we calculated the relative extinction rate under a constant-rate diversification model for the original 29-tip phylogeny, using the bd.shifts.optim function implemented in the R package TreePar ( Stadler 2011 ), with the global sampling fraction set to 0.01%. Second, for several major clades in the backbone phylogeny (tribes, species-rich genera) (Fig. S1), we calculated the net rate of diversification using the method-of-moments (MoM) estimator ( Magallón & Sanderson 2001 ), implemented in the function bd.ms in the R package geiger ( Harmon et al. 2008 ). This function allows estimating the net diversification rate conditioning it to the crown-age of the clade (retrieved from the BEAST MCC tree generated above), the number of extant species (obtained from Pinto & Bologna 1999 ; Bologna & Pinto 2002 ; Campos-Soldini et al. 2018 ), and the relative extinction rate; the latter was fixed to the value estimated by TreePar for the backbone phylogeny. Third, for each of the major clades mentioned above (Fig. S1), we simulated 100 subtrees under a constant-rate model with the function tess.sim.taxa.age included in the R package TESS ( Höhna 2013 ). Simulations used the MoM estimators of diversification rates for the clade and the TreePar turnover value, and were conditioned to clades’ crown age and number of extant species. Fourth, we randomly selected one of the 100 subtrees and bound it to the most inclusive node subtending the clade in the backbone phylogeny, using the bind.tree function in the R package ape ( Paradis et al. 2004 ). We only simulated subtrees for clades/nodes for which the time of origin was estimated using two species at least. For example, we simulated a subtree with 26 species of Eupomphini ( López-Estrada et al. 2019 ), which was bound to the node connecting the two representative genera Megetra and Tegrodera in our phylogen y ; for Nemognathinae, we simulated a 523 species subtree and bound it to the node which was the MRCA of the six species of this subfamily in our backbone phylogeny. Tribes represented by only one tip could not be simulated (e.g. Pyrotini). For some tribes, we simulated subtrees at the generic level, as well, to represent the most speciose genera, for example Hycleus in tribe Mylabrini. Figure S1 indicates those nodes in the backbone tree, for which we generated and bound subtrees with associated taxonomic richness. This procedure allowed us to construct a phylogeny with 2150 tips out of approximately 3000 species described for Meloidae, representing 71% of species diversity. It is important to note that these simulated subtrees have no effect on the representation of character states in the phylogeny, as all simulated clades were homogeneous for the trait in question, i.e. each subtree has exactly the same state for all simulated tips. Our aim here was to incorporate the observed differences in species diversity among clades, while decreasing the potential bias introduced by low and uneven sampling when estimating of trait-dependent diversification rates (i.e. subtrees were simulated under clade-specific diversification rates). To test whether host jump has been a driver of diversification in Meloinae, we used the Binary State-dependent Speciation-Extinction (BiSSE) model ( Maddison et al. 2007 ), allowing for speciation and extinction rates to differ between two states within a character ( Maddison et al. 2007 ; FitzJohn 2010 ; 2012 ; Beaulieu & O’Meara 2016 ; Sanmartín & Meseguer 2016 ). We used the Bayesian MCMC implementation of this model in the open software RevBayes ( Höhna et al. 2016 ). The hierarchical Bayesian approach to SSE models in RevBayes, represented as directional acyclic graphs (DAGs, Höhna et al. 2014 ), allows estimation of the marginal posterior probabilities for transition rates and ancestral states through the use of hyperpriors, and is therefore efficient in integrating uncertainty in parameter values ( Freyman & Höhna 2019 ). We coded terminals for two states: parasitoids of bees (0) and grasshopper specialists (1). We followed similar priors to Freyman & Höhna (2019) . Speciation and extinction rates for the two states were modeled with a log-normal prior with an expectation centered in the total extant diversity under a constant rate diversification model (rate_mean= ln(Total number of species/2.0)/crown root age); the log-standard deviation was set to a value of = 0.587405 (rate_sd= 0.587405), which assigns a 95% credibility interval that spans one order of magnitude around the mean. Two indirect parameters were also estimated, diversification rate, speciation minus extinction, and relative extinction rate, the ratio of extinction to speciation. Transitions rates were drawn from an exponential distribution with a mean of ten character state transitions over the tree length. Root state probabilities were modeled with a Dirichlet prior with mean= 1. We set the global sampling fraction to 0.71 to account for incomplete taxon sampling in our empirical backbone + simulated subtrees phylogeny. The analysis was run with a chain length of 40000 generations. We summarized ancestral states as nodal marginal posterior probabilities using the maximum a posteriori tree and code provided in RevBayes ( https://revbayes.github.io/tutorials/morph/morph_more.html ). Finally, we employed a heuristic approximation to stochastic character mapping that does not require a rejection-sampling step ( Freyman & Höhna 2019 ), to estimate the number and timing of transition events between states. Stochastic character mapping was run in RevBayes with 500 time slices. The script to run this analysis is provided in Appendix S2-1. The frequency of type I error, detecting a signal of trait-dependency when there is not, has prompted some authors to criticize SSE models because of the simplistic, unrealistic null model used for comparison, a lineage evolving with trait-independent but constant diversification rate ( Rabosky & Goldberg 2015 ; Beaulieu & O’Meara 2016 ). This null model is regarded as a “straw man” because heterogeneity in net diversification rates is to be expected in any phylogeny, except when dealing with very recent times, albeit not necessarily correlated with the focal trait ( Rabosky & Goldberg 2015 ). This issue spurred the development of “hidden-trait” SSE models, such as the Hidden State Speciation Extinction (HiSSE) model ( Beaulieu & O’Meara 2016 ), which allows for the detected variation in diversification rates to be unrelated to the focal character but instead explained by other unobserved character. The hidden, uninformative trait interacts with the character of interest and acts as a “stand-in” for the unknown evolutionary factor that is leading the observed heterogeneity in diversification rates ( Beaulieu & O’Meara 2016 ). Our HiSSE model contained two hidden states (A, B) within each of the two observed states of “host-type”; in total, four character states: parasitoids of bees (0) with hidden states A and B (resulting in 0A and 0B), and grasshopper specialists (1) with hidden states A and B (1A and 1B). If the differences in speciation and extinction rates detected by BiSSE for the observed character states 0 and 1 holds for the two hidden states (A and B), we can conclude that the focal character drives diversification rate heterogeneity among clades. When this is not the case, the observed differences in species richness among clades cannot be associated to the focal character, but to a hidden trait to which diversification is related ( Freyman and Höhna 2019 ). Initial runs with the same priors as BiSSE resulted in poor mixing. We used instead alternative priors suggested in the RevBayes tutorial ( https://revbayes.github.io/tutorials/sse/hisse.html ). Speciation and extinction rates for the observed states were modeled as an identical broad log-uniform distribution with bounds “1E-6 and 1E2”; hidden states were modeled through a discretized log-normal distribution, with the number of categories (quantiles) equal to the number of hidden states; the mean of the log-normal for each hidden rate was set to 1.0, to make them relative, and the standard deviation set to an exponential distribution with mean= 0.587405. Transitions rates between observed states were allowed to differ and modeled through an exponential distribution prior with lambda 10, as in the BiSSE model; transition rates between hidden sates were constrained to be equal. All other prior settings followed BiSSE. The analysis was run for 40000 generations, and results summarized as above. Stochastic character mapping was used in RevBayes to obtain marginal probabilities for the number and timing of transitions events between the four-joint observed*hidden states, using 500 time slices. The script to run this analysis is provided in Appendix S2-2. To test for the potential effect of phoresy and host-type on the diversification rates of hypermetamorphic Meloidae, we ran a second analysis using the Multiple State-dependent Speciation and Extinction (MuSSE) model ( FitzJohn 2012 ). We considered three character states: (0) “bee-by-phoresy”: phoretic parasitoids on bees, i.e. lineages that feed on eggs, larvae and provisions of the bees and use adult bees as a mean of transport to the nest; (1) “bee-by-crawling”: non-phoretic parasitoids of bees, i.e. lineages that feed on eggs, larvae and provisions of the bees but use crawling on the ground to find bee nests; (2) “grasshopper-by-crawling”: non-phoretic lineages that feed on acridid eggs and use crawling to find the egg pod. Since phoresy only occurs in bee-parasitoids, it is not possible to study the joint evolution of phoresy and host-jump as two independent binary traits: i.e. there are no lineages that are “phoretic grasshopper specialists”. Settings for this analysis were similar to the BiSSE analysis, with log-normal priors for speciation and extinction rates and transitions between observed states modeled as an exponential distribution centered on a mean of 10 transitions over the tree length. A pre-burnin step of 5000 generations was used for parameter auto-tuning before a final MCMC chain length of 20000 generations. We also performed stochastic character mapping and estimation of marginal probabilities for ancestral states in the maximum a posteriori tree using the same procedure as BiSSE. The script to run this analysis is provided in Appendix S2-3. Second, we tested the robustness of the MuSSE results to Type I error or “false positives”, i.e. a trait different from changes in phoresy and host-type is leading rate heterogeneity among lineages. We ran a multiple-state HiSSE analysis (MuHiSSE), with two hidden states (A, B) associated to each of the three observed states in MuSSE (0, 1, 2); in total, the model included six states (0A, 1A, 2A, 0B, 1B, 2B). We used the same priors as in the HiSSE model above, with a log-uniform distribution for observed speciation and extinction rates, and a discretized log-normal distribution for the hidden speciation and extinction rates. A pre-burnin step of 10000 generations was used for auto-tuning before running the final MCMC chain with 40000 generations. Because mixing was difficult, we ran two analyses in parallel and merged the MCMC posterior probabilities of the two runs to increase sample size (i.e. reported results are based on this combined run). The script to run this analysis is provided in Appendix S2-4. Both standard SSE models (BiSSE, MuSSE) and hidden SSE models (HiSSE, MuHiSSE) assume that the heterogeneity in diversification rates among clades is structured in the phylogeny as an evolving trait, i.e. as discrete shifts. When this is not the case, for example, rates of diversification increasing continuously in one clade and decreasing in another clade, as a response to abiotic factors, HiSSE models can be affected by “false positives”, as well (i.e. accepting trait-dependence when none exists ( Rabosky & Goldberg 2017 ). To provide a different null model to compare against the MuSSE model, we ran a Character-Independent Diversification (CID) model ( Caetano et al. 2018 ). In the CID model, speciation and extinction rates among the states of the focal trait are constrained to be the same, but they are allowed to vary among the hidden states, in other words, rate-heterogeneity is part of the model but is unlinked to the focal observed character ( Caetano et al. 2018 ). Our CID model had two hidden states (A and B) within each of the observed character traits (0, 1 and 2) (CID2). Thus, the model included six transition rates, as well as two speciation and two relative extinction rates, corresponding to the two hidden states (those of the observed focal states are assumed to be equal): 0A=0B, 1A=1B, 2A=2B. Figures S2-S6 represent each of the models ran in this study (BiSSE, HiSSE, MuSSE, MuHiSSE, and CID2) as DAGs, indicating the parameter dependencies and prior distributions; DAGs were plotted with the graphical software GraphViz ( Ellson et al. 2004 ). We compared the fit of each model to the data using Bayes Factor comparisons of the model marginal likelihood, which were estimated via path sampling and stepping-stone sampling using parallel power posterior analyses in RevBayes ( Höhna et al. 2016 ). We ran 100 power posteriors, with a pre-burnin of 10000 generations, and a 1000-generation chain-length for each power posterior run. The script to run this analysis is provided in Appendix S2-5. Finally, we tested the robustness of our results against phylogenetic uncertainty, that is, the effect of choosing a particular random simulated subtree to represent clade diversity in the empirical-simulated phylogeny. First, using an R loop script (Appendix S2-6), we generated 100 additional empirical-simulated phylogenies, with alternative random subtrees selected and bounded for clades in the backbone phylogeny. We then checked that all simulated phylogenies belonged to the same “congruence diversification class” ( Louca & Pennell 2020 ) by calculating their “pulled speciation rate” (PSR) with the fit_hbd_psr_on_grid function, included in the R package castor ( Louca & Doebeli 2018 ); we confirmed that all simulated trees shared similar PSR mean values, allowing for a mean relative deviation of 0.5% ( Louca & Pennell 2020 ). Second, we compared posterior estimates of state-dependent speciation and extinction rates across the simulated phylogenies to ensure differences were not dependent on the topology/branch lengths of the selected random subtrees. We ran a MuSSE analysis, with the same settings as above, over the distribution of 100 simulated phylogenies; we used the multiple-processor mpi version of RevBayes ( Höhna et al. 2016 ) and a Unix bash shell script to perform a loop across phylogenies and to record the results. We then performed pairwise comparisons of the posterior distributions of the speciation, extinction, net diversification, and relative extinction across the 100 simulated-empirical trees. If the estimated difference was centered on a mean of 0, we concluded that there was no significant difference in our MuSSE estimates due to the random choice of simulated subtrees, i.e. there is not phylogenetic uncertainty effect. State-dependent analyses were performed on the Hydra supercomputer provided by the facilities of the Laboratories of Analytical Biology (LAB) of the National Museum of Natural History, Smithsonian Institution. Appendix S1 and S2 are available at https://github.com/isabelsanmartin/Trait-dependent-analyses-Meloinae RESULTS Sequencing, Assembly and Mitogenome Organization A total of 14 complete and 7 partial mitogenomes of 21 species of Meloidae were obtained. The number of reads, mean coverage, and length of each mitogenome are provided in Table S1 and Appendix S1. Genome organization was shared across all complete sequenced mitogenomes, and followed the molecular organization described by Du et al. (2016 ; 2017 ). The circular genome, represented in Figure S7, encoded for 13 protein-coding genes, and 2 rRNAs and 22 tRNA genes, and also contained a putative control region. Major strand encodes genes: cox1 , cox2 , cox3 , cytB , ATP6 , ATP8 , NAD2 , NAD3 , and NAD6 ; as well as the following tRNAs: L2 , K , D , G , A , R , N , S1 , E , T , S2 , I , M and W . Minor strand encodes genes: NAD1 , NAD4 , NAD4L , NAD1 and ribosomal 16S and 12S ; as well as the following tRNAs: F , H , P , L1 , V , Q , C and Y . Phylogenetic Inference and Divergence time estimation Phylogenetic inference was obtained using the NT-matrix or the AA+rNT datasets, and different inference methods: Maximum Likelihood in RAxML under site-homogeneous models (Fig. S8), Bayesian Inference in MrBayes under site-homogeneous models (Fig. S9), and Bayesian inference under relaxed clock models in BEAST ( Fig. 4 ). The topology exhibiting higher posterior probability values (PP) was obtained with Bayesian Inference in MrBayes and the AA+rNT dataset (88% of nodes exhibited values equal to 1) (Fig. S9). Bayesian inference with site-heterogeneous models in PhyloBayes did not converge for the more complex CAT-GTR model, due to poor mixing (results not shown). PhyloBayes results under the simpler CAT-Poisson analysis are shown in Fig S10; this tree was largely congruent with the ML, MrBayes and BEAST topologies above, recovering all major clades; however, resolution was considerably lower, especially for intertribal relationships (PP < 0.5). Subfamilies Nemognathinae and Meloinae were recovered as reciprocally monophyletic (PP =1/ BS =100) ( Fig. 4 ). Within Nemognathinae, the genus Cissites , representing the tribe Horiini, is placed as sister to all included representatives of Nemognathini (1/100). The first splitting event within Meloinae separates the tribe Mylabrini, recovered as monophyletic (1/100), as sister to a clade that includes all remaining taxa (1/89). Within this larger clade, relationships among tribes and genera were not fully resolved (0.6/28, Fig. 4 ). Tribes Eupomphini and Epicautini are also recovered as monophyletic (1/100). Meloini was found to be non-monophyletic. Species belonging to Meloe and Physomeloe form a clade (Meloini I) (1/100), closely related to Lagorina sericea (Lyttini II), and sequentially, to a clade that includes Oenas fusicornis and Lytta caragenae (LyttiniI) (1/100). In contrast, Spastonyx nemognathoides (Meloini II) was more closely related to Pyrotini and Epicautini than to Meloini I. These two subclades, Meloini I-Lyttini I-Lyttini II and Epicautini-Pyrotini-Meloini II, received high support in the Bayesian analyses, but not in the maximum likelihood analyses (1/43; 1/40, respectively). Lyttini was also recovered as non-monophyletic. Representative species and genera were arranged in three non-sister lineages: Lyttini I and Lyttini II, as mentioned above, and Berberomeloe payoyo (Lyttini III), sister to the subclade Epicautini-Pyrotini-Meloini II, albeit with low clade support ( Fig. 3 ). Bayes Factor comparison between the marginal likelihoods of the unconstrained and constrained analyses rejected the monophyly hypothesis ( H 0 ) for Epicautini plus Mylabrini in favor of the alternative hypothesis ( H 1 ), with 2lnBF= 2*((-93393.74)-(-93407.19))= 26.9, which, according to the scale given in Kass and Raftery (1995) , can be interpreted as very strong support against a monophyletic grasshopper specialists’ group (Epicautini plus Mylabrini). Estimates of lineage divergence times with relaxed molecular clocks in BEAST (MCC tree, Fig. 4 ) dated the crown-age of the family Meloidae in the Eocene (Mean 37.97 Ma, 95% HPD 29.99–55.81 Ma). An Eocene-Oligocene origin was estimated for the MRCA of the subfamily Meloinae (Mean 33.19 Ma, 26.57–48.94 Ma), whereas the MRCA of Nemognathinae was inferred as Late Oligocene-Early Miocene (Mean 24.16 Ma, 18.64–35.43 Ma). Mylabrini originated in the Late Oligocene (Mean 25.25 Ma, 19.72–37.19 Ma; the MRCAs of the remaining tribes ranged between 20 and 14 Ma) (Fig. S11). State Dependent Diversification State-dependent Speciation and Extinction Diversification analysis with BiSSE ( Fig. 5 ) supported higher speciation and extinction rates in bee-parasitoids (state 0) compared to those shown by clades feeding on grasshopper eggs (1) ( Fig. 5 ). Net diversification rates were also lower, and relative extinction rates higher, for bee parasitoids than for grasshopper specialists. There was some overlap between the marginal posterior distributions of speciation and extinction rates ( Fig. 5 ). However, pairwise comparisons of values across the MCMC posterior set produced a distribution of differences (s0–s1) in which the 95% credibility interval was larger or smaller than zero (i.e. versus overlapping zero), indicating significant differences between speciation rates or between extinction rates (bee vs grasshopper specialists) (Fig. S12). Bayesian reconstruction of ancestral states in the maximum a posteriori tree ( Fig. 5 ) showed that the ancestral condition for the MRCA of Nemognathinae and Meloinae was bee-parasitoid. Two independent events of host-jump from bees towards grasshoppers, in tribes Mylabrini and Epicautini, were recovered by stochastic character mapping analysis (Fig. S13). Download figure Open in new tab Figure 5. Maximum a posteriori reconstruction of host choice evolution in Meloidae and trait-dependent posterior distributions of diversification rates estimated through BiSSE. (A) Host choice evolution simulated under Bayesian stochastic character mapping; divergence times in millions of years are indicated by the axis at the bottom of the tree; branch colors denote different host; transitions between character states are indicated by changes in color along the branches; note that the state 0 (bee-parasitoid) is reconstructed as the ancestral state of the hypermetamorphic Meloidae and also as the ancestral state of each family. (B) Posterior densities of speciation (λ), extinction (μ), relative extinction (μ/λ) and net-diversification (λ−μ) rates. Colors correspond to the posterior probabilities for a given state; changes in host choice from state 0 to 1 are associated with diversification rate heterogeneity. Bee parasitoids lineages show higher speciation and extinction rates than grasshopper specialists, thus grasshopper specialists’ lineages are associated with the higher diversification rates and lower relative extinction rates than bee parasitoids. The HiSSE model, allowing for the existence of correlated hidden traits ( Fig. 6 ) supported a somewhat different pattern: differences in rates of speciation were still present between the two observed states, 0 and 1, but not for extinction. As in the BiSSE analysis, the net diversification rate and the rate of relative extinction were lower and higher, respectively, for bee-nests parasitoids than for grasshopper specialists, and these differences were maintained within each hidden state ( Fig. 6 ). However, differences for these rates were larger for the hidden states than for the observed states, with state B showing higher speciation rates than state A ( Fig. 6 ). In fact, for some parameters there was overlap between the marginal posterior distributions of opposite observed*hidden states, such as 1A and 0B, indicating a strong effect of the hidden trait ( Fig. 6 ). Reconstruction of ancestral states (Fig. S14a) and stochastic character mapping of transition events in the maximum a posteriori tree ( Fig. 6 , Fig. S14b) suggested that the “hidden trait” was phoresy: lineages that behave as phoretic parasitoids of bees were reconstructed as 0B (e.g. Nemognathinae, Meloe ), while non-phoretic parasitoids of bees were inferred as 0A (e.g. Pyrotini, Eupomphini). The state of the most recent common ancestor (MRCA) of subfamilies Nemognathinae and Meloinae was reconstructed as 0A, as well as the MRCA of each subfamily ( Fig. 6 ). Intriguingly, the ancestor of genus Cissites , a phoretic parasitoid on bees, as well as the ancestor of the remaining Nemognathinae, were reconstructed as 0A, with this state changing to 0B in the terminal part of the branch ( Fig. 6 ). However, marginal posterior probabilities of character states were low on the longest branches, such as the one subtending Cissites , indicating uncertainty in the reconstruction (Fig. S13b). Download figure Open in new tab Figure 6. Maximum a posteriori reconstruction of host choice evolution in Meloidae and trait-dependent posterior distributions of diversification rates estimated through HiSSE. (A) Host choice evolution simulated under Bayesian stochastic character mapping; divergence times in millions of years are indicated by the axis at the bottom of the tree; branch colors denote the four different states, being 0 and 1 the observed states (bee parasitoid and grasshopper specialists, respectively) and A and B the hidden states; transitions between character states are indicated by changes in color along the branches; note that lineages reconstructed as 0B coincide in most cases with the phoretic lineages, such as Nemognathinae and Meloini II, while non-phoretic parasitoids of bees, were reconstructed as 0A such as Lyttini II and Pyrotini; inset panel is signaling that the phoretic genus Cissites was reconstructed as 0A while the rest of the subfamily Nemognathinae was reconstructed as 0B. (B) Posterior densities of speciation (λ), extinction (μ), relative extinction (μ/λ) and net-diversification (λ−μ) rates. Colors correspond to the posterior probabilities for a given state; note that the posterior densities of speciation, relative extinction and diversification rates are in partial agreement with BiSSE results, however, the overlapping of the marginal posterior distributions for states 1A and 0B indicates that the background rate changes are unassociated with the trait in question: host-type. The importance of phoresy to explain differences in diversity related to host specificity in Meloidae was confirmed by the MuSSE analysis ( Fig. 7 ). Significant differences in speciation and extinction rates were found between bee-parasitoids and grasshopper specialists, but only for the non-phoretic lineages. Non-phoretic parasitoids of bees (state 1) exhibited marginally higher speciation rates and significantly higher extinction rates than grasshopper specialists (state 2), and also than phoretic parasitoids of bees (state 0), resulting in significantly lower diversification rates and higher relative extinction rates for state 1 than for states 0 and 2 ( Fig. 7 ). In contrast, no significant differences were found between phoretic bee-parasitoids and grasshopper specialists’ lineages, with overlapping marginal posterior distributions for states 0 and 2 across all parameters ( Fig. 7 ). Reconstruction of ancestral states (Fig. S15) and stochastic character mapping ( Fig. 7 ) supported non-phoretic bee-parasitoid as the ancestral state for the MRCA of Nemognathinae and Meloinae. Transition events from non-phoretic bee-parasitoids towards phoresy in bee-parasitoids, and towards a new grasshopper-host, took place in hypermetamorphic Meloidae in at least five independent events ( Fig. 7 ). Introducing hidden traits in this model (MuHiSSE) did not change these results ( Fig. 8 ), with hidden state B showing larger differences among the focal states than hidden state A. Likewise, our inferences were robust to the choice of subtrees for the empirical backbone-simulated phylogeny in the MuSSE analysis; pairwise comparisons across the 100 empirical-simulated phylogenies show a distribution of values with a mean centered on 0 (Fig. S16). Download figure Open in new tab Figure 7. Maximum a posteriori reconstruction of host choice evolution in Meloidae and trait-dependent posterior distributions of diversification rates estimated through MuSSE. (A) Life strategy evolution simulated under Bayesian stochastic character mapping; divergence times in millions of years are indicated by the axis at the bottom of the tree; branch colors denote different life strategies; transitions between character states are indicated by changes in color along the branches; note that a non-phoretic parasitoids of bee nest is reconstructed as the ancestral state of the hypermetamorphic Meloidae and also as the ancestral state of each family. (B) Posterior densities of speciation (λ), extinction (μ), relative extinction (μ/λ) and net-diversification (λ−μ) rates. Colors correspond to the posterior probabilities for a given state; changes in life strategy are associated with diversification rate heterogeneity. The ancestral non-phoretic parasitoids of bees’ nests showed the lowest diversification rates and the highest relative extinction rates, while there was no significant difference neither in diversification rates nor relative extinction rates between phoretic parasitoids of bee nest and parasitoids of grashoppers. Download figure Open in new tab Figure 8. Trait-dependent posterior distributions of diversification rates estimated through MuHiSSE. Posterior densities of speciation (λ), extinction (μ), net-diversification (λ−μ) rates and relative extinction (μ/λ). Colors correspond to the posterior probabilities for a given state; note that differences in posterior distributions are larger between observed states (0, 1, 2) than between hidden states (A, B). All SSE analyses ran in this study reached convergence, with traceplots showing adequate mixing and ESS values larger than 200 for all parameters (for many, this value was > 1000). The only exception was MuHiSSE, where ESS values were < 100 or 200 for many of the speciation and extinction parameters, even after combining the two runs, suggesting a possible lack of power. For this reason, the MuHiSSE model was excluded from model comparison using Bayes Factor power posteriors. The results from the Bayes Factor power posteriors, i.e. the marginal likelihood estimation for each model, are shown in Table 2 ; stepping-stone and path-sampling gave very similar values. The best-fitting model was HiSSE (ps = -16904.29), followed by MuSSE (ps=-16925.32) and BiSSE (ps = -16950.15). The CID2 model showed a worse fit to the data than any of the other SSE models (ps = -16956.41). View this table: View inline View popup Download powerpoint Table 2. Log-marginal likelihood values estimated with the Path Sampling (PS) and Stepping-stone (SS) methods for each trait-dependent diversification model applied in this study. DISCUSSION Life Strategy Evolution in Blister Beetles Our phylogenomic hypothesis supports previous works that considered parasitizing bees as the primitive life strategy of the hypermetamorphic clade of Meloidae ( Bologna & Pinto 2001 ; Bologna et al. 2008 ). The sister-taxon relationship between the reciprocally monophyletic Meloinae and Nemognathinae suggests that the origin of this strategy dates back to the Eocene ( Fig. 4 ). This general strategy has been retained by all but two of the descendent tribes, the non-sister taxa Mylabrini and Epicautini, which were able to experiment two separated non-simultaneous – homoplastic – host jump events from bees to grasshopper eggs ( Fig. 4 ). Divergence time estimation and stochastic character mapping ( Figs. 4 , 5 ; Figs. S13, S14) indicate that the first host jumping event occurred in the Early Miocene along the stem-branch of Mylabrini, while the second, took place ten million years later, in the Mid-Miocene, along the stem-branch of Epicautini. Phylogenetic evidence for the independence of the two host jumping events is also provided by previous phylogenies ( Bologna & Pinto 2001 ; Bologna et al. 2008 ). Though there is ample evidence that host changes are common in nature ( Sorenson et al. 2003 ; Wolfe et al. 2007 ; Giraud et al. 2010 ; Johnson et al. 2011 ; Forbes et al. 2017 ; Zhang et al. 2020 ), “dramatic” host-jumps, i.e. jumps to a new host species in a different family, order or phyla, are evolutionarily “rare”, compared to host shifts towards species closely related to the original host ( Engelstädter & Fortuna 2019 ; Foster 2019 ; Braga et al. 2020 ). This is because establishing a sustainable relationship with a new host species represents an important challenge for parasites, which might require new morphological or physiological adaptations ( Engelstädter & Fortuna 2019 ). Although the precise mechanism that facilitates a parasite jumping from one host to another is often unknown ( Eppinger et al. 2006 ; Lowder et al. 2009 ; Shi et al. 2014 ), jumping probability success seems limited by the phylogenetic distance between the original and the new host ( Engelstädter & Fortuna 2019 ; Foster 2019 ; Braga et al. 2020 ). Host change in Mylabrini and Epicautini represents an example of “dramatic” host-jump: a shift between insect orders, from parasitizing hymenopteran larvae to feeding on orthopteran eggs. Hymenoptera and Orthoptera are phylogenetically distant, with their respective ancestors separated by more than 350 Ma ( Misof et al. 2014 ; Song et al. 2015 ). Additionally, grasshoppers and bees exhibit very different development strategies, and their eggs differ markedly in chemical composition and properties ( Hilker & Meiners 2008 ). Grasshoppers of the family Acrididae are a dominant component of biodiversity in grassland ecosystems ( Baldi & Kisbenedek 1997 ; Badenhausser et al. 2009 ), where species of Mylabrini and Epicautini are also abundant ( Pinto 1991 ; Pinto & Bologna 1999 ; Bologna & Pinto 2002 ). Given that hosts and parasitoids share the same biome, and that food resource is abundant (i.e. grasshopper egg-pods), one could envisage continuous attempts by the parasitoid to shift to a new host, until a successful jump was achieved, not once, but twice independently over the evolutionary history of Meloinae. Other equally abundant insects in this biome are beetles of the families Tenebrionidae and Chrysomelidae. However, the highly complex hypermetamorphic life cycle of Meloidae imposes severe evolutionary constraints, limiting a potential shift towards a non-orthopteran or hymenopteran host. A regular life cycle of a species in Nemognathinae and Meloinae spans a minimum of one year ( Figs. 1 , 2 ), but some have been recorded up to six years ( Horsfall 1941 ; MacSwain 1956 ; Selander & Mathieu 1964 ; Selander & Weddle 1969 ). Stable temperature and soil moisture are the main factors governing the survival of the different larval stages ( Erickson & Werner 1974a ; Zhu et al. 2006 ). As in the case of bee nests, acridid egg-pods meet the requirements of a stable environment with constant temperature and soil humidity in which the successive larval types of the blister beetle may complete their development. This “ecological match” between the old and the new host, seems to be a common requirement for host-jump in other parasitoid lineages, such as velvet ants (Hymenoptera: Mutillidae) that were able to shift hosts from Hymenoptera to ant-nest dwelling Coleoptera with enclosed larvae (Chrysomelidae: Clythrini) ( Brothers et al. 2000 ). Intriguingly, there are reports of larvae of Cyaneolytta (Meloinae) being phoretic on Carabidae, although their feeding habits are unknown ( Di Giulio et al. 2003 ). But, host-jump is not the only relevant change that occurred along the evolution of hypermetamorphic Meloidae. A change in the way first instar larvae reaches the bee-host from active crawling on the ground, to phoresy, occurred multiple times on the history of Meloidae ( Fig. 7 , Fig. S15). This inference agrees well with Bologna & Pinto (2001) and Bologna et al. (2008) , who separately pointed out that bee-parasitism and non-phoresy were the ancestral states of hypermetamorphic Meloidae. The homoplastic nature of phoresy in Meloidae translates in that traditionally regarded as natural tribes on the basis of morphology, such as Lyttini (non-phoretic), or Meloini (phoretic), become non-monophyletic in our mitogenomic phylogeny, including now genera with different life-history strategies. Lyttini, as currently defined, is composed of three non-sister lineages at least ( Fig. 4 ): one lineage (Lyttini I) formed by Lytta ( L . caraganae ) and Oenas ( O . fusicornis ); a second lineage (Lyttini II) formed by Berbermeloe ( B. payoyo ), and a third (Lyttini III) lineage represented by Lagorina sericea . Bologna et al. (2008) already questioned the monophyly of Lyttini, and considered this tribe as a “non-phoretic taxonomic entity” where lineages with phylogenetic uncertainty were placed. Unlike Lyttini, Meloini remains monophyletic if the genus Spastonyx is excluded, but again, Meloini is not characterized by a common life-history strategy: Physomeloe is a non-phoretic bee-parasitoid, which inclusion in Meloini, supported by our phylogenetic hypothesis ( Fig. 4 ), was already suggested by Bologna et al. (2008) . The phylogenetic position of the phoretic Spastonyx has long been a mystery. This genus comprises only two small-sized species that inhabit the arid areas of northern Mexico and southern United States ( López-Estrada et al. 2018 ). Our phylogeny suggests that Spastonyx is more closely related to Pyrotini (non-phoretic) than to the rest of Meloini or Lyttini, where it has been doubtfully included ( Pinto & Selander 1970 ; Bologna & Di Giulio 2011 ). Bologna et al. (2008) suggested that phoresy was evolutionarily advantageous in Meloidae because it enhances the ability of the first instar larvae to reach the host, thus ensuring the food resource availability. Actively latching to the host, rather than wandering around to locate the nest, can also be seen as a strategy to save energy ( Baumann et al., 2018 ). In other organisms, phoresy acts as a dispersal strategy that ensures genetic exchange between populations and homogenizes them ( Athias-Binche et al. 1993 ; DiBlasi et al. 2018 ; Opatova & Št’áhlavský 2018 ). Conversely, the ancestral “bee-by-crawling” strategy would not be advantageous in a similar scenario. The first instar larva of the hypermetamorphic Meloidae is highly mobile, but unlikely able to cover great distances ( Erickson & Werner 1974a ; Selander & Weddle 1969 ). Activity periods, host searching behavior and longevity of these larvae are highly influenced by the environmental temperatures and soil moisture ( Erickson & Werner 1974b ). Without a clear preference for oviposition site ( Erickson & Werner 1974a ; Selander & Weddle 1969 ), with limited dispersal capabilities and under adverse climatic conditions, first instar larvae are forced to reach the host as quickly as possible. Under this scenario, a “wandering”, “bee-by-crawling” strategy seems the least suitable. Some studies, especially in mites ( Brown & Wilson 1992 ; Athias-Binche et al. 1993 ), suggest that phoresy may also induce speciation by “host-phoretic specialization”, in which the parasite modifies certain traits to ensure successful latching to a specific host. Some lineages such as Meloe exhibit large morphological diversity in larval traits related to phoresy ( Bologna 1983 ; 1991 ; Bologna & Pinto 2001 ), including variable development of a pygopod to crawl on vertical surfaces, or variable degree of abdominal sclerotization ( Cros 1940 ; MacSwain 1956 ; Kaszab 1969 ; Selander 1964 ; Pinto & Selander 1970 ; Bologna & Pinto 2001 ; Bologna et al. 2008 ). Meeting thus, some of the evolutionary requirements involved in a host-phoretic specialization process. Host-Jump and Phoresy Triggered Diversification in Meloidae Independent host jumping events seem to have triggered an increase in the rate of diversification in Meloinae. BiSSE and MuSSE indicate that the grasshopper specialist tribes Mylabrini and Epicautini, which are the most speciose, exhibit higher speciation rates than those of their sister taxa, non-phoretic, bee-parasitoid lineages ( Figs. 5 , 7 , Fig. S16). Evidence that host specialization can be a powerful driver of diversification comes from many different organisms, including phytophagous insects ( Ehrlich & Raven 1964 ; Forbes et al. 2017 ), pathogenic fungi ( Giraud et al. 2010 ), helminth worms ( Zietara & Lumme 2002 ), avian brood parasites ( Sorenson et al. 2003 ), and ectoparasitic arthropods ( Johnson et al. 2011 ) among others. Adaptation to the new habitat (host) can induce reproductive barriers in a relatively small number of generations ( Hendry et al. 2007 ), and thus accelerate the rate of speciation events, sometimes leading to patterns consistent with adaptive radiations ( Zietara & Lumme 2002 ; Farrel & Sequeira 2004; Fordyce 2010 ; Karvonen & Seehausen 2012 ; Forbes et al. 2017 ; Bush et al. 2019 ). Host specialization is often associated with major changes in morphology, physiology, anatomy, or reproductive mechanisms ( Barnhart et al. 2008 ; Schulte et al. 2010 ; Turrisi & Vilhelmsen 2010 ). But, so far, no set of traits shared between the larvae of grasshopper specialists Epicautini and Mylabrini has been found ( Bologna & Pinto 2001 ; Bologna et al. 2008 ). This suggests that no prior morphological adaptations or key innovations were involved in the evolution of the grasshopper specialization strategy. This is in accordance with the concept of “ecological fitting” ( Janzen 1985 ), in which the formation of a new interaction does not require the evolution of new traits; instead, it is based on traits developed for previous host-parasite interactions and “co-opted” for a new interaction given the right conditions ( Agosta 2006 ). In addition to host-jump, phoresy is another driver of diversification in hypermetamorphic blister beetles. Our SSE analyses indicate that phoretic bee-parasitoid lineages exhibit higher net diversification rates and lower relative extinction rates than non-phoretic bee-parasitoid lineages, being the latter the ancestral condition ( Figs. 6 - 8 ; Fig. S15). Although phoresy is a diversification driver, it might be acting either through host-phoretic specificity, or as a strategy to ensure panmixia and food availability ( Opatova & Št’áhlavský 2018 ), however, these alternative scenarios have not been statistically tested. Our SSE analyses demonstrate that phoresy is a life strategy as effective as a change of host to increase diversification rates in hypermetamorphic Meloidae. In fact, no significant differences in net diversification or relative extinction rates were found between phoretic bee-parasitoids and non-phoretic grasshopper specialists ( Figs. 7 , S16). Historical Contingency and Evolutionary Dead Ends in Hypermetamorphic Meloidae While trying to understand why certain groups diversify more than others, the idea of a single key factor promoting elevated diversification rates (e.g., the evolution of a morphological innovation, the invasion of a new isolated environment, or the effect of a mass extinction event depleting extant diversity) has been a dominant one in the literature ( Hodges & Arnold 1995 ; Hunter & Jernvall 1995 ; De Queiroz 2002 ; Donoghue 2005 ). Yet, despite the numerous recent studies testing for an association between trait evolution and diversification, few of them have found evidence of a single trait driving a shift in diversification rates ( Lagomarsino et al. 2017 ; Condamine et al. 2018 ; Moharrek et al. 2019 ). Instead, the dominant pattern is one in which bursts of diversification are explained by the confluence of multiple factors, sometimes acting at unison ( Donoghue & Sanderson 2015 ), sometimes in a sequence ( Donoghue 2005 ), or contingent upon one another (Givnish et al. 2015). Our study of diversification in blister beetles supports this idea of multiple interacting factors. Although phoresy and host-jump may both induce speciation by host specialization, neither host-type nor phoresy themselves can explain the pronounced differences in species richness observed among hypermetamorphic Meloidae tribes and genera. Instead, it is the parallel innovation brought by these two traits that directly impacts the diversification dynamics of Nemonagthinae and Meloinae. Though BiSSE detected a signal of trait-diversification dependency with host-type ( Fig. 5 ), HiSSE and MuSSE ( Figs. 6 , 7 ) pointed out that host-type alone cannot explain the observed differences in species richness; these results were robust against the introduction of hidden factors ( Fig. 8 ) or phylogenetic uncertainty (Fig. S16). On the other hand, these two traits have not acted simultaneously or synergistically as “synnovations” ( Donoghue 2005 ; Donoghue & Sanderson 2015 ). Non-phoretic lineages can be either bee-parasitoids or grasshopper specialists, but all phoretic lineages are bee-parasitoids (i.e. there are not phoretic grasshopper specialists). The fact that adult grasshoppers never return to the egg-pods once the oviposition takes place, makes phoresy in grasshopper specialist lineages highly unlikely. In contrast, the social behavior of bees, characterized by long-term offspring rearing, makes phoresy a highly efficient strategy to ensure food resources for lineages that are bee parasitoids ( Danforth 2007 ). Thus, in our analysis, we could not treat these two traits as independent binary characters and use an extended four-state SSE model to examine their joint evolution ( Nakov et al. 2019 ). Instead, a three-state MuSSE model ( Fig. 6 ) was used to examine how phoresy interacts with host-type in explaining rate heterogeneity between bee and grasshopper parasitoids, as detected by BiSSE ( Fig. 5 ). Our results show that phoresy and grasshopper specialization acted as two independent causal forces behind the diversification rate heterogeneity in blister beetles. Parasitoidism on bees could be seen as a peculiar example of “historical contingency” in evolution ( Gould 1989 ; Losos et al. 1998 ; Vermeij 2001 ; De Queiroz 2002 ). Beatty (2006) proposed a distinction between “causal contingency”: i.e. a change in one character that limits the evolutionary outcome of a different character, and “unpredictable contingency” ( Gould 1989 ). In our study, parasitoidism on bees fits the definition of “causal contingency” ( Beatty 2006 ; Losos et al. 1998 ): phoresy could not be achieved if hypermetamorphic Meloidae never acquired bee-parasitism as a life strategy. On the other hand, a host-jump towards feeding on orthopteran eggs would represent an example of “unpredictable contingency” ( Gould 1989 ). In other words, the onset of the “bee-by-crawling parasitoidism” in the hypermetamorphic Meloidae acted as a “historical constraint” in two different ways: as a causal factor enabling further specialization (phoresy) or by leaving an open window for unpredictability, allowing for a dramatic innovation (a host-jump) along non-sister lineages. These two “derived” life strategies in Meloinae, phoresy and grasshopper specialization, could represent adaptive radiation drivers ( Erwin 1992 ; Gavrilets & Losos 2009 ; Givnish et al. 2014 ; Simões et al. 2016 ) – seen as an “ecological opportunity” (i.e. the adaptation to a new environment such as the grasshopper egg-pods in Mylabrini and Epicautin), or as a “key innovation” (the evolution of a phoretic strategy in bee-parasitoid lineages). On the contrary, the life-history strategy of the ancestors of Meloinae and Nemonagthinae is associated with a lower net diversification rate and a significantly higher relative extinction rate than the “bee-by-phoresy” or “grasshoppers-by-crawling” strategies. Our MuSSE reconstruction ( Fig. 7 ) indicates that transitions to the latter strategies occurred at least five times independently over the evolution of the hypermetamorphic Meloidae, always accompanied by an increase in diversification rate. Under these circumstances the “bee-by-crawling” or “wandering” strategy could be tending towards an “evolutionary dead end”. Interestingly, stochastic character mapping (Fig. S15b) suggests that the “bee-by-crawling” strategy might have been ancestrally present in Cissites , a species-poor genus subtended by a long stem-branch, belonging to the species-poor tribe Horiini (15 species; Pinto & Bologna 1999 ; Bologna & Pinto 2002 ), which is the sister-group of the more speciose tribe Nemognathini ( Table 2 ). It is possible that the “bee-by-crawling” strategy was initially present within Nemognathinae, but because of the high extinction rate associated with this condition, no lineage exhibiting it presently survived in the subfamily. “Bee-by-crawling” lineages within Meloinae might also suffer a similar fate. This supports the idea that the non-phoretic bee-parasitoid lineages of Meloidae act as “depauperons”, the “flip side” of evolutionary radiations, species-poor clades with old divergence times that could be slowly deriving towards an evolutionary dead end ( Donoghue & Sanderson 2015 ). In other words, the host-jump to grasshopper eggs and the acquisition of the phoretic behavior in bee parasitoids could be evolutionary pathways by which the ancestral hypermetamorphic Meloidae lineages “escaped” extinction. The rushing forward of blister beetles to escape possible extinction determined by their ancestral complex life strategy is remarkable. The low net levels of diversification for the bee-by-crawling strategy suggest that lineages retaining the ancestral condition (bee-by-crawling parasitoids) would be slowly disappearing until extinction through progressive depauperation ( Donoghue & Sanderson 2015 ). However, the ability of blister beetles to explore new evolutionary scenarios ( López-Estrada et al. 2019 ) seems to constantly open new pathways to escape extinction (phoresy, host-jumps). It does not seem likely that selective pressures driving homoplasy acted in two characters (host-jump and phoresy) that do not set the conditions for the success of the other, unless blister beetles were constantly exploring alternative life strategies, breaking the evolutionary restrictions imposed by their phylogenetic history. In this line, López-Estrada et al. (2019) , suggested that the facility to explore different major axis of change (in morphological body ground plan) allowed a tribe of Meloidae to escape climatic extinction later in the evolution of the group. Following Vermeij (2006) reasoning, macroevolutionary change associated to life-history strategies in Meloidae is largely an effect of contingency, but since host-jump and phoresy occurred multiple times following speciation events (therefore somewhat predictable), replicate independent ecological analysis at an accessible (e.g. recent) temporal scale, are attainable. This situation opens a window to explore at a microevolutionary scale the mechanisms that originate large-scale macroevolutionary change. Hidden-SSE models are typically used to discard Type I error or “false positives” when examining a causal association between trait evolution and diversification-rate heterogeneity ( Condamine et al. 2018 ; Fernandez et al. 2018 ; Gajdzik et al. 2019 ; Nakov et al. 2019 ). Little attention is paid to negative results ( Donoghue, 2005 ), for example when the HiSSE model does not support the association between the focal trait and increased diversification rates ( Fig. 6 ). In our study, we showed that HiSSE may be useful to identify the actual “hidden” trait that is driving diversification dynamics in the group. In hypermetamorphic blister beetles, rather than host-jump, as originally hypothesized, diversification rate heterogeneity is the product of changes in life strategy, either by developing phoresy to secure the resource (food) or by a dramatic host-jump to a different order (Orthoptera). SUPPLEMENTARY MATERIAL The following supplementary figures can be found within the online supplementary material. Figure S1. The backbone ultrametric mitogenomic tree, obtained using BEAST using the AA+rNT-matrix showing the location of nodes for which we simulated subtrees with a number of taxa equal to the current taxonomic species richness (right inset). Subtrees were simulated under an integrative approach (see main text) using functions from TreePar , geiger , TESS and ape R packages. Figures S2-S6. Directional Acyclic Graphs (DAGs) showing parameter dependency and priors used in the five hierarchical Bayesian Diversification models examined here: S2: BiSSE, S3: HiSSE, S4: MuSSE, S5: MuHiSSE, S6: CD2. More details on models can be found in the main text of the article. DAGs were plotted using GraphViz. Figure S7. Gene arrangement and general structure of the mitogenome of hypermetamorphic Meloidae (subfamilies Nemognathinae and Meloinae). Figure S8. Phylogenetic hypothesis obtained using Maximum Likelihood Estimate (MLE) with RAXML on (A) NT and (B) AA+rbNT matrices. Bootstrap support values are shown above nodes. Figure S9 . Bayesian Majority-rule consensus tree obtained with MrBayes using (A) NT and (B) AA+rbNT matrices. Posterior probability (pp) clade support values are shown above nodes (∼88% of the nodes are supported with a pp value equal to 1; see Fig. 3 for additional information). Figure S10. Bayesian Majority-rule consensus tree obtained with PhyloBayes under the CT-Poisson model using (A) NT and (B) AA matrices. Posterior probability clade support values are shown above nodes. Figure S11. Lineage divergence times as estimated in BEAST using Bayesian relaxed clocks using the AA+rbNT matrix. Mean ages and 95% High Posterior Density (HPD) values are shown near to each node; violet horizontal bars represent the range of the HPD. Figure S12 . Plot of pairwise differences of speciation (A) and extinction (B) rates between states 0 and 1 estimated by BiSSE. The histogram shows the distribution of differences across the MCMC posterior distribution; the red line indicates the 0 value (no differences); black dotted lines indicate the confidence intervals. Figure S13. Results from the BiSSE model. (A) Maximum A Posteriori (MAP) tree showing reconstructed ancestral states, indicated with colors; size of circles represents marginal posterior probabilities. (B) MAP tree showing the number and timing of transition events between states reconstructed along branches using stochastic character mapping; inset to the left indicates marginal posterior probability for the inferences depicted as colors. Figure S14. Results from the HiSSE model. (A) Maximum A Posteriori (MAP) tree showing reconstructed ancestral states, indicated with colors; size of circles represents marginal posterior probabilities. (B) MAP tree showing the number and timing of transition events between states reconstructed along branches using stochastic character mapping; inset to the left indicates marginal posterior probability for the inferences depicted as colors. Figure S15. Results from the MuSSE model. (A) Maximum A Posteriori (MAP) tree showing reconstructed ancestral states, indicated with colors; size of circles represents marginal posterior probabilities. (B) MAP tree showing the number and timing of transition events between states reconstructed along branches using stochastic character mapping; inset to the left indicates marginal posterior probability for the inferences depicted as colors. Figure S16. Phylogenetic uncertainty analysis. Plot of pairwise differences of diversification rates among states (0, 1, 2) estimated by the MuSSE model. The histogram represents the distribution of these values across the 100 backbone+simulated subtrees phylogenies. The distribution is centered in 0 (represented by the red line) for states 1 and 2 (i.e. there are no differences in diversification rates across the simulated phylogenies). However, there are significant differences for the other two pairwise comparisons, confirming the results reported in Figure 6 . Table S1. List of specimens included in this study, locality data, ID number of the Entomological Collection of the Museo Nacional de Ciencias Naturales (MNCN-CSIC) and GenBank (GB) accession numbers for the mitochondrial (mt) genomes. DECLARATIONS OF INTEREST None. FUNDING This study was supported by the Spanish government Ministerio de Ciencia, Innovación y Universidades and the European Fund for Regional Development (FEDER), under grant PID2019-108109GB-I00 to IS and CGL2015-66571-P (collecting and Museum visits) and PID2019-110243GB-I00 (mitogenomic analyses) to MGP. AUTHOR CONTRIBUTIONS EKLE, IS and MGP conceived the idea. EKLE and MGP carried out the fieldwork. JEU and SA generated the molecular data. EKLE, JEU and SA assembled and annotated the mitogenomes, and performed phylogenetic inference with help from IS. EKLE and IS designed and conducted the diversification analyses. EKLE wrote the manuscript together with IS and MGP, with contributions from SA and JEU. ACKNOWLEDGMENTS We are thankful to Alberto Sánchez-Vialas, José Luis Ruiz, Ernesto Recuero, Jorge Gutiérrez-Rodríguez, Paloma Mas-Peinado, Nohemí Percino-Daniel, and Gonzalo García for their help with fieldwork; Yolanda Jiménez Ruiz for assistance in the laboratory, and Edna Gabriela López Estrada for help with R scripting. We especially thank Rogelio Delgado Román for generously providing us with additional computational resources to conduct this study. We thank the Escobar Biological Station (Cazalla de la Sierra, Spain) and all members of the Biodiversity Screening Laboratory for support during the COVID19 confinement. EKLE is supported by a doctoral scholarship from CONACyT–Mexico (330519/472100). SA has been supported by the Spanish Ministry of Science, Industry and Innovation, (BES-2014-069575) and the SponGES H2020 grant (BG-01-2015.2, agreement number 679849-2). JEU was supported by Peter Buck Postdoctoral Fellowship Program from the Smithsonian Institution (2017–2019) and is currently supported by the Atracción Talento de la Comunidad de Madrid Fellowship Program (REFF 2019-T2/ AMB-13166). Footnotes ↵ ¶ Co-senior authors ↵ § Equal contributions https://github.com/isabelsanmartin/Trait-dependent-analyses-Meloinae REFERENCES ↵ Abalde S. , Tenorio M.J. , Afonso C.M. , Uribe J.E. , Echeverry A.M. & Zardoya R . 2017 . Phylogenetic relationships of cone snails endemic to Cabo Verde based on mitochondrial genomes . BMC Evol. Biol . 17 ( 1 ): 231 . OpenUrl ↵ Abascal F. , Zardoya R. & Telford M.J . 2010 . TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations . Nucleic Acids Res . 38 ( suppl_2 ): W7 – W13 . OpenUrl CrossRef PubMed Web of Science ↵ Agosta S.J . 2006 . On ecological fitting, plant–insect associations, herbivore host shifts, and host plant selection . Oikos 114 ( 3 ): 556 – 565 . OpenUrl CrossRef Web of Science ↵ Allio R. , Scornavacca C. , Nabholz B. , Clamens A.L. , Sperling F.A. & Condamine F.L . 2020 . Whole genome shotgun phylogenomics resolves the pattern and timing of swallowtail butterfly evolution . Syst. Biol . 69 ( 1 ): 38 – 60 . OpenUrl ↵ Amor Mayor F . 1860 . Memoria sobre los insectos epispásticos de algunas provincias de España: Presentada al Colegio de Farmacéuticos de Madrid . Imprenta de Manuel Álvarez, Madrid . 36 pp. ↵ Athias-Binche F. , Schwarz H.H. & Meierhofer I . 1993 . Phoretic association of Neoseius novus (Ouds., 1902) (Acari: Uropodina) with Nicrophorus spp. (Coleoptera: Silphidae): a case of sympatric speciation? Int. J. Acarol . 19 ( 1 ): 75 – 86 . OpenUrl ↵ Badenhausser I. , Amouroux P. , Lerin J. & Bretagnolle V . 2009 . Acridid (Orthoptera: Acrididae) abundance in Western European grasslands: sampling methodology and temporal fluctuations . J. Appl. Entomol . 133 ( 9-10 ): 720 – 732 . OpenUrl Bai Y. , Li C. , Yang M. & Liang S . 2018 . Complete mitochondrial genome of the dark mealworm Tenebrio obscurus Fabricius (Insecta: Coleoptera: Tenebrionidae). Mitochondrial DNA , B , 3 ( 1 ): 171 – 172 . OpenUrl ↵ Baldi A. & Kisbenedek T . 1997 . Orthopteran assemblages as indicators of grassland naturalness in Hungary . Agr. Ecosyst. Environ . 66 : 121 – 129 . OpenUrl ↵ Barnhart M.C. , Haag W.R. & Roston W.N . 2008 . Adaptations to host infection and larval parasitism in Unionoida . J. N. Am. Benthol. Soc . 27 ( 2 ): 370 – 394 . OpenUrl CrossRef ↵ Barraclough T.G. , Vogler A.P. & Harvey P.H . 1998 . Revealing the factors that promote speciation . Philos. T. R. Soc. B . 353 ( 1366 ): 241 – 249 . OpenUrl CrossRef Web of Science Barrett C.F. , Baker W.J. , Comer J.R. , Conran J.G. , Lahmeyer S.C. , Leebens-Mack J.H. , Li J. , Lim G.S. , Mayfield-Jones D.R. , Perez L. , Medina J. , Pires J.C. , Santos C. , Stevenson D.W. , Zomlefer W.B. & Davis J.I . 2016 . Plastid genomes reveal support for deep phylogenetic relationships and extensive rate variation among palms and other commelinid monocots . New Phytol . 209 ( 2 ): 855 – 870 . OpenUrl CrossRef PubMed ↵ Baumann J. , Ferragut F. & Šimić S . 2018 . Lazy hitchhikers? Preliminary evidence for within-habitat phoresy in pygmephoroid mites (Acari , Scutacaridae). Soil Org . 90 ( 3 ): 9 – 599 . OpenUrl ↵ Beatty J . 2006 . Replaying life’s tape . J. Philos . 103 : 336 – 362 . OpenUrl CrossRef Web of Science ↵ Beaulieu J.M . 2020 . The problem with clade-specific sampling fractions . Retrieved from https://rdrr.io/cran/hisse/f/inst/doc/Clade-specific-sampling.pdf (27-09-2020) ↵ Beaulieu J.M. & O’Meara B.C . 2016 . Detecting hidden diversification shifts in models of trait-dependent speciation and extinction . Syst. Biol . 65 ( 4 ): 583 – 601 . OpenUrl CrossRef PubMed ↵ Bergsten J. , Nilsson A.N. & Ronquist F . 2013 . Bayesian tests of topology hypotheses with an example from diving beetles . Syst. Biol . 62 ( 5 ): 660 – 673 . OpenUrl CrossRef PubMed ↵ Bernt M. , Donath A. , Jühling F. , Externbrink F. , Florentz C. , Fritzsch G. , Pütz J. , Middendorf M. & Stadler P.F . 2013 . MITOS: improved de novo metazoan mitochondrial genome annotation . Mol. Phylogenet. Evol . 69 ( 2 ): 313 – 319 . OpenUrl CrossRef PubMed ↵ Bertaux B. , Prost C. , Heslan M. & Dubertret L . 1988 . Cantharide acantholysis: endogenous protease activation leading to desmosomal plaque dissolution. Brit . J. Dermatol . 118 ( 2 ): 157 – 165 . OpenUrl ↵ Blount Z.D. , Borland C.Z. & Lenski R.E . 2008 . Historical contingency and the evolution of a key innovation in an experimental population of Escherichia coli . Proc. Natl. Acad. Sci . 105 ( 23 ): 7899 – 7906 . OpenUrl Abstract / FREE Full Text ↵ Bologna M.A . 1983 . Utilizzazione dei dati biologici nella sistematica dei Meloidae (Coleoptera) . Atti XII Congresso Nazionale Italiano di Entomologia, Roma (1989) , 2 : 21 – 36 . OpenUrl ↵ Bologna M.A . 1991 . Fauna de Italia . XXVIII. Coleoptera Meloidae. Calderini, Bologna . i– xiv, 1 – 54 . ↵ Bologna M.A . 2009 . Taxonomic and biogeographical review of the Afrotropical tribe Morphozonitini (Coleoptera, Meloidae , Eleticinae) with the description of three new taxa and a key to the genera. Afr. Entomol . 17 ( 1 ): 34 – 42 . OpenUrl ↵ Bologna M.A. & Di Giulio A. 2011 . Biological and morphological adaptations in the pre-imaginal phases of the beetle family Meloidae . Atti Accad. Naz. Ital. Ent . 59 : 141 – 152 . OpenUrl ↵ Bologna M.A. , Di Giulio A. & Pinto J.D. 2002 . Review of the genus Stenodera with a description of the first instar larva of S . puncticollis (Coleoptera: Meloidae) . Eur. J. Entomol . 99 ( 3 ): 299 – 314 . OpenUrl ↵ Bologna M.A. , Oliverio M. , Pitzalis M. & Mariottini P . 2008 . Phylogeny and evolutionary history of the blister beetles (Coleoptera , Meloidae). Mol. Phylogenet. Evol . 48 ( 2 ): 679 – 693 . OpenUrl ↵ Bologna M.A. & Pinto J.D . 2001 . Phylogenetic studies of Meloidae (Coleoptera), with emphasis on the evolution of phoresy . Syst. Entomol . 26 ( 1 ): 33 – 72 . OpenUrl CrossRef ↵ Bologna M.A. & Pinto J.D . 2002 . The Old-World genera of Meloidae (Coleoptera): a key and synopsis . J. Nat. Hist . 36 : 2013 – 2102 . OpenUrl ↵ Bonett R.M. & Chippindale P.T . 2004 . Speciation, phylogeography and evolution of life history and morphology in plethodontid salamanders of the Eurycea multiplicata complex . Mol. Ecol . 13 ( 5 ): 1189 – 1203 . OpenUrl PubMed ↵ Boore J.L. , Macey J.R. & Medina M . 2005 . Sequencing and comparing whole mitochondrial genomes of animals . Methods Enzymol . 395 : 311 – 348 . OpenUrl CrossRef PubMed Web of Science ↵ Braga M.P. , Landis M.J. , Nylin S. , Janz N. & Ronquist F . 2020 . Bayesian inference of ancestral host-parasite interactions under a phylogenetic model of host repertoire evolution . Syst. Biol . doi: https://doi.org/10.1093/sysbio/syaa019 ↵ Bravo C. , Bautista L.M. , García-París M. , Blanco G. , Alonso J.C . 2014 . Males of a strongly polygynous species consume more poisonous food than females . PLOS One 9 : e111057 . OpenUrl ↵ Bravo C. , Mas-Peinado P. , Bautista L.M. , Blanco G. , Alonso J.C. & García-París M . 2017 . Cantharidin is conserved across phylogeographic lineages and present in both morphs of Iberian Berberomeloe blister beetles (Coleoptera , Meloidae). Zool. J. Linn. Soc . 180 ( 4 ): 790 – 804 . OpenUrl ↵ Brooks D.R. , León-Règagnon V. , McLennan D.A. & Zelmer D . 2006 . Ecological fitting as a determinant of the community structure of platyhelminth parasites of anurans . Ecology 87 ( sp7 ): S76 – S85 . OpenUrl CrossRef PubMed Web of Science ↵ Brothers D.J. , Tschuch G. & Burger F . 2000 . Associations of mutillid wasps (Hymenoptera, Mutillidae) with eusocial insects. Insectes sociaux 47 ( 3 ): 201 – 211 . OpenUrl ↵ Brown J.M. & Wilson D.S . 1992 . Local specialization of phoretic mites on sympatric carrion beetle hosts . Ecology 73 ( 2 ): 463 – 478 . OpenUrl ↵ Bush S.E. , Villa S.M. , Altuna J.C. , Johnson K.P. , Shapiro M.D. & Clayton D.H . 2019 . Host defense triggers rapid adaptive radiation in experimentally evolving parasites . Evolution 3 ( 2 ): 120 – 128 . OpenUrl ↵ Caetano D.S. , O’Meara B.C. & Beaulieu J.M . 2018 . Hidden state models improve state-dependent diversification approaches, including biogeographical models . Evolution 72 ( 11 ): 2308 – 2324 . OpenUrl CrossRef PubMed ↵ Cai L. , Xi Z. , Lemmon E.M. , Lemmon A.R. , Mast A. , Buddenhagen C.E. , Liang L. & Davis C.C . 2020 . The perfect storm: Gene tree estimation error, incomplete lineage sorting, and ancient gene flow explain the most recalcitrant ancient angiosperm clade , Malpighiales. bioRxiv (preprint) . doi: https://doi.org/10.1101/2020.05.26.112318 ↵ Calatayud J. , Hórreo J.L. , Madrigal-González J. , Migeon A. , Rodríguez M.Á. , Magalhães S. & Hortal J . 2016 . Geography and major host evolutionary transitions shape the resource use of plant parasites . Proc. Natl. Acad. Sci . 113 ( 35 ): 9840 – 9845 . OpenUrl Abstract / FREE Full Text Cameron S.L. , Sullivan J. , Song H. , Miller K.B. & Whiting M.F . 2009 . A mitochondrial genome phylogeny of the Neuropterida (lace-wings, alderflies and snakeflies) and their relationship to the other holometabolous insect orders . Zool. Scripta 38 ( 6 ): 575 – 590 . OpenUrl CrossRef ↵ Campos-Soldini M.P. , Safenraiter M.E. , Wagner L.S. , Fernández E.N. , & Sequín C.J . 2018 . Checklist of Epicauta Dejean from America (Meloidae, Meloinae , Epicautini). Zookeys 807 : 14 – 124 OpenUrl ↵ Carrel J.E. & Eisner T . 1974 . Cantharidin: potent feeding deterrent to insects . Science 183 : 755 – 757 . OpenUrl Abstract / FREE Full Text ↵ Castresana J . 2000 . Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis . Mol. Biol. Evol . 17 ( 4 ): 540 – 552 . OpenUrl CrossRef PubMed Web of Science ↵ Chan P.P. & Lowe T.M . 2019 . tRNAscan-SE: Searching for tRNA Genes in Genomic Sequences . Methods Mol. Biol . 1962 : 1 – 14 . OpenUrl CrossRef PubMed ↵ Condamine F.L. , Rolland J. , Höhna S. , Sperling F.A. , & Sanmartín I . 2018 . Testing the role of the Red Queen and Court Jester as drivers of the macroevolution of Apollo butterflies . Syst. Biol . 67 ( 6 ): 940 – 964 . OpenUrl ↵ Cros A . 1940 . Essai de classification des Meloides algériens. VI Congreso Internacional de Entomologia, Madrid (1935) . pp. 312–338. Laboratorio de Entomología del Museo Nacional de Ciencias Naturales, Madrid . 1 – 961 . ↵ Culshaw V. , Stadler T. & Sanmartín I . 2019 . Exploring the power of Bayesian birth-death skyline models to detect mass extinction events from phylogenies with only extant taxa . Evolution 73 ( 6 ): 1133 – 1150 . OpenUrl ↵ Danforth B . 2007 . Bees . Curr. Biol . 17 ( 5 ): R156 – R161 . OpenUrl CrossRef PubMed Web of Science ↵ De Queiroz A. 2002 . Contingent predictability in evolution: key traits and diversification . Syst. Biol . 51 ( 6 ): 917 – 929 . OpenUrl CrossRef PubMed Web of Science ↵ Denier P.C.L . 1935 . Coleopterorum americanorum familiae meloidarum . Enumeratio synonymica. Rev. Soc. Ent. Arg . 7 : 139 – 176 . OpenUrl ↵ Di Giulio A. , Aberlenc H.P. , Vigna Taglianti A. & Bologna M.A. 2003 . Definition and description of larval types of Cyaneolytta (Coleoptera Meloidae) and new records of their phoretic association with Carabidae (Coleoptera) . Trop. Zool . 16 ( 2 ): 165 – 187 . OpenUrl ↵ DiBlasi E. , Johnson K.P. , Stringham S.A. , Hansen A.N. , Beach A.B. , Clayton D.H. & Bush S.E . 2018 . Phoretic dispersal influences parasite population genetic structure . Mol. Ecol . 27 ( 12 ): 2770 – 2779 . OpenUrl ↵ Dioscorides P. 1636 . Acerca de la materia medicinal y de los venenos mortíferos. Traduzido de la lengua griega, en la vulgar castellana & ilustrado con claras y sustanciales anotaciones, y con las figuras de innumerables platas exquisitas y raras, por el Doctor Andres de Laguna, Médico de Iulio Tercero Pont. Max. Miguel Sorolla, Valencia . 11 + 616 + 26 pp. ↵ Donoghue M.J . 2005 . Key innovations, convergence, and success: macroevolutionary lessons from plant phylogeny . Paleobiology 31 ( 2 ): 77 – 93 . OpenUrl Abstract / FREE Full Text ↵ Donoghue M.J. & Edwards E.J . 2014 . Biome shifts and niche evolution in plants . Annu. Rev. Ecol. Evol. Syst . 45 : 547 – 572 . OpenUrl CrossRef ↵ Donoghue M.J. & Sanderson M.J . 2015 . Confluence, synnovation, and depauperons in plant diversification . New Phytol . 207 ( 2 ): 260 – 274 . OpenUrl CrossRef PubMed ↵ Du C. , He S. , Song X. , Liao Q. , Zhang X. & Yue B. 2016 . The complete mitochondrial genome of Epicauta chinensis (Coleoptera: Meloidae) and phylogenetic analysis among coleopteran insects . Gene 578 ( 2 ): 274 – 280 . OpenUrl ↵ Du C. , Zhang L. , Lu T. , Ma J. , Zeng C. , Yue B. & Zhang X. 2017 . Mitochondrial genomes of blister beetles (Coleoptera , Meloidae) and two large intergenic spacers in Hycleus genera. BMC Genomics 18 ( 1 ): 698 . OpenUrl ↵ Drummond A.J. & Rambaut A . 2007 . BEAST: Bayesian evolutionary analysis by sampling trees . BMC Evol. Biol . 7 : 214 . OpenUrl CrossRef PubMed ↵ Drummond C.S. , Eastwood R.J. , Miotto S.T. & Hughes C.E . 2012 . Multiple continental radiations and correlates of diversification in Lupinus (Leguminosae): testing for key innovation with incomplete taxon sampling . Syst. Biol . 61 ( 3 ): 443 – 460 . OpenUrl CrossRef PubMed ↵ Ehrlich P.R. , Raven P.H . 1964 . Butterflies and plants: a study in coevolution . Evolution 18 : 586 – 608 . OpenUrl CrossRef Web of Science ↵ Mutzel P. Jünger M Ellson J. , Gansner E.R. , Koutsofios E. , North S.C. & Woodhull G . 2004 . Graphviz and dynagraph—static and dynamic graph drawing tools . Pp. 127 – 148 . In: Mutzel P. & Jünger M . (eds.). Graph drawing software . Springer , Berlin, Heidelberg . 364 pp. ↵ Engelstädter J. & Fortuna N.Z . 2019 . The dynamics of preferential host switching: Host phylogeny as a key predictor of parasite distribution . Evolution 73 ( 7 ): 1330 – 1340 . OpenUrl ↵ Enns W.R. 1956 . A revision of the genera Nemognatha, Zonitis and Pseudozonitis (Coleoptera, Meloidae) in America North of Mexico, with a proposed new genus . Univ. Kan. Sci. Bull . 37 : 685 – 909 . OpenUrl ↵ Eppinger M. , Baar C. , Linz B. , Raddatz G. , Lanz C. , Keller H. , Morelli G. , Gressmann H. , Achtman M. & Schuster S.C . 2006 . Who ate whom? Adaptive Helicobacter genomic changes that accompanied a host jump from early humans to large felines . PLoS Genet . 2 ( 7 ): e120 . OpenUrl CrossRef PubMed ↵ Erickson E.H. & Werner F.G . 1974a . Bionomics of Nearctic bee-associated Meloidae (Coleoptera); life histories and nutrition of certain Meloinae . Ann. Entomol. Soc. Am . 67 ( 3 ): 394 – 400 . OpenUrl CrossRef ↵ Erickson E.H. & Werner F.G . 1974b . Bionomics of Nearctic bee-associated Meloidae (Coleoptera). A comparative analysis of larval host-seeking behavior among the Meloinae and Nemognathinae . Ann. Entomol. Soc. Am . 67 ( 6 ): 903 – 908 . OpenUrl CrossRef ↵ Erwin D.H . 1992 . A preliminary classification of evolutionary radiations . Hist. Biol . 6 ( 2 ): 133 – 147 . OpenUrl Farrell B.D. & Sequeira A.S . 2004 . Evolutionary rates in the adaptive radiation of beetles on plants . Evolution 58 ( 9 ): 1984 – 2001 . OpenUrl PubMed Web of Science ↵ Fernandez R. , Kallal R.J. , Dimitrov D. , Ballesteros J.A. , Arnedo M.A. , Giribet G. , & Hormiga G . 2018 . Phylogenomics, diversification dynamics, and comparative transcriptomics across the spider tree of life . Curr. Biol . 28 ( 9 ): 1489 – 1497 . OpenUrl CrossRef ↵ Finstermeier K. , Zinner D. , Brameier M. , Meyer M. , Kreuz E. , Hofreiter M. & Roos C . 2013 . A mitogenomic phylogeny of living primates . PloS One . 8 ( 7 ): e69504 . OpenUrl CrossRef PubMed ↵ Fischer J.B. 1827 . Tentamen conspectus Cantharidiarum. Disertatio inauguralis, quam pro summis in Medicina et Chirurgia honoribus legitime obtinendis eruditorum examini subjicit. Lindauer, Monachii . 6 + 26 pp. ↵ FitzJohn R.G . 2010 . Quantitative traits and diversification . Syst. Biol . 59 ( 6 ): 619 – 633 . OpenUrl CrossRef PubMed Web of Science ↵ FitzJohn R.G . 2012 . Diversitree: comparative phylogenetic analyses of diversification in R . Methods Ecol. Evol . 3 : 1084 – 1092 . OpenUrl CrossRef ↵ FitzJohn R.G. , Maddison W.P. & Otto S.P . 2009 . Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies . Syst. Biol . 58 ( 6 ): 595 – 611 . OpenUrl CrossRef PubMed Web of Science ↵ Forbes A.A. , Devine S.N. , Hippee A.C. , Tvedte E.S. , Ward A.K. , Widmayer H.A. & Wilson C.J . 2017 . Revisiting the particular role of host shifts in initiating insect speciation . Evolution 71 ( 5 ): 1126 – 1137 . OpenUrl CrossRef ↵ Fordyce J.A . 2010 . Host shifts and evolutionary radiations of butterflies . P. Roy. Soc. B-Biol. Sci . 277 ( 1701 ): 3735 – 3743 . OpenUrl ↵ Foster C.S . 2019 . Digest: The phylogenetic distance effect: Understanding parasite host switching . Evolution 73 ( 7 ): 1494 – 1495 . OpenUrl ↵ Freyman W.A. & Höhna S . 2019 . Stochastic character mapping of state-dependent diversification reveals the tempo of evolutionary decline in self-compatible Onagraceae lineages . Syst. Biol . 68 ( 3 ): 505 – 519 . OpenUrl ↵ Gajdzik L. , Aguilar-Medrano R. & Frédérich B . 2019 . Diversification and functional evolution of reef fish feeding guilds . Ecol. Lett . 22 ( 4 ): 572 – 582 . OpenUrl ↵ Gandon S. & Michalakis Y . 2002 . Local adaptation, evolutionary potential and host– parasite coevolution: interactions between migration, mutation, population size and generation time . J. Evol. Biol . 15 ( 3 ): 451 – 462 . OpenUrl CrossRef Web of Science ↵ Gavrilets S. & Losos J.B . 2009 . Adaptive radiation: contrasting theory with data . Science 323 : 732 – 737 . OpenUrl Abstract / FREE Full Text ↵ Gibb G.C. , Condamine F.L. , Kuch M. , Enk J. , Moraes-Barros N. , Superina M. , Poinar H.N. & Delsuc F . 2016 . Shotgun mitogenomics provides a reference phylogenetic framework and timescale for living xenarthrans . Mol. Biol. Evol . 33 ( 3 ): 621 – 642 . OpenUrl CrossRef PubMed ↵ Giraud T. , Gladieux P. & Gavrilets S . 2010 . Linking the emergence of fungal plant diseases with ecological speciation . Trends Ecol. Evol . 25 : 387 – 395 . OpenUrl CrossRef PubMed Web of Science Givnish T.J . 2015 . Adaptive radiation versus “radiation” and “explosive diversification”: why conceptual distinctions are fundamental to understanding evolution . New Phytol . 207 ( 2 ): 297 – 303 . OpenUrl CrossRef PubMed ↵ Givnish T.J. , Barfuss M.H. , Van Ee B. , Riina R. , Schulte K. , Horres R. , Gonsiska P.A. , Jabaily R.S. , Crayn D.M. , Smith J.A.C. , Winter K. , Brown G.K. , Evans T.M. , Holst B.K. , Luther H. , Till W. , Zizka G. , Berry P.E. & Sytsma K.J. 2014 . Adaptive radiation, correlated and contingent evolution, and net species diversification in Bromeliaceae . Mol. Phylogenet. Evol . 71 : 55 – 78 . OpenUrl CrossRef PubMed ↵ Göker M. , Riethmüller A. , Voglmayr H. , Weiss M. & Oberwinkler F . 2004 . Phylogeny of Hyaloperonospora based on nuclear ribosomal internal transcribed spacer sequences . Mycol. Prog . 3 : 83 – 94 . OpenUrl Goldberg E.E. & Igić B . 2008 . On phylogenetic tests of irreversible evolution . Evolution . 62 ( 11 ): 2727 – 2741 . OpenUrl CrossRef PubMed Web of Science ↵ Goldberg E.E. , Otto S.P. , Vamosi J.C. , Mayrose I. , Sabath N. , Ming R. & Ashman T.L . 2017 . Macroevolutionary synthesis of flowering plant sexual systems . Evolution 71 ( 4 ): 898 – 912 . OpenUrl ↵ Gould S.J . 1989 . Wonderful Life. Norton, New York , London . 352 pp. ↵ Hafernik J. & Saul-Gershenz L . 2000 . Beetle larvae cooperate to mimic bees . Nature 405 : 35 – 36 . OpenUrl PubMed ↵ Hardy N.B. & Otto S.P . 2014 . Specialization and generalization in the diversification of phytophagous insects: tests of the musical chairs and oscillation hypotheses . P. Roy. Soc. B-Biol. Sci . 281 ( 1795 ): 2013 – 2960 . OpenUrl ↵ Harmon L.J. , Weir J.T. , Brock C. , Glor R.E. & Challenger W.E . 2008 . GEIGER: Investigating evolutionary radiations . Bioinformatics 24 : 129 – 131 . OpenUrl CrossRef PubMed Web of Science ↵ Hasegawa M. , Kishino H. & Yano T.A . 1985 . Dating of the human-ape splitting by a molecular clock of mitochondrial DNA . J. Mol. Evol . 22 ( 2 ): 160 – 174 . OpenUrl CrossRef PubMed Web of Science ↵ Hendry A.P. , Nosil P. & Rieseberg L.H . 2007 . The speed of ecological speciation . Funct. Ecol . 21 : 455 – 464 . OpenUrl ↵ Herrera-Alsina L. , van Els P. & Etienne R.S. 2019 . Detecting the dependence of diversification on multiple traits from phylogenetic trees and trait data . Syst. Biol . 68 ( 2 ): 317 – 328 . OpenUrl ↵ Hilker M. & Meiners T . 2008 . Chemoecology of insect eggs and egg deposition. Blackwell Publishing , Oxford, Malden, Carlton Victoria, Berlin . 416 pp. ↵ Hodges S.A. & Arnold M.L . 1995 . Spurring plant diversification: are floral nectar spurs a key innovation? P. Roy. Soc. B-Biol. Sci . 262 ( 1365 ): 343 – 348 . OpenUrl ↵ Höhna S . 2013 . Fast simulation of reconstructed phylogenies under global time-dependent birth-death processes . Bioinformatics 29 ( 11 ): 1367 – 1374 . OpenUrl CrossRef PubMed Web of Science ↵ Höhna S . 2014 . Likelihood inference of non-constant diversification rates with incomplete taxon sampling . PLoS One . 9 ( 1 ): e84184 . OpenUrl CrossRef PubMed ↵ Höhna S. , Freyman W.A. , Nolen Z. , Huelsenbeck J.P. , May M.R. & Moore B.R . 2019 . A Bayesian approach for estimating branch-specific speciation and extinction rates . bioRxiv (preprint) doi: https://doi.org/10.1101/555805 ↵ Höhna S. , Heath T.A. , Boussau B. , Landis M.J. , Ronquist F. & Huelsenbeck J.P . 2014 . Probabilistic graphical model representation in phylogenetics . Syst. Biol . 63 ( 5 ): 753 – 771 . OpenUrl CrossRef PubMed ↵ Höhna S. , Landis M.J. , Heath T.A. , Boussau B. , Lartillot N. , Moore B.R. , Huelsenbeck J.P. & Ronquist F . 2016 . RevBayes: Bayesian phylogenetic inference using graphical models and an interactive model-specification language . Syst. Biol . 65 : 726 – 736 . OpenUrl CrossRef PubMed ↵ Horsfall W.R . 1941 . Biology of the black blister beetle (Coleoptera: Meloidae) . Ann. Entomol. Soc. Am . 34 ( 1 ): 114 – 126 . OpenUrl CrossRef ↵ Hunter J.P. & Jernvall J . 1995 . The hypocone as a key innovation in mammalian evolution . Proc. Natl. Acad. Sci . 92 ( 23 ): 10718 – 10722 . OpenUrl Abstract / FREE Full Text ↵ Irisarri I. , Uribe J.E. , Eernisse D.J. & Zardoya R . 2020 . A mitogenomic phylogeny of chitons (Mollusca: Polyplacophora) . BMC Evol. Biol . 20 ( 1 ): 1 – 15 . OpenUrl Iturralde-Vincent M.A. & Mac-Phee R.D.E . 1996 . Age and Paleogeographic origin of Dominican amber . Science 273 : 1850 – 1852 . OpenUrl Abstract / FREE Full Text ↵ Janzen D.H . 1985 . On ecological fitting . Oikos 45 ( 3 ): 308 – 310 . OpenUrl CrossRef PubMed Web of Science ↵ Johnson K.P. , Weckstein J.D. , Meyer M.J. & Clayton D.H . 2011 . There and back again: switching between host orders by avian body lice (Ischnocera: Goniodidae) . Biol. J. Linn. Soc . 102 : 614 – 625 . OpenUrl ↵ Kaiser E. & Michl H. 1958 . Die Biochemie der Tierischen Gifte . F. Deuticke , Wisconsin . 257 pp. ↵ Karvonen A. & Seehausen O . 2012 . The role of parasitism in adaptive radiations–when might parasites promote and when might they constrain ecological speciation? Int . J. Ecol . 2012 : 1 – 20 . OpenUrl ↵ Kass R.E. & Raftery A.E . 1995 . Bayes factors . J. Am. Stat. Assoc . 90 ( 430 ): 773 – 795 . OpenUrl CrossRef PubMed Web of Science ↵ Kaszab Z. , 1969 . The system of the Meloidae (Coleoptera) . Mem. Soc. Entomol. Ital . 48 : 241 – 248 . OpenUrl ↵ Katoh K. , Kuma K. , Toh H. & Miyata T . 2005 . MAFFT version 5: improvement in accuracy of multiple sequence alignment . Nucleic Acids Res . 33 ( 2 ): 511 – 518 . OpenUrl CrossRef PubMed Web of Science ↵ Katoh K. , Rozewicki J. & Yamada K.D . 2017 . MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization . Brief. Bioinform . bbx108 ↵ Kriebel R. , Drew B. , González-Gallegos J.G. , Celep F. , Heeg L. , Mahdjoub M.M. & Sytsma K.J . 2020 . Pollinator shifts, contingent evolution, and evolutionary constraint drive floral disparity in Salvia (Lamiaceae): evidence from morphometrics and phylogenetic comparative methods . Evolution 74 ( 7 ): 1335 – 1355 . OpenUrl ↵ Lagomarsino L.P. , Forrestel E.J. , Muchhala N. & Davis C.C . 2017 . Repeated evolution of vertebrate pollination syndromes in a recently diverged Andean plant clade . Evolution 71 ( 8 ): 1970 – 1985 . OpenUrl ↵ Lanfear R. , Frandsen P.B. , Wright A.M. , Senfeld T. & Calcott B . 2016 . PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses . Mol. Biol. Evol . 34 ( 3 ): 772 – 773 . OpenUrl ↵ Lartillot N. , Lepage T. & Blanquart S . 2009 . PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating . Bioinformatics 25 : 2286 – 2288 OpenUrl CrossRef PubMed Web of Science ↵ Lartillot N. & Philippe H . 2004 . A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process . Mol. Biol. Evol . 21 ( 6 ): 1095 – 1109 . OpenUrl CrossRef PubMed Web of Science Li-Na L. & Cheng-Ye W . 2014 . Complete mitochondrial genome of yellow meal worm ( Tenebrio molitor ) . Zool. Res . 35 ( 6 ): 537 – 545 . OpenUrl ↵ Liu Y. , Li H. , Song F. , Zhao Y. , Wilson J.J. & Cai W . 2019 . Higher-level phylogeny and evolutionary history of Pentatomomorpha (Hemiptera: Heteroptera) inferred from mitochondrial genome sequences . Syst. Entomol . 44 ( 4 ): 810 – 819 . OpenUrl ↵ López-Estrada E.K. , Pérez-Flores O. , Zaldívar-Riverón A. & García-París M . 2018 . Notes on the geographic distribution of Spastonyx nemognathoides Selander, 1954 (Coleoptera: Meloidae) in North America (Mexico and USA) . Pan-Pac Entomol . 94 ( 2 ): 55 – 58 . OpenUrl ↵ López-Estrada E.K. , Sanmartín I. , García-París M. & Zaldívar-Riverón A . 2019 . High extinction rates and non-adaptive radiation explains patterns of low diversity and extreme morphological disparity in North American blister beetles (Coleoptera , Meloidae). Mol. Phylogenet. Evol . 130 : 156 – 168 . OpenUrl ↵ Losos J.B. , Jackman T.R. , Larson A. , de Queiroz K. & Rodríguez-Schettino L. 1998 . Contingency and determinism in replicated adaptive radiations of island lizards . Science 279 : 2115 – 2118 . OpenUrl Abstract / FREE Full Text ↵ Losos J.B. & Miles D.B . 2002 . Testing the hypothesis that a clade has adaptively radiated: iguanid lizard clades as a case study . Am. Nat . 160 ( 2 ): 147 – 157 . OpenUrl CrossRef PubMed Web of Science ↵ Louca S. & Doebeli M . 2018 . Efficient comparative phylogenetics on large trees . Bioinformatics 34 ( 6 ): 1053 – 1055 . OpenUrl ↵ Louca S. & Pennell M.W . 2020 . Extant timetrees are consistent with a myriad of diversification histories . Nature 580 ( 7804 ): 502 – 505 . OpenUrl CrossRef PubMed ↵ Lowder B.V. , Guinane C.M. , Zakour N.L.B. , Weinert L.A. , Conway-Morris A. , Cartwright R.A. , Simpson A.J. , Rambaut A. , Nübel U. & Fitzgerald J.R . 2009 . Recent human-to-poultry host jump, adaptation, and pandemic spread of Staphylococcus aureus . Proc. Natl. Acad. Sci . 106 ( 46 ): 19545 – 19550 . OpenUrl Abstract / FREE Full Text ↵ MacSwain J.W . 1956 . A classification of the first instar larvae of the Meloidae (Coleoptera) . Univ. Calif. Publ. Entomol . 12 : 1 – 182 . OpenUrl ↵ Maddison W.P. & FitzJohn R.G . 2015 . The unsolved challenge to phylogenetic correlation tests for categorical characters . Syst. Biol . 64 ( 1 ): 127 – 136 . OpenUrl CrossRef PubMed ↵ Maddison W.P. , Midford P.E. & Otto S.P . 2007 . Estimating a binary character’s effect on speciation and extinction . Syst. Biol . 56 ( 5 ): 701 – 710 . OpenUrl CrossRef PubMed Web of Science ↵ Magallón S. & Sanderson M.J . 2001 . Absolute diversification rates in angiosperm clades . Evolution 55 ( 9 ): 1762 – 1780 . OpenUrl CrossRef PubMed Web of Science ↵ May M.R. , Höhna S. & Moore B.R . 2016 . A Bayesian approach for detecting the impact of mass-extinction events on molecular phylogenies when rates of lineage diversification may vary . Methods Ecol. Evol . 7 ( 8 ): 947 – 959 . OpenUrl CrossRef PubMed ↵ May M.R. & Moore B.R . 2020 . A Bayesian Approach for Inferring the Impact of a Discrete Character on Rates of Continuous-Character Evolution in the Presence of Background-Rate Variation . Syst. Biol . 69 ( 3 ): 530 – 544 . OpenUrl ↵ Mayr E Schu-z E Miller A.H. 1949 . Some ecologic and morphologic considerations in the evolution of higher taxonomic categories. Pp 84–89 . In: Mayr E . & Schu-z E . (eds.). Ornithologie als Biologische Wissenschaft . Carl Winter, Cornell University . 291 pp. ↵ Miller M.A. , Pfeiffer W. & Schwartz T. 2010 . Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Pp. 45–52 . In: 2010 Gateway Computing Environments Workshop (GCE 2010) . Institute of Electrical and Electronics Engineers (IEEE) , New Orleans . 106 pp. ↵ Misof B. , Liu S. , Meusemann K. , Peters R.S. , Donath A. , Mayer C. , Frandsen P.B. , Ware J. , Flouri T. , Beutel R.G. , Niehuis O. , Petersen M. , Izquierdo-Carrasco F. , Wappler T. , Rust J. , Aberer A.J. , Aspöck U. , Aspöck H. , Bartel D. , Blanke A. , Berger S. , Böhm A. , Buckley T.R. , Calcott B. , Chen J. , Friedrich F. , Fukui M. , Fujita M. , Greve C. , Grobe P. , Gu S. , Huang Y. , Jermiin L.S. , Kawahara A.Y. , Krogmann L. , Kubiak M. , Lanfear R. , Letsch H. , Li Y. , Li Z. , Li J. , Lu H. , Machida R. , Mashimo Y. , Kapli P. , McKenna D.D. , Meng G. , Nakagaki Y. , Navarrete-Heredia J.L. , Ott M. , Ou Y. , Pass G. , Podsiadlowski L. , Pohl H. , Reumont von B.M. , Schütte K. , Sekiya K. , Shimizu S. , Slipinski A. , Stamatakis A. , Song W. , Su X. , Szucsich N.U. , Tan M. , Tan X. , Tang M. , Tang J. , Timelthaler G. , Tomizuka S. , Trautwein M. , Tong X. , Uchifune T. , Walzl M.G. , Wiegmann B.M. , Wilbrandt J. , Wipfler B. , Wong T.K.F. , Wu Q. , Wu G. , Xie Y. , Yang S. , Yang Q. , Yeates D.K. , Yoshizawa K. , Zhang Q. , Zhang R. , Zhang W. , Zhang Y. , Zhao J. , Zhou C. , Zhou L. , Ziesmann T. , Zou S. , Li Y. , Xu X. , Zhang Y. , Yang H. , Wang J. , Wang J. , Kjer K.M. , Zhou X . 2014 . Phylogenomics resolves the timing and pattern of insect evolution . Science 346 : 763 – 767 . OpenUrl Abstract / FREE Full Text ↵ Moharrek F. , Sanmartín I. , Kazempour-Osaloo S. & Nieto Feliner G . 2019 . Morphological innovations and vast extensions of mountain habitats triggered rapid diversification within the species-rich Irano-Turanian genus Acantholimon (Plumbaginaceae) . Front. Genet . 9 : 698 . OpenUrl ↵ Mooers O. & Heard S.B . 1997 . Inferring Evolutionary process from phylogenetic tree shape . Q. Rev. Biol . 71 ( 1 ): 31 – 54 . OpenUrl ↵ Moore B.R. , Höhna S. , May M.R. , Rannala B. & Huelsenbeck J.P . 2016 . Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures . Proc. Natl. Acad. Sci . 113 ( 34 ): 9569 – 9574 . OpenUrl Abstract / FREE Full Text ↵ Morlon H . 2014 . Phylogenetic approaches for studying diversification . Ecol. Lett . 17 ( 4 ): 508 – 525 . OpenUrl CrossRef PubMed ↵ Morlon H. , Hartig F. & Robin S . 2020 . Prior hypotheses or regularization allow inference of diversification histories from extant timetrees . bioRxiv . doi: https://doi.org/10.1101/2020.07.03.185074 ↵ Morlon H. , Parsons T.L. & Plotkin J.B . 2011 . Reconciling molecular phylogenies with the fossil record . Proc. Natl. Acad. Sci . 108 ( 39 ): 16327 – 16332 . OpenUrl Abstract / FREE Full Text ↵ Nakov T. , Beaulieu J.M. & Alverson A.J . 2019 . Diatoms diversify and turn over faster in freshwater than marine environments . Evolution 73 ( 12 ): 2497 – 2511 . OpenUrl Nie R. , Vogler A.P. , Yang X.K. & Lin M . 2020 . Higher-level phylogeny of longhorn beetles (Coleoptera: Chrysomeloidea) inferred from mitochondrial genomes . Syst. Entomol . doi: https://doi.org/10.1111/syen.12447 ↵ Nylin S. , Agosta S. , Bensch S. , Boeger W.A. , Braga M.P. , Brooks D.R. , Forister M.L. , Hämback P.A. , Hoberg E.P. , Nyman T. , Schäpers A. , Stigall A.L. , Wheat C.W. , Österling M. & Janz N . 2018 . Embracing colonizations: a new paradigm for species association dynamics . Trends Ecol. Evol . 33 ( 1 ): 4 – 14 . OpenUrl ↵ Opatova V. & Št’áhlavský F . 2018 . Phoretic or not? Phylogeography of the pseudoscorpion Chernes hahnii (Pseudoscorpiones: Chernetidae) . J. Arachnol . 46 ( 1 ): 104 – 113 . OpenUrl ↵ Paradis E. , Bolker B. & Strimmer K . 2004 . APE: analysis of phylogenetics and evolution in R language . Bioinformatics 20 ( 2 ): 289 – 290 . OpenUrl CrossRef PubMed Web of Science ↵ Percino-Daniel N. , Buckley D. & García-París M . 2013 . Pharmacological properties of blister beetles (Coleoptera: Meloidae) promoted their integration into the cultural heritage of native rural Spain as inferred by vernacular names diversity, traditions, and mitochondrial DNA . J. Ethnopharmacol . 147 ( 3 ): 570 – 583 . OpenUrl ↵ Philippe H. , Brinkmann H. , Lavrov D.V. , Littlewood D.T.J. , Manuel M. , Wörheide G. & Baurain D . 2011 . Resolving difficult phylogenetic questions: why more sequences are not enough . PLoS Biology 9 ( 3 ): e1000602 . OpenUrl CrossRef PubMed ↵ Pinto J.D. 1991 . The taxonomy of North American Epicauta (Coleoptera: Meloidae), with a revision of the nominate subgenus and a survey of courtship behavior . University of California Press , California . 372 pp. ↵ Pinto J.D. & Bologna M.A . 1999 . The New World genera of Meloidae (Coleoptera): a key and synopsis . J. Nat. Hist . 33 : 569 – 620 . OpenUrl ↵ Pinto J.D. , Bologna M.A. & Bouseman J.K . 1996 . First instar larvae, courtship and oviposition in Eletica : amending the definition of the Meloidae (Coleoptera: Tenebrionoidea) . Syst. Entomol . 21 : 63 – 74 . OpenUrl ↵ Pinto J.D. & Selander R.B . 1970 . The bionomics of blister beetles of the genus Meloe and a classification of the New World species . Ill. Biol. Monogr . 42 : 1 – 222 . OpenUrl ↵ Poinar G.O . 2009 . Meloe dominicanus n. sp. (Coleoptera: Meloidae) phoretic on the bee Proplebia dominicana (Hymenoptera: Apidae) in Dominican amber . Proc. Entomol. Soc. Wash . 111 ( 1 ): 145 – 151 . OpenUrl ↵ Rabosky D.L . 2010 . Extinction rates should not be estimated from molecular phylogenies . Evolution 64 ( 6 ): 1816 – 1824 . OpenUrl CrossRef PubMed Web of Science Rabosky D.L . 2014 . Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees . PLoS One 9 : 389543 . OpenUrl ↵ Rabosky D.L. , Donnellan S.C. , Grundler M. & Lovette I.J . 2014 . Analysis and visualization of complex macroevolutionary dynamics: an example from Australian scincid lizards . Syst. Biol . 63 ( 4 ): 610 – 627 . OpenUrl CrossRef PubMed ↵ Rabosky D.L. & Goldberg E.E . 2015 . Model inadequacy and mistaken inferences of trait-dependent speciation . Syst. Biol . 64 ( 2 ): 340 – 355 . OpenUrl CrossRef PubMed ↵ Rabosky D.L. & Goldberg E.E . 2017 . FiSSE: A simple nonparametric test for the effects of a binary character on lineage diversification rates . Evolution 71 ( 6 ): 1432 – 1442 . OpenUrl CrossRef PubMed ↵ Rabosky D.L. , Santini F. , Eastman J. , Smith S.A. , Sidlauskas B. , Chang J. & Alfaro M.E . 2013 . Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation . Nat. Commun . 4 ( 1 ): 1 – 8 . OpenUrl CrossRef PubMed ↵ Ricklefs R.E . 2007 . Estimating diversification rates from phylogenetic information . Trends Ecol. Evol . 22 : 601 – 610 . OpenUrl CrossRef PubMed Web of Science ↵ Ricklefs R.E. & Fallon S.M . 2002 . Diversification and host switching in avian malaria parasites . P. Roy. Soc. B-Biol. Sci . 269 ( 1494 ): 885 – 892 . OpenUrl Rider Jr S.D . 2016 . The complete mitochondrial genome of the desert darkling beetle Asbolus verrucosus (Coleoptera , Tenebrionidae). Mitochondrial DNA A 27 ( 4 ): 2447 – 2449 . OpenUrl ↵ Ronquist F. , Teslenko M. , van der Mark P. , Ayres D.L. , Darling A. , Höhna S. , Larget B. , Liu L. , Suchard L.A. & Huelsenbeck J.P. 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Syst. Biol . 61 ( 3 ): 539 – 542 . OpenUrl CrossRef PubMed ↵ Runge F. & Thines M . 2012 . Reevaluation of the host specificity of the closely related species Pseudoperonospora humuli and P. cubensis . Plant Dis . 96 : 55 – 61 . OpenUrl ↵ Sanderson M.J. & Donoghue M.J . 1996 . Reconstructing shifts in diversification rates on phylogenetic trees . Trends Ecol. Evol . 11 ( 1 ): 15 – 20 . OpenUrl CrossRef PubMed Web of Science ↵ Sanmartín I. & Meseguer A.S . 2016 . Extinction in phylogenetics and biogeography: from timetrees to patterns of biotic assemblage . Front. Genet . 7 : 35 . OpenUrl ↵ Schlee D . 1990 . Das Bernstein-Kabinett. Begleitheft zur Bernsteinausstellung im Museum am LoCwentor, Stuttgart . Stuttg. Beitr. Naturkd, Ser. C 28 : 1 – 100 . OpenUrl ↵ Schulte R.D. , Makus C. , Hasert B. , Michiels N.K. & Schulenburg H . 2010 . Multiple reciprocal adaptations and rapid genetic change upon experimental coevolution of an animal host and its microbial parasite . Proc. Natl. Acad. Sci . 107 ( 16 ): 7359 – 7364 . OpenUrl Abstract / FREE Full Text ↵ Schwarz G . 1978 . Estimating the dimension of a model . Ann. Stat . 6 ( 2 ): 461 – 464 . OpenUrl CrossRef Web of Science ↵ Selander R.B . 1964 . Sexual Behavior in Blister Beetles (Coleoptera: Meloidae): I. The Genus Pyrota . Can. Entomol . 96 ( 8 ): 1037 – 1082 . OpenUrl ↵ Selander R.B. & Mathieu J.M . 1964 . The ontogeny of blister beetles (Coleoptera, Meloidae) I. A study of three species of the genus Pyrota . Ann. Entomol. Soc. Am . 57 ( 6 ): 711 – 732 . OpenUrl CrossRef ↵ Selander R.B. & Weddle R.C . 1969 . The ontogeny of blister beetles (Coleoptera, Meloidae). II. The effects of age of triungulin larvae at feeding and temperature on development in Epicauta segmenta . Ann. Entomol. Soc. Am . 62 ( 1 ): 27 – 39 . OpenUrl CrossRef ↵ Shi Y. , Wu Y. , Zhang W. , Qi J. & Gao G.F . 2014 . Enabling the “host jump”: structural determinants of receptor-binding specificity in influenza A viruses . Nat. Rev. Microbiol . 12 ( 12 ): 822 . OpenUrl CrossRef PubMed ↵ Silva D.N. , Talhinhas P. , Cai L. , Manuel L. , Gichuru E.K. , Loureiro A. , Várzea V. , Paulo O.S. & Batista D . 2012 . Host-jump drives rapid and recent ecological speciation of the emergent fungal pathogen Colletotrichum kahawae . Mol. Ecol . 21 ( 11 ): 2655 – 2670 . OpenUrl CrossRef PubMed ↵ Simões M. , Breitkreuz L. , Alvarado M. , Baca S. , Cooper J.C. , Heins L. , Herzog K. , & Lieberman B.S . 2016 . The Evolving Theory of Evolutionary Radiations . Trends. Ecol. Evol . 31 ( 1 ): 27 – 34 . OpenUrl CrossRef PubMed ↵ Song H. , Amédégnato C. , Cigliano M.M. , Desutter-Grandcolas L. , Heads S.W. , Huang Y. , Otte D. & Whiting M.F . 2015 . 300 million years of diversification: elucidating the patterns of orthopteran evolution based on comprehensive taxon and gene sampling . Cladistics 31 ( 6 ): 621 – 651 . OpenUrl CrossRef ↵ Sorenson M.D. , Sefc K.M. & Payne R.B . 2003 . Speciation by host switch in brood parasitic indigobirds . Nature 424 : 928 – 931 OpenUrl CrossRef PubMed Web of Science ↵ Stadler T . 2009 . On incomplete sampling under birth-death models and connections to the sampling-based coalescent . J. Theor. Biol . 261 ( 1 ): 58 – 66 . OpenUrl CrossRef PubMed Web of Science ↵ Stadler T . 2011 . Mammalian phylogeny reveals recent diversification rate shifts . Proc. Natl. Acad. Sci . 108 ( 15 ): 6187 – 6192 . OpenUrl Abstract / FREE Full Text ↵ Stamatakis A . 2006 . RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models . Bioinformatics 22 ( 21 ): 2688 – 2690 . OpenUrl CrossRef PubMed Web of Science ↵ Stanley S.M . 1975 . A theory of evolution above the species level . Proc. Natl. Acad. Sci . 72 ( 2 ): 646 – 650 . OpenUrl Abstract / FREE Full Text ↵ Talavera G. & Castresana J . 2007 . Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments . Syst. Biol . 56 : 564 – 577 . OpenUrl CrossRef PubMed Web of Science ↵ Tavaré S . 1986 . Some probabilistic and statistical problems on the analysis of DNA sequences . Lect. Math. Life Sci . 17 : 57 – 86 OpenUrl ↵ Thines M . 2019 . An evolutionary framework for host shifts–jumping ships for survival . New Phytol . 224 ( 2 ): 605 – 617 . OpenUrl CrossRef ↵ Thode V.A. , Lohmann L.G. & Sanmartín I . 2020 . Evaluating character partitioning and molecular models in plastid phylogenomics at low taxonomic levels: A case study using Amphilophium (Bignonieae , Bignoniaceae). J. Syst. Evol . doi: https://doi.org/10.1111/jse.12579 ↵ Timmermans M.J. , Barton C. , Haran J. , Ahrens D. , Culverwell C.L. , Ollikainen A. , Dodsworth S. , Foster P.G. , Bocak L. & Vogler A. P . 2015 . Family-level sampling of mitochondrial genomes in Coleoptera: compositional heterogeneity and phylogenetics . Genome Biol. Evol . 8 ( 1 ): 161 – 175 . OpenUrl PubMed Timmermans M.J. , Dodsworth S. , Culverwell C.L. , Bocak L. , Ahrens D. , Littlewood D.T. , Pons J. & Vogler A.P . 2010 . Why barcode? High-throughput multiplex sequencing of mitochondrial genomes for molecular systematics . Nucleic Acids Res . 38 ( 21 ): e197 – e197 . OpenUrl CrossRef PubMed ↵ Turrisi G.F. & Vilhelmsen L . 2010 . Into the wood and back: morphological adaptations to the wood-boring parasitoid lifestyle in adult aulacid wasps (Hymenoptera: Aulacidae) . J. Hym. Res . 19 ( 2 ): 244 – 258 . OpenUrl ↵ Uribe J.E. , Puillandre N. , Zardoya R . 2017 . Beyond Conus: phylogenetic relationships of Conidae based on complete mitochondrial genomes . Mol. Phylogenet. Evol . 107 : 142 – 51 . OpenUrl ↵ Vamosi J.C. , Otto S.P. & Barrett S.C . 2003 . Phylogenetic analysis of the ecological correlates of dioecy in angiosperms . J. Evol. Biol . 16 ( 5 ): 1006 – 1018 . OpenUrl CrossRef PubMed Web of Science ↵ Vermeij G.J . 2001 . Innovation and evolution at the edge: origins and fates of gastropods with a labral tooth . Biol. J. Linn. Soc . 72 ( 4 ): 461 – 508 . OpenUrl CrossRef GeoRef Web of Science ↵ Vermeij G.J . 2006 . Historical contingency and the purported uniqueness of evolutionary innovations . Proc. Natl. Acad. Sci . 103 ( 6 ): 1804 – 1809 . OpenUrl Abstract / FREE Full Text ↵ Villaverde T. , Pokorny L. , Olsson S. , Rincón-Barrado M. , Johnson M.G. , Gardner E.M. , Wickett N.J. , Molero J. , Rinna R. & Sanmartín I . 2018 . Bridging the micro-and macroevolutionary levels in phylogenomics: Hyb-Seq solves relationships from populations to species and above . New Phytol . 220 ( 2 ): 63 – 6650 . OpenUrl ↵ Whelan S. & Goldman N . 2001 . A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach . Mol. Biol. Evol . 18 ( 5 ): 691 – 699 . OpenUrl CrossRef PubMed Web of Science ↵ Wolfe N.D. , Dunavan C.P. & Diamond J . 2007 . Origins of major human infectious diseases . Nature 447 ( 7142 ): 279 – 283 . OpenUrl CrossRef PubMed Web of Science ↵ Xie W. , Lewis P.O. , Fan Y. , Kuo L. , Chen M.-H . 2011 . Improving marginal likelihood estimation for Bayesian phylogenetic model selection . Syst. Biol . 60 : 150 – 160 . OpenUrl CrossRef PubMed Web of Science ↵ Yan L. , Pape T. , Elgar M.A. , Gao Y. & Zhang D . 2019 . Evolutionary history of stomach bot flies in the light of mitogenomics . Syst. Entomol . 44 ( 4 ): 797 – 809 . OpenUrl ↵ Yang Z . 1994 . Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods . J. Mol. Evol . 39 ( 3 ): 306 – 314 . OpenUrl CrossRef PubMed Web of Science ↵ Young A.D. & Gillung J.P . 2020 . Phylogenomics–principles, opportunities and pitfalls of big-data phylogenetics . Syst. Entomol . 45 ( 2 ): 225 – 247 OpenUrl ↵ Yuan M.L. , Zhang Q.L. , Zhang L. , Guo Z.L. , Liu Y.J. , Shen Y.Y. & Shao R . 2016 . High-level phylogeny of the Coleoptera inferred with mitochondrial genome sequences . Mol. Phylogenet. Evol . 104 : 99 – 111 . OpenUrl ↵ Zhang T. , Wu Q. & Zhang Z . 2020 . Probable pangolin origin of SARS-CoV-2 associated with the COVID-19 outbreak . Curr. Biol . 30 ( 7 ): 1346 – 1351 . OpenUrl CrossRef PubMed ↵ Zhu F. , Xue F. & Lei C . 2006 . The effect of environmental conditions on diapause in the blister beetle, Mylabris phalerata (Coleoptera: Meloidae) . Eur. J. Entomol . 103 ( 3 ): 531 OpenUrl ↵ Zietara M.S. & Lumme J . 2002 . Speciation by host switch and adaptive radiation in a fish parasite genus Gyrodactylus (Monogenea, Gyrodactylidae) . Evolution 56 ( 12 ): 2445 – 2458 . OpenUrl CrossRef PubMed Web of Science Back to top Previous Next Posted January 04, 2021. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction? Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction? E.K. López-Estrada , I. Sanmartín , J.E. Uribe , S. Abalde , M. García-París bioRxiv 2021.01.04.425192; doi: https://doi.org/10.1101/2021.01.04.425192 Share This Article: Copy Citation Tools Diversification dynamics of hypermetamorphic blister beetles (Meloidae): Are homoplastic host shifts and phoresy key factors of a rushing forward strategy to escape extinction? E.K. López-Estrada , I. Sanmartín , J.E. Uribe , S. Abalde , M. García-París bioRxiv 2021.01.04.425192; doi: https://doi.org/10.1101/2021.01.04.425192 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7969) Biochemistry (18634) Bioengineering (14770) Bioinformatics (44159) Biophysics (22460) Cancer Biology (19594) Cell Biology (26747) Clinical Trials (138) Developmental Biology (13903) Ecology (20890) Epidemiology (2067) Evolutionary Biology (25317) Genetics (16100) Genomics (23402) Immunology (18609) Microbiology (42247) Molecular Biology (17948) Neuroscience (92896) Paleontology (693) Pathology (2969) Pharmacology and Toxicology (5063) Physiology (8066) Plant Biology (15909) Scientific Communication and Education (2091) Synthetic Biology (4538) Systems Biology (10188) Zoology (2376) window.__CF$cv$params={r:'a37c764faf144fa7',t:'MTc4ODg1NDg3Mw==',u:'01a0800f3b367e42a391697b7db441d6',ut:'ecKBK3OhDn3ntFYNkE7839b9F6tKFacwc_qUHBfi2Ww-1788854876-1.2.1.1-Gi6qv0HGvvFKdnau6xKtdez5IFqGdXE9JCCEaS.2tB9UaM5vCY5T0mwePJqaJgJoQ6IcE4XkG8dHPAePorTPcy8LLa4PqEWho8BGZiyUOpQ',i:60};(function(){if(!document.body)return;var s=document.createElement('script');s.src='/cdn-cgi/challenge-platform/scripts/precursor/main.js';document.head.appendChild(s);})();