Frédéric Legendre
email: frederic.legendre@mnhn.fr

Fabien L Condamine
email: fabien.condamine@gmail.com.

When Darwin's Special Difficulty Promotes Diversification in Insects

Keywords: Aggregate trait, emergent trait, eusocial insects, macroevolution, species selection, termites

Eusociality, Darwin's special difficulty, has been widely investigated but remains a topic of great debate in organismal biology. Eusocial species challenge existing theories, and the impact of highly integrated societies on diversification dynamics is controversial with opposing assertions and hypotheses in the literature. Here, using phylogenetic approaches in termites -the first group that has evolved eusociality -we assessed the fundamental prediction that eusocial lineages have higher diversification rates than non-eusocial clades.

We found multiple lines of evidence that eusociality provided higher diversification as compared to non-eusociality. This is particularly exacerbated for eusocial species with 'true' workers as compared to species with 'false' workers. Because most species with 'true' workers have an entirely prokaryotic microbiota, the latter feature is also related to higher diversification rates, but it should be investigated further, notably in relation to angiosperm diversification. Overall, this study suggests that societies with 'true' workers are not only more successful at ecological timescales but also over millions of years, which further implies that both organism-and species-level traits act on species selection.

Eusociality is the most integrated organizational level of animal sociality and is defined by cooperative brood care, overlapping generations within a colony of adults, and division of labour into reproductive and non-reproductive groups [START_REF] Wilson | Eusociality: origin and consequences[END_REF]. Although the origin of eusociality is still widely investigated and debated (e.g. [START_REF] Nowak | The evolution of eusociality[END_REF][START_REF] Nonacs | Kinship, greenbeards, and runaway social selection in the evolution of social insect cooperation[END_REF][START_REF] Rousset | Much ado about nothing: Nowak et al.'s charge against inclusive fitness theory[END_REF][START_REF] Johnstone | Sex-biased dispersal, haplodiploidy and the evolution of helping in social insects[END_REF], this rare but widespread behavioural characteristic has independently appeared in insects (at least [START_REF] Ma | Sponge@dwelling$ snapping$shrimps$$ Origineusociality$supportedmolecular[END_REF][8][9][10][11][12][13][14][15][16][17][18][19][20] times depending on the definition of eusociality), crustaceans (at least thrice), and mammals (at least twice) (Fig. 1; [START_REF] Bourke | Principles of Social Evolution[END_REF]. Within insects, eusocial species have colonies with caste differences: queens and reproductive males take the roles of sole reproducers, while soldiers defending the colony and workers foraging and maintaining resources contribute to creating a living situation favourable for the colony [START_REF] Wilson | The Insect Societies[END_REF]. In his theory of natural selection, [START_REF] Darwin | On the Origin of Species by Means of Natural Selection[END_REF] famously portrayed the eusocial insects as a "special difficulty, which at first appeared to [him] insuperable and actually fatal to [his] theory". Indeed, eusociality presents an apparent paradox: if adaptive evolution is mediated by differential reproduction of individuals on which natural selection acts, how can individuals incapable of passing on their genes evolve and persist? Figure 1. Origin and timeline of eusocial lineages for mammals, crustaceans and insects (bees, wasps, ants, and termites). Bars represent temporal ranges with solid lines based on the fossil record, whereas broken lines indicate possible eusociality supported by molecular dating studies (ants: [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF]bees: Cardinal and Danforth 2011, wasps: Hines et al. 2007, termites: Ware et al. 2010;Legendre et al. 2015b) but not supported by fossil evidence.

The multilevel selection theory could settle this paradox. This theory has a long and troubled history [START_REF] Okasha | The levels of selection debate: philosophical issues[END_REF]) but is increasingly accepted and supported, especially at the species level [START_REF] Jablonski | Species selection: theory and data[END_REF]. Species selection results from differences in speciation and extinction rates among species within a clade, and these differences are not merely related to the fitness at the level of individuals [START_REF] Jablonski | Species selection: theory and data[END_REF][START_REF] Rabosky | Reinventing species selection with molecular phylogenies[END_REF]. Data and theory suggest that both organism-level and species-level traits (called aggregate and emergent traits, respectively) act on species selection. Aggregate traits are characteristics of individual organisms (e.g. body size), whereas emergent traits are species properties that are not reducible to organismal traits such as eusociality. Extending this vocabulary, authors often refer to the emergent fitness of species within clades when discussing species selection [START_REF] Jablonski | Species selection: theory and data[END_REF].

At the macroevolutionary timescale, the origin and evolution of eusociality are increasingly well-understood due to recent phylogenetic analyses, along with estimates of divergence times and diversification rates notably for arthropod groups (ants: [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF][START_REF] Pie | Relictual ant lineages (Hymenoptera: Formicidae) and their evolutionary implications[END_REF], bees: Cardinal and Danforth 2011, beetles: Farrell et al. 2001, wasps: Hines et al. 2007, spiders: [START_REF] Agnarsson | Sociality in theridiid spiders: repeated origins of an evolutionary dead end[END_REF], termites: Ware et al. 2010). These studies have investigated the number of origins of eusociality and the putative factors behind the extraordinary diversification of eusocial groups (e.g. the Cretaceous radiation of Angiosperms: [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF][START_REF] Cardinal | Bees diversified in the age of eudicots[END_REF]. But another challenge surrounds eusocial organism evolution: the impact of highly integrated societies on the dynamics of diversification in these organisms [START_REF] Davis | Eusociality and the success of the termites: insights from a supertree of dictyopteran families[END_REF][START_REF] Ware | The effects of fossil placement and calibration on divergence times and rates: an example from the termites (Insecta: Isoptera)[END_REF]. This topic has received less attention probably due to a lack of methods to handle large phylogenies that include both non-eusocial and eusocial lineages. Although the role of eusociality has not been formally studied with state-of-the-art diversification methods, the literature is filled with opposing assertions and hypotheses. In terms of ecological and evolutionary success, modern eusocial insects can either be seen as species-rich groups such as ants, bees, or termites [START_REF] Chapman | The influence of sociality on the conservation biology of social insects[END_REF][START_REF] Wilson | Success and Dominance in Ecosystems: the Case of the Social Insect[END_REF] or as a group with lower species diversity than non-eusocial relatives such as ambrosia beetles, or poplar spiral gall aphids [START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF]. Although the success of eusocial groups could naively lead us to think that eusociality was a trigger of diversification, [START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF] argued that social insects (contrary to social vertebrates) show lower diversification (because they are "less speciose") than non-eusocial insects. Based on studies of population genetics and molecular evolution, sociality reduces effective population size and possibly reduces levels of genetic variation, which would in turn increase the risk of extinction of social organisms [START_REF] Lanfear | Population size and the rate of evolution[END_REF].

Hence, eusociality could indeed be seen as an impediment to diversification. The fossil record shows that taxa went extinct within eusocial groups (ants: LaPolla et al. 2013, termites: Krishna et al. 2010) but, so far, eusocial insect groups have been successful at the macroevolutionary scale. For the most part, they have persisted since the origins of these highly integrated behaviours, leading [START_REF] Wilson | Success and Dominance in Ecosystems: the Case of the Social Insect[END_REF] to state that "eusociality conveys evolutionary long life to social insects". In reality, several emergent traits influence both extinction and speciation rates, resulting in unpredictable trade-offs [START_REF] Jablonski | Species selection: theory and data[END_REF].

The ensuing questions 'does eusociality confer higher diversification rates than noneusocial species?' and 'what is the dynamic of diversification of social insects and what factor could explain their success and potential shift in diversification rates?' are long-standing research topics [START_REF] Darwin | On the Origin of Species by Means of Natural Selection[END_REF][START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF][START_REF] Wilson | Eusociality: origin and consequences[END_REF]. Biotic and abiotic factors are important in diversification processes [START_REF] Benton | The Red Queen and the Court Jester: species diversity and the role of biotic and abiotic factors through time[END_REF] and range from climatic/geologic events (e.g. [START_REF] Condamine | Macroevolutionary perspectives to environmental change[END_REF]), to competition (e.g. Liow et al. 2015;Silvestro et al. 2015a), and key innovations (e.g. [START_REF] Rainford | Phylogenetic distribution of extant richness suggests metamorphosis is a key innovation driving diversification in insects[END_REF][START_REF] Sánchez-García | Is the switch to an ectomycorrhizal state an evolutionary key innovation in mushroom-forming fungi? A case study in the Tricholomatineae (Agaricales)[END_REF]. Eusociality can be classified as a key innovation, but has surprisingly been understudied when one considers the ecological success that eusocial insect species have experienced [START_REF] Wilson | The Insect Societies[END_REF][START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF][START_REF] Chapman | The influence of sociality on the conservation biology of social insects[END_REF]. This ecological dominance is mainly explained by one property of eusociality: the reproductive division of labour. In other words, some individuals do not reproduce (i.e. they belong to sterile castes), and this has been qualified as the most important feature of eusociality [START_REF] Wilson | Success and Dominance in Ecosystems: the Case of the Social Insect[END_REF].

Within insects, termites (Dictyoptera, order Blattodea, infraorder Isoptera) contain only eusocial species, and it is thought that they were the first modern eusocial animals to evolve, sometime in the Late Jurassic (ca. 155 Ma) (Legendre et al. 2015b;[START_REF] Engel | Morphologically specialized termite castes and advanced sociality in the Early Cretaceous[END_REF].

A termite colony is usually differentiated into reproductive, worker and soldier castes (Fig. 2a). The latter two castes, which are sterile and perform highly specialized tasks, derive from different stages depending on the species [START_REF] Noirot | Ontogenetic development and evolution of the worker caste in termites[END_REF]. Whereas soldiers likely have a single origin, the evolution of workers is more controversial (Watson and Sewell 1981;[START_REF] Noirot | Ontogenetic development and evolution of the worker caste in termites[END_REF][START_REF] Thompson | Phylogenetic evidence for a single, ancestral origin of a "true" worker caste in termites[END_REF]Grandcolas and D'Haese 2002;[START_REF] Inward | A comprehensive phylogenetic analysis of termites (Isoptera) illuminates key aspects of their evolutionary biology[END_REF][START_REF] Legendre | The phylogeny of termites (Dictyoptera: Isoptera) based on mitochondrial and nuclear markers: Implications for the evolution of the worker and pseudergate castes, and foraging behaviors[END_REF], 2013;[START_REF] Roisin | Social Organisation and the Status of Workers in Termites[END_REF]. Importantly, two types of workers were distinguished early: Some species have a 'true' worker caste, while others have a 'false' worker caste (also called pseudergates; Fig. 2a). 'True' workers differ from pseudergates notably because they diverge early and irreversibly from the imaginal line [START_REF] Noirot | Ontogenetic development and evolution of the worker caste in termites[END_REF], while pseudergates do not. Societies with 'true' workers have been qualified as more 'socially' highly integrated, which might have consequences for their ecological and evolutionary success given that >80% of all termites have 'true' workers [START_REF] Korb | Termites: an alternative road to eusociality and the importance of group benefits in social insects[END_REF][START_REF] Brune | Symbiotic digestion of lignocellulose in termite guts[END_REF]. illustrates that the origin of eusociality is nested within the radiation of cockroaches (termites and cockroaches together forming the Blattodea). The phylogeny also highlights important events in the evolution of the digestive symbiosis: the presence of cellulolytic flagellates, acquired by a common ancestor of termites and Cryptocercidae, is associated to their woodfeeding lifestyle; the loss of flagellates in Termitidae gave rise to an enormous dietary diversification (chronogram adapted from Legendre et al. 2015b). Dashed lines for the common ancestor of cockroaches denote that the group is older than 201 Ma. The asterisk (*) specifies that some families or subfamilies are paraphyletic according to the latest studies.

The study of termite digestive systems might also bring important insights on the ecological and evolutionary successes of these insects (Fig. 2b). Although digestive strategies and the gut microbiota of termites remain to be fully deciphered, important progress has been made and two groups can be distinguished by their primary cellulolytic partners and their hindgut compartmentation [START_REF] Brune | Symbiotic digestion of lignocellulose in termite guts[END_REF]. The first group comprises termite species that harbour cellulolytic flagellates in their gut; they group all termite species except those belonging to Termitidae. In the second group (i.e. Termitidae), these flagellates are absent; instead Termitidae have an entirely prokaryotic microbiota. They also show a higher compartmentation of the hindgut (except for Macrotermitinae) and a large dietary diversification [START_REF] Brune | The gut microbiota of termites: Digesting the diversity in the light of ecology and evolution[END_REF]. Because most termite species belong to the family Termitidae [START_REF] Beccaloni | Order Blattodea[END_REF], their peculiar digestive system might have been a trigger of diversification. A parallel can be made with cockroach species; even though the gut microbiota of most cockroach species has not been investigated, flagellates have been found in some wood-eating species, including, but not restricted to (e.g. [START_REF] Pellens | A new and independently evolved case of xylophagy and the presence of intestinal flagellates in the cockroach Parasphaeria boleiriana (Dictyoptera, Blaberidae, Zetoborinae) from the remnants of the Brazilian Atlantic forest[END_REF], Cryptocercus, the alleged sister group of termites [START_REF] Cleveland | The wood feeding roach Cryptocercus, its protozoa, and the symbiosis between protozoa and roach[END_REF][START_REF] Klass | Wood-feeding cockroaches as models for termite evolution (Insecta: Dictyoptera): Cryptocercus vs. Parasphaeria boleiriana[END_REF].

Because of characteristics shared with 'lower' termites [START_REF] Klass | Wood-feeding cockroaches as models for termite evolution (Insecta: Dictyoptera): Cryptocercus vs. Parasphaeria boleiriana[END_REF], it is then generally assumed that Cryptocercus belongs to the group with cellulolytic flagellates in their gut [START_REF] Brune | The gut microbiota of termites: Digesting the diversity in the light of ecology and evolution[END_REF], contrary to other supposedly or asserted wood-eating species (but see [START_REF] Pellens | The evolutionary transition from subsocial to eusocial behaviour in Dictyoptera: Phylogenetic evidence for modification of the "shift-in-dependent-care" hypothesis with a new subsocial cockroach[END_REF] and the section Trait-dependent diversification for more details and alternative codings).

Our knowledge on the role of eusociality in the dynamic of diversification is limited.

Only one previous study has investigated the diversification of termites [START_REF] Davis | Eusociality and the success of the termites: insights from a supertree of dictyopteran families[END_REF], suggesting increases in diversification early in termite evolution. However, this study was based on a family-level supertree, and methods to estimate diversification were still in their infancy. Besides, this study did not differentiate the potential role of 'true' workers or that of gut microbiota composition in termite diversification. To better understand the triggers of diversification in Dictyoptera (mantises, cockroaches and termites), we analyze the most comprehensive time-calibrated phylogeny of the group. We investigate the diversification dynamic of termites compared to their dictyopteran counterparts using a battery of birth-death methods to corroborate and strengthen these results. More specifically, we tested the competing hypotheses in which eusociality is associated with lower [START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF] or higher [START_REF] Davis | Eusociality and the success of the termites: insights from a supertree of dictyopteran families[END_REF]) net diversification rate compared to non-eusociality. We also tested whether the presence of 'true' workers in termite societies can be considered a key innovation relative to societies with pseudergates, by assessing the putative role on diversification of these castes on termite evolution. Similarly, we tested whether an entirely prokaryotic gut microbiota could also be a key innovation for termite diversification.

MATERIALS AND METHODS

Analytical Pipeline

We used a recently published time-calibrated phylogeny of Dictyoptera including 762 species representatives of living dictyopteran diversity (ca. 10,000 known species, Legendre et al. 2015b). This phylogeny resulted from a molecular dataset comprising four mitochondrial and two nuclear markers. Divergence times were estimated using a relaxedclock approach coupled with 17 fossils, which were used as minimum age constraints on different nodes and a maximum age for the tree root (Legendre et al. 2015b).

Three putatively key innovations were investigated: eusociality, societies with 'true' workers and entirely prokaryotic gut microbiota. To test the role of these innovations, we used three approaches: (i) the maximum likelihood (ML) approach of time-dependent diversification [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF], implemented in the R-package RPANDA v.1.2 [START_REF] Morlon | RPANDA: An R package for macroevolutionary analyses on phylogenetic trees[END_REF], (ii) the Bayesian analysis of macroevolutionary mixture (BAMM v.2.5, Rabosky et al. 2013), implemented in the R-package BAMMtools v.2.1.4 (Rabosky et al. 2014), and (iii) the ML approach of trait-dependent diversification (Madisson et al. 2007;[START_REF] Fitzjohn | Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies[END_REF], implemented in the R-package diversitree v.0.9-8 (FitzJohn 2012).

Each method is designed to estimate speciation and extinction rates, and the three methods are used to cross-test hypotheses and corroborate results. Nonetheless, it is worth mentioning that each method differs at several points in the way speciation and extinction rates are estimated. For instance, trait-dependent birth-death models (diversitree) estimate constant speciation and extinction rates, while time-dependent birth-death models (BAMM and RPANDA) estimate the speciation and extinction rates and their variation through time.

Therefore we expect some differences in the values of the estimated diversification rates that are inherent to each approach. Our diversification analyses should be seen as complementary for the inferred diversification trend rather than corroborative on the values and magnitude of the speciation and extinction rates.

Across-Clade and Time-Variation Diversification

The ML approach of [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF] is a birth-death method that extends previous birth-death methods such that speciation and/or extinction rates may change continuously through time, and subclades may have different speciation and extinction rates. This method has the advantage of not assuming constant extinction rate over time (unlike BAMM, Rabosky et al. 2013), and allows clades to have declining diversity because extinction can exceed speciation, meaning that diversification rates can be negative [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF].

We designed six nested diversification models to test with this approach: (i) a Yule First, the Dictyoptera tree was analyzed as a whole using this approach and taking into account missing species (f=0.08). To test the hypotheses that eusociality, societies with 'true' workers and gut microbiota composition act on diversification rates, we compared the likelihoods of models that allow for different patterns of rate variation in different clades.

More specifically, we first allowed for a rate shift along the branch corresponding to the apparition of eusociality (the termite crown node), and we sequentially tested for a rate shift between the crown node of termites and the node corresponding to the family Termitidae (prokaryotic microbiota). Therefore, we partitioned the Dictyoptera phylogeny into seven evolutionary scenarios: from a backbone without the termites, and the termites as a subtree; to a backbone without the Termitidae, and the Termitidae as a subtree (with all five possibilities in between; for a schematic illustration see Fig. S1 in the Supplementary Material that accompanies this article available on Dryad at http://dx.doi.org/10.5061/dryad.2tg34). After creating these backbones and corresponding subtrees, and accounting for the missing species in each of them, we fitted the same diversification models as detailed above but independently on each subtree and its corresponding backbone tree (Fig. S1 on Dryad). The global loglikelihood was calculated as the sum of both subclade and backbone log-likelihoods as determined by the corresponding best-fitting models. This approach allowed us to compare the global log-likelihood (under the assumption of a single speciation rate and extinction rate) to the log-likelihood of one shift at the termite crown (eusociality), and to the log-likelihood of one shift at the crown of Termitidae (prokaryotic microbiota). We directly compare the seven scenarios using AICc.

We also used BAMM to estimate speciation and extinction rates through time and among/within clades (Rabosky et al. 2013). BAMM was constructed to study complex evolutionary processes on phylogenetic trees, potentially shaped by a heterogeneous mixture of distinct evolutionary dynamics of speciation and extinction across clades. BAMM can automatically detect rate shifts and sample distinct evolutionary dynamics that explain the diversification dynamics of a clade without a priori hypotheses on how many and where these shifts might occur. Evolutionary dynamics can involve time-variable diversification rates; in BAMM, speciation is allowed to vary exponentially through time while extinction is maintained constant: subclades in a tree may diversify faster (or slower) than others. This Bayesian approach can be useful in detecting shifts of diversification potentially associated with key innovations [START_REF] Rabosky | Automatic detection of key innovations, rate shifts, and diversitydependence on phylogenetic trees[END_REF].

We ran BAMM by setting four Markov chain Monte Carlo running (MCMC) for 20 million generations and sampled every 2,000 generations. A compound Poisson process is implemented for the prior probability of a rate shift along any branch. We used a gradient of prior values ranging from 0.1 to 50 to test the sensitivity to the prior, because it has been shown that BAMM can be affected by the prior [START_REF] Moore | Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures[END_REF] but see [START_REF] Rabosky | Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models[END_REF]. We accounted for non-random incomplete taxon sampling using the implemented analytical correction; we set a sampling fraction for mantises (f=0.118), cockroaches (f=0.042), and termites (f=0.094). We performed independent runs (with a 15% burn-in) using different seeds to assess the convergence of the runs with effective sample size. We processed the output data using the BAMMtools [START_REF] Rabosky | BAMMtools: An R package for the analysis of evolutionary dynamics on phylogenetic trees[END_REF]) by estimating (i) the mean global rates of diversification through time, (ii) the configuration of the diversification rate shifts evaluating alternative diversification models as compared by posterior probabilities, and (iii) the clade-specific rates through time when a distinct macroevolutionary regime is identified.

Trait-Dependent Diversification

We used trait-dependent diversification models to simultaneously model trait evolution and its impact on diversification [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF]). In those models, species are characterized by an evolving trait, and their diversification follows a birth-death process in which speciation and extinction rates may depend on the trait state. We created three datasets of traits. First, we categorized all dictyopteran species as being non-eusocial or eusocial. This two-trait scheme distinguishes the termites (all species are eusocial) from the rest of dictyopterans (all non-eusocial). Second, we made a three-trait scheme by keeping the noneusocial category, and dividing the eusocial trait into two: eusocial with pseudergates ('false' workers) and eusocial with 'true' workers. Third, we made another three-trait scheme to test the role of the primary cellulolytic partners in diversification. We followed [START_REF] Brune | The gut microbiota of termites: Digesting the diversity in the light of ecology and evolution[END_REF] to distinguish Termitidae (entirely prokaryotic microbiota), Cryptocercus + other termites (cellulolytic flagellates), and the remaining Dictyoptera (no specialized microbiota for lignocellulose digestion). All these coding schemes are imperfect because, like any broad categorization process, they might oversimplify the reality found in nature (Legendre et al. 2015a;[START_REF] Goutte | Environmental constraints and call evolution in torrent dwelling frogs[END_REF]. The category 'remaining Dictyoptera' in the latter coding scheme, for instance, is imperfect because it is delineated by default and because of a lack of information for some cockroach species. However, this three-state coding reflects our current understanding of the primary cellulolytic partners in Dictyoptera as reviewed in [START_REF] Brune | The gut microbiota of termites: Digesting the diversity in the light of ecology and evolution[END_REF]. Anyway, to take into account our lack of knowledge on cockroach microbiota and the existence of species that, arguably, could not fit well in these three categories -while limiting the number of categories to save statistical power, we created two supplementary trait datasets. In the first alternative coding, the wood-feeding non-eusocial lineages were coded with cellulolytic flagellates (same as 'lower' termites) even though the existence and role of cellulolytic flagellates have not always been investigated and proved.

These lineages represent several genera of the subfamily Panesthiinae, and the genera Cyrtotria, Lauraesilpha, Paramuzoa, Parasphaeria, and Colapteroblatta [START_REF] Grandcolas | Le genre Paramuzoa Roth, 1973 : sa répartition et un cas de xylophagie chez les Nyctiborinae (Dictyoptera, Blattaria)[END_REF]Grandcolas et al. 2002;[START_REF] Pellens | A new and independently evolved case of xylophagy and the presence of intestinal flagellates in the cockroach Parasphaeria boleiriana (Dictyoptera, Blaberidae, Zetoborinae) from the remnants of the Brazilian Atlantic forest[END_REF]. In the second alternative coding, the Macrotermitinae were not coded with the Termitidae (i.e. initial coding supported by their entirely prokaryotic gut microbiota) but with the other termites (coding supported by their relatively simple hindgut structure).

We applied the Binary State Speciation and Extinction model (BiSSE, [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF]) with constant diversification rate for the two-trait dataset. We then performed the Multi-State Speciation and Extinction model (MuSSE, FitzJohn et al. 2009) on the three-trait datasets. Both BiSSE and MuSSE models account for incomplete taxon sampling, which is informed as a sampling fraction of species at present having a given trait [START_REF] Fitzjohn | Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies[END_REF]. Simulation studies have shown that a large, well-sampled tree is required by these methods, whereas trees containing fewer than 300 species may lack sufficient phylogenetic signal to produce enough statistical power [START_REF] Davis | Exploring power and parameter estimation of the BiSSE method for analyzing species diversification[END_REF]).

The BiSSE model has six distinct parameters: two speciation rates without character change (i.e. in situ speciation) associated with non-eusocial species (λ N ) and eusocial species (λ E ), two extinction rates associated with non-eusocial (µ N ) and eusocial species (µ E ), and two transition rates (i.e. anagenetic change) with one from non-eusocial to eusocial (q N-E ), and from eusocial to non-eusocial (q E-N ). Using the same rationale, but applied to three traits, the MuSSE model has 12 distinct parameters (three for speciation, three for extinction, and six for transitions).

Here we used 100 dictyopteran phylogenies (sampled from the Bayesian analysis in Legendre et al. 2015b) and combined them with our three datasets of traits. Analyses were performed using the R-package diversitree (FitzJohn 2012) using the functions make.bisse and make.musse to construct the likelihood functions for each model based on the data, and the functions constrain and find.mle to apply different diversification scenarios.

For each model, we computed the AICc based on the log-likelihood and the number of parameters. We checked support for the selected model against all models using the difference between AICc (∆AIC) and the Akaike weight (AICω). The scenario supported with the lowest AICc was considered the best when ∆AIC>2 against other models and with AICω>0.5 (otherwise the model with less parameter was instead considered the best). Finally, we used the consensus tree and MCMC simulations with the best model to examine the confidence interval of the parameter estimates. We used an exponential prior and started the chain with the parameters obtained by maximum likelihood (FitzJohn 2012). We ran 20,000 MCMC steps and applied a 10% burn-in. Net diversification rates were then computed. SSE models have recently been criticized due to type I error [START_REF] Rabosky | Model inadequacy and mistaken inferences of traitdependent speciation[END_REF]. To test whether our diversification results were biased, we estimated the difference of fit (ΔAIC) between the best model and a null model (i.e. with no state dependence) and compared this with the difference between the same models as estimated from 100 simulated trait datasets for each SSE model.

RESULTS

Diversification Across-Clade and Time

The BAMM analyses supported a model with six evolutionary regimes (i.e. five rate shifts; Fig. S2 and Table S1 available on Dryad). The post burn-in posterior distribution indicated that five shifts occur with a posterior probability (PP) of 0.54, and PP=0.33 and 0.10 for four and six shifts, respectively (less or more shifts occur with PP<0.05; Table S1 available on Dryad). The 95% credible set of shift configurations and marginal shift probabilities strongly favoured a model that includes shifts along branches within termites (Figs S3-S4 available on Dryad). The shift configurations within this credible set contained also shifts at internal nodes in mantises and at a more derived position within cockroaches (Blaberidae: Panesthiinae). The different shift configurations in the credible set allowed quantifying uncertainty in placement of a termite shift (Fig. S4 available on Dryad). The shift configurations sampled at the highest frequency always contained a shift within termites (cumulative posterior probability of 0.62). The best configuration shift retained five core shifts, one located within the termites (crown Termitidae), two within Mantodea and two others within cockroaches (Fig. S5 available on Dryad). We repeated the analyses by changing the Poisson process for the prior probability of a rate shift along any branch: these analyses supported a similar diversification pattern, although we found an additional cladespecific macroevolutionary regime when a higher prior probability of rate shift was used (Fig.

S6 and Table S1 available on Dryad).

We found evidence for low and constant-rate diversification through time for Dictyoptera punctuated by few shifts. When diversification changed, speciation rates increased for the termites in particular at the base of Termitidae (Fig. 3). Although other increases of speciation were found, the most elevated speciation rate (and also net diversification rate) occurred within termites. Interestingly, when termite macroevolutionary regime was studied in isolation, the net diversification rate of the clade increased towards the present (Fig. 3). On the contrary, the other macroevolutionary regimes showed a more common pattern of diversification slowdown over time.

Termite Key Innovations and Diversification

The RPANDA analyses refined the location of the rate shift within termites, and allowed testing the effect of three termite key innovations: eusociality, societies with 'true' workers, and entirely prokaryotic gut microbiota, which have evolved sequentially from the origin of termites to the origin of Termitidae. Applying a series of time-dependent birth-death models (including a rate shift at the crown node of a termite subclade), we compared the fit of alternative scenarios to locate the best change of diversification rates. We first found that any scenario including a rate shift within termites is better than a scenario having no shift at all (Table 1). Comparing all diversification scenarios including a shift in the termite tree, we found that the best-fitting scenario includes a shift along a branch supporting termites that only contain species with 'true' workers in their societies (scenario 6, Table 1). There is strong support for this scenario (lowest ∆AIC=12.7) as compared to competing hypotheses of other locations in the termite tree such as the evolution of eusociality (scenario 1), or the change in gut microbiota composition (scenario 7).

BAMM identified an important rate shift within the mantises (166.8 Ma, Fig. 3), which corresponds to an increase of diversification. Using RPANDA, we tested whether a model including this shift performs better than the model with the best-fitting shift in termites.

We reproduced the approach explained above but we isolated the mantis clade as a subtree, from the rest of the Dictyoptera (this time including termites). We analyzed all the same models and found that the shift within the Mantodea does not improve the likelihood as compared with the model including the shift within termites (AICc for the shift in mantises=7640.83; resulting in a ∆AIC=327 between the best scenario including a shift within termites and this scenario, Table 1).

RPANDA analyses endorsed BAMM results: The best-fit scenario indicated a speciation rate increasing through time for both the termite tree and the backbone.

Nonetheless, we estimated an elevated speciation rate of 0.193 lineage/Myr for termites, strikingly contrasting with the speciation rate of the remaining dictyopteran lineages (0.058 lineage/Myr). Moreover, we found the termite clade diversified with low or zero extinction rates, again contrasting with the extinction rate for the rest of the tree (0.023 lineage/Myr).

Consequently, the net diversification rate of the termite tree is 5.5 times higher than the remaining part of the dictyopteran tree (scenario 7, Table 1).

We corroborated the BAMM and RPANDA analyses with BiSSE and MuSSE applied to a two-trait dataset and two three-trait datasets, respectively. We found that the best-fitting BiSSE model was the one with different speciation and extinction, but with equal transition rates. Simpler or more complex models were not supported (∆AIC=4.13 with the second bestfit BiSSE model, which has one extra parameter, Table 2a, Table S2 available on Dryad). In the best-fit model, both speciation and extinction rates were higher for eusocial species (Table 2a), and the net diversification rate was two-fold higher for eusocial lineages (Fig. 3).

Bayesian MCMC analysis showed that speciation, extinction, and net diversification rates were significantly different between the two traits (Fig. S7 available on Dryad).

We further investigated whether difference of diversification rates occurred within the termites by splitting the eusociality trait into two (societies with pseudergates or societies with 'true' workers). We found that the best-fitting MuSSE model supported the most complex model, in which all parameters are free. Simpler models were not supported (∆AICc=15.8 with the second best-fit MuSSE model, which has one less parameter, Table 2b, Table S2 available on Dryad). In the best-fit model, speciation and extinction rates were higher for eusocial lineages with 'true' workers. Speciation and extinction rates were both low for the two other traits (non-eusocial, and eusocial with pseudergates). The net diversification rate was higher for eusocial lineages with 'true' workers (Fig. 3). The trait-dependent rates of speciation, extinction, and net diversification for the eusocial species with 'true' workers were significantly different from the rates of the two other traits, and the two other traits were not significantly different from each other (Fig. S8 available on Dryad).

As for the primary cellulolytic partners, we found that the net diversification rate was higher for the Termitidae, the lineage with an entirely prokaryotic gut microbiota (Fig. 3; Table 2c; Fig. S9 and Table S2 available on Dryad). The two other traits were not significantly different in terms of diversification even though the lineages with cellulolytic flagellates showed higher speciation and extinction rates than mantises and cockroaches (except Cryptocercus that harbours flagellates).

Randomization analyses of the three datasets of traits showed that all the SSE-based results were robust to type-I error (Fig. S10 available on Dryad) indicating that our inferences are not biased. We performed additional MuSSE analyses by re-coding either the woodfeeding non-eusocial lineages or Macrotermitinae with Cryptocercus spp. and all the termites but the other Termitidae, and we found that the results were not sensitive to this effect: eusocial lineages with a prokaryotic microbiota or with a complex hindgut compartmentation diversified significantly faster than their counterparts (Figs. S11-S12 available on Dryad).

DISCUSSION

The origin of eusociality has been widely investigated at the organismal level. For instance, [START_REF] Hall | The effects of kin selection on rates of molecular evolution in social insects[END_REF] explored the differences in rate of substitution for queen-and worker-selected loci, whereas [START_REF] Rehan | Climbing the social ladder: the molecular evolution of sociality[END_REF] reviewed several hypotheses about molecular evolution of eusociality and tried to integrate them in a synthetic framework. Above the species level, the origin and evolution of eusociality has been also addressed with molecular dated phylogenies of eusocial groups [START_REF] Farrell | The evolution of agriculture in beetles (Curculionidae: Scolytinae and Platypodinae)[END_REF][START_REF] Agnarsson | Sociality in theridiid spiders: repeated origins of an evolutionary dead end[END_REF][START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF][START_REF] Hines | Multigene phylogeny reveals eusociality evolved twice in vespid wasps[END_REF]Cardinal andDanforth 2011, 2013;[START_REF] Pie | Relictual ant lineages (Hymenoptera: Formicidae) and their evolutionary implications[END_REF]. These studies have shown that eusociality is a labile trait that has appeared independently in many clades. However, the consequences of eusociality as a driver of species diversification have received less attention [START_REF] Davis | Eusociality and the success of the termites: insights from a supertree of dictyopteran families[END_REF][START_REF] Ware | The effects of fossil placement and calibration on divergence times and rates: an example from the termites (Insecta: Isoptera)[END_REF].

What role does eusociality play at this macroevolutionary scale? In social insects, only two studies mentioned this aspect and suggested contradictory hypotheses. Edward O. [START_REF] Wilson | The effects of complex social life on evolution and biodiversity[END_REF] proposed that social insects would have a lower diversification rate than non-eusocial insects, while [START_REF] Davis | Eusociality and the success of the termites: insights from a supertree of dictyopteran families[END_REF] suggested that termites would have a higher diversification rate than their close relatives; the study of [START_REF] Ware | The effects of fossil placement and calibration on divergence times and rates: an example from the termites (Insecta: Isoptera)[END_REF] did not compare termite diversification with the non-eusocial clades. These opposite predictions are in fact not so surprising because, for numerous aggregate and emergent traits, both low speciation and extinction rates are predicted [START_REF] Jablonski | Species selection: theory and data[END_REF]). These traits would then lead either to increase or decrease in diversification depending on the intensity of speciation and extinction.

Speciation rate supposedly increases with molecular rate [START_REF] Webster | Molecular phylogenies link rates of evolution and speciation[END_REF][START_REF] Lanfear | Mutation rate is linked to diversification in birds[END_REF]; but see [START_REF] Rabosky | Macroevolutionary speciation rates are decoupled from the evolution of intrinsic reproductive isolation in Drosophila and birds[END_REF], which in turn depends on several features including population size, population abundance, and genetic population. [START_REF] Harvey | Positive association between population genetic differentiation and speciation rates in New World birds[END_REF] showed, for instance, a link between population differentiation within species and speciation rates inferred from phylogenies in New World birds. All the aforementioned features are emergent traits connected to eusociality. In termites, the archetype society mode involves large and abundant populations headed by pairs of reproductives, resulting in colonies with strong genetic structures and genetically isolated by distance [START_REF] Goodisman | Population and colony genetic structure of the primitive termite Mastotermes darwiniensis[END_REF][START_REF] Thompson | Molecular-genetic analyses of dispersal and breeding behaviour in the Australian termite Coptotermes lacteus: Evidence for non-random mating in a swarm-dispersal mating system[END_REF][START_REF] Dupont | Genetic variability of the soil-feeding termite Labiotermes labralis (Termitidae, Nasutitermitinae) in the Amazonian primary forest and remnant patches[END_REF][START_REF] Eggleton | An introduction to termites: biology, taxonomy and functional morphology[END_REF]Vargo and Hussedener 2011).

According to classical predictions, high abundance would imply low extinction and low speciation rates, while highly structured populations would imply both high extinction and speciation rates. As for large population size, predictions suggest low extinction rates and are contradictory about speciation rates [START_REF] Jablonski | Species selection: theory and data[END_REF]. For Dictyoptera, we found here higher extinction and speciation rates for eusocial lineages when compared to non-eusocial lineages (using three different approaches). Because speciation is much higher than extinction, the net diversification rate is higher for eusocial species than for non-eusocial species, resulting in a strong positive link between eusociality and net diversification.

These results suggest that a high level of genetic structure could be the main factor acting on Dictyoptera diversification. However, the classical vision of termites is outdated.

Population structure in termites probably shows a large spectrum of variation that should be investigated further [START_REF] Vargo | Genetic structure of termite colonies and populations[END_REF]) and compared to cockroach population structures that should also vary between solitary, gregarious and subsocial species [START_REF] Bell | Cockroaches: Ecology, Behavior, and Natural History[END_REF], or flying and wingless species. Population size -and size-dependent properties -is also poorly known [START_REF] Korb | Termites: an alternative road to eusociality and the importance of group benefits in social insects[END_REF]) but range from a few hundreds to millions of individuals.

Even within the family Termitidae comprising only species with 'true' workers, there is a large variance with nest population of a few thousands to a few millions of individuals [START_REF] Lepage | Population dynamics of termites[END_REF]. But, even though the underlying factors remain to be investigated further, eusociality seems to act as a key innovation favouring diversification in Dictyoptera.

Key innovations are alluded to whenever an evolutionary novelty is thought to have favoured the diversification of a clade. Three main categories of key innovations can be distinguished [START_REF] Heard | Key evolutionary innovations and their ecological mechanisms[END_REF]: (i) innovations opening new adaptive zones sensu [START_REF] Simpson | Tempo and Mode in Evolution[END_REF]; (ii) innovations increasing fitness and allowing to outcompete other species; and (iii) innovations increasing the propensity for reproductive or ecological specialization. These three categories are not necessarily exclusive. For termites, eusociality unlikely opened new adaptive zones (definition 1). Termites act mainly as decomposers, a niche already existing and occupied before the emergence of this clade [START_REF] Raymond | Rates and processes of terrestrial nutrient cycling in the Paleozoic: the world before beetles, termites, and flies[END_REF]. However, a highly integrated society, especially with 'true' workers, might have allowed termites to outcompete other species and increased their propensity for reproductive and ecological specialization (definitions 2 and 3). This scenario is classically put forward at the ecological timescales, wherein the higher dominance and the large ecological success of termites are known, and could be here translated at the macroevolutionary timescale. Interestingly, all the diversification analyses (three different types of models) show that both speciation and extinction rates were higher for eusocial lineages with 'true' workers resulting in a higher net diversification rate for these lineages. Accordingly, our results suggest that societies with 'true' workers are not only more successful at ecological timescales but also over millions of years, which further implies that both organism-and species-level traits (aggregate and emergent traits, respectively) act on species selection.

Other factors are also probably at play in Dictyoptera diversification and they should be examined in future research [START_REF] Jablonski | Species selection: theory and data[END_REF][START_REF] Benton | The Red Queen and the Court Jester: species diversity and the role of biotic and abiotic factors through time[END_REF]. The diversification of angiosperms in the Cretaceous is often proposed to have played an important role in insect diversification, especially social insects [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF][START_REF] Cardinal | Bees diversified in the age of eudicots[END_REF].

Interestingly, molecular dating phylogenies of termites concur to show that termites appeared in the very Late Jurassic -Early Cretaceous [START_REF] Ware | The effects of fossil placement and calibration on divergence times and rates: an example from the termites (Insecta: Isoptera)[END_REF][START_REF] Bourguignon | The evolutionary history of termites as inferred from 66 mitochondrial genomes[END_REF]Legendre et al. 2015b). Although there is a debate on the origin of angiosperms, from the Early Cretaceous to Triassic (e.g. [START_REF] Beaulieu | Heterogeneous rates of molecular evolution and diversification could explain the Triassic age estimate for angiosperms[END_REF][START_REF] Magallón | A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity[END_REF]Silvestro et al. 2015b;Foster et al. 2017), the main period of diversification always occurred in the Cretaceous and coincides well with the rise of many extant insect clades. It is nonetheless questioned whether the angiosperm radiation has boosted insect diversification [START_REF] Rainford | Diet evolution and clade richness in Hexapoda: a phylogenetic study of higher taxa[END_REF][START_REF] Condamine | Global patterns of insect diversification: towards a reconciliation of fossil and molecular evidence?[END_REF]. The increase in termite diversification in the Cretaceous can be associated with the rise to dominance of angiosperms; a similar scenario is proposed for another important eusocial insect clade, the ants [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF]. Other dramatic landscape modifications that happened in the last 300 million years might have promoted shifts in diversification in termites. For instance, global climate changes or the break-up of ancient continents [START_REF] Bourguignon | The evolutionary history of termites as inferred from 66 mitochondrial genomes[END_REF] could have induced different dynamics of diversification between non-eusocial and eusocial lineages. Future studies could investigate the impact of such important environmental changes using both phylogenies and fossil data. Also, to better understand these intricate mechanisms, approaches integrating various factors such as caste and colony size but also foraging types, dispersal abilities, diets and digestive capacity [START_REF] Abe | Evolution of life types in termites[END_REF][START_REF] Jeanson | Positive feedback, convergent collective patterns and social transitions in arthropods[END_REF][START_REF] Brune | The gut microbiota of termites: Digesting the diversity in the light of ecology and evolution[END_REF], should be used to sort out the effect of these factors. Regarding diet specialization, which is strongly related to gut anatomy and microbiota composition [START_REF] Bignell | On the elevated intestinal pH of higher termites (Isoptera : Termitidae)[END_REF][START_REF] Köhler | High-resolution analysis of gut environment and bacterial microbiota reveals functional compartmentation of the gut in wood-feeding higher termites (Nasutitermes spp.)[END_REF][START_REF] Mikaelyan | Diet is the primary determinant of bacterial community structure in the guts of higher termites[END_REF], the concomitance of termite and angiosperm diversifications brings credit to its importance in termite diversification because changes in intestinal anatomy and microbiota complexity would have opened new niches to diversify in. [START_REF] Rabosky | Ecological limits and diversification rate: alternative paradigms to explain the variation in species richness among clades and regions[END_REF] underlined the possible change in species "carrying capacity" for some clades in a given area. Instead of playing on diversification rates per se, may eusociality impact ecological limits in a given environment, which would in turn increase termite diversification? This remains a possibility that could be linked to our side-result that net diversification rates within termites increase towards the present, which is unusual, but the underlying mechanisms are still to be deciphered. The apparent increased diversification rate in termites at the present can be caused by a phenomenon, called the pull of the recent, which is due to a more complete sampling of recent and still extant species than in the past (Jablonsky 2008). This is slightly different from the pull of the present that is due to constant birth-death models showing an upward turn in species number towards the present [START_REF] Nee | Birth-death models in macroevolution[END_REF]. Most of the phylogenies are better explained by a slowdown of speciation rates towards the present (e.g. [START_REF] Phillimore | Density-dependent cladogenesis in birds[END_REF][START_REF] Morlon | Inferring the dynamics of diversification: a coalescent approach[END_REF]). An increase in speciation through time is rarely inferred, thus the increase in termite speciation goes against the classical pattern. We do not think the increase in speciation can be attributed to the pull of the present/recent because the increase started tens of million years ago, while the effect of the pull of the present/recent is usually considered to be important in the last million years [START_REF] Nee | Birth-death models in macroevolution[END_REF]. The species sampling is relatively homogeneous across the phylogenetic tree, except for the Panesthiinae which are more sampled than the rest, but more importantly the termite subtree is not better sampled than the other groups (Legendre et al. 2015b). Therefore, we think the impact of the pull of the present/recent is limited here. Finally, even though there were probably multiple origins of the 'true' worker caste in termites [START_REF] Legendre | The phylogeny of termites (Dictyoptera: Isoptera) based on mitochondrial and nuclear markers: Implications for the evolution of the worker and pseudergate castes, and foraging behaviors[END_REF][START_REF] Legendre | Phylogenetic analyses of termite postembryonic sequences illuminate caste and developmental pathway evolution[END_REF], Dictyoptera offers a single replicate to investigate the role of eusociality in clade diversification. We expect advances in our understanding of the role of eusociality at the species selection level, when similar studies would be conducted with other eusocial organisms, especially in the insect order Hymenoptera [START_REF] Peters | Evolutionary history of the Hymenoptera[END_REF] where ants, bees and wasps constitute conspicuous elements of ecosystems and have evolved eusociality independently.

SUPPLEMENTARY MATERIAL

Supplementary material can be found in the Dryad Digital Repository at http://datadryad.org, http://dx.doi.org/10.5061/dryad.2tg34. The table reports the models applied, their number of parameters (NP), the means for the log-likelihood (logL) and the corrected Akaike Information Criterion (AICc), the difference of AICc between the best model (lowest AIC) and a given model (ΔAICc), the Akaike weight (AICω) and the values for each parameter of the corresponding model based on 100 dated trees. The best model is highlighted in bold. Standard errors for each parameter estimate are reported in Table S2. 0 = non-eusocial and 1 = eusocial; λ = speciation rate (one for each character state, λ1 is the speciation rate for eusocial species); µ = extinction rate (one for each character state, µ1 is the extinction rate for eusocial species); and q = transition rate between character states. Notes: 1 = gut with no specialized microbiota for lignocellulose digestion, 2* = gut with cellulolytic flagellates, and 3* = gut with entirely prokaryotic microbiota; λ = speciation rate (one for each character state, λ1 is the speciation rate for species without microbiota); µ = extinction rate (one for each character state, µ1 is the extinction rate for species without microbiota); and q = transition rate between character states. * Results of alternative coding strategies (see main text) are provided in Figs. S11-S12.

Figure 2 .

 2 Figure 2. Eusociality in termites: life cycle (a) and (b) evolutionary origin associated to

  model, where speciation is constant and extinction is null; (ii) a constant birth-death model, where speciation and extinction rates are constant; (iii) a variable speciation rate model without extinction; (iv) a variable speciation rate model with constant extinction; (v) a rateconstant speciation and variable extinction rate model; and (vi) a model in which both speciation and extinction rates vary. Models were compared by computing the ML score of each model and the resulting Akaike information criterion corrected by sample size (AICc).

Figure 3 .

 3 Figure 3. Diversification and eusociality in Dictyoptera. (a) Time-calibrated phylogeny of

TABLE 1 .

 1 Support for an increase of diversification rates for termites, especially for lineages including 'true' workers. The table reports the log-likelihood (logL) and the corrected Akaike Information Criterion (AICc) calculated for the whole phylogeny of Dictyoptera (assuming no shift of diversification) and seven evolutionary scenarios including a shift of diversification (as specified in the table). For each model, the parameter values of the best-fitting model are reported (rates are in events/Myr). The 'backbone' rates represent the estimated diversification parameters for the tree containing all Dictyoptera species except the termite subclades. The termite subclades are analyzed aside. The logL, the number of parameters (NP), and AICc are thus the sum of the logL/parameters/AICc of the backbone tree plus the logL/parameters/AICc of the termite subclade. The best-fitting scenario is determined by ∆AIC and AIC and includes a shift within the termites that only include species with 'true' workers in their societies (scenario 7, highlighted in bold).

					Scenario with a rate shift located at:
	Models	Whole phylogeny	the MRCA of Isoptera	the MRCA of Euisoptera	the MRCA of Kalotermitidae and Neoisoptera	the MRCA of Neoisoptera	the MRCA of Neoisoptera minus Rhinotermitinae	the MRCA of Termitidae, Coptotermitinae and Heterotermitinae	the MRCA of Termitidae
	NP	3	4	4	5	4	5	5	4
	logL	-3889.247	-3709.158	-3706.857	-3693.107	-3658.84	-3659.554	-3651.496	-3691.284
	AICc	7784.526	7426.385	7421.781	7396.308	7325.754	7329.206	7313.092	7390.653
	∆AIC	471.434	113.293	108.689	83.216	12.662	16.114	0	77.561
	AICω	0	0	0	0	≈0.002	<0.0001	0.998	0
	Subclade speciation rate	0.191	0.177	0.178	0.185	0.186	0.191	0.193	0.197
	Subclade extinction rate	0.175	0	0	0	0	0	0	0
	Subclade net div. rate	0.016	0.177	0.178	0.185	0.186	0.191	0.193	0.197
	Backbone speciation rate	-	0.044	0.044	0.057	0.044	0.058	0.058	0.084
	Backbone extinction rate	-	0	0	0.024	0	0.024	0.023	0.06
	Backbone net. div. rate	-	0.044	0.044	0.033	0.044	0.034	0.035	0.024
	Notes:

TABLE 3 .

 3 Supports for eusociality with true workers as drivers of diversificationNotes: 1 = non-eusocial or societies without worker, 2 = societies with pseudergates, and 3 = societies with true workers; λ = speciation rate (one for each character state, λ1 is the speciation rate for species without worker); µ = extinction rate (one for each character state, µ1 is the extinction rate for species without worker); and q = transition rate between character states.

	Models	NP	logL	AICc	∆AIC AICω	λ1	λ2	λ3	µ1	µ2	µ3	q12	q13	q21	q23	q31	q32
	Null model	3 -3889.534 7785.100 373.94	0		0.1719			0.1525				4.21E-05
	all λ are varying	5 -3717.017 7444.114 32.95	0	0.0683	0.0727	0.1699		0.0438				2.43E-05
	all µ are varying	5 -3751.789 7513.657 102.49	0		0.1214		0.1020 0.0974	0.0000			1.43E-05
	all q are varying	8 -3879.999 7776.191 365.03	0		0.1713			0.1519		1.29E-05	6.67E-07	1.67E-06	0.00161	4.64E-06	0.00023
	all λ and µ are varying 7 -3706.089 7426.327 15.16 ≈0.005 0.0637	0.0580	0.3198 0.0384 0.0266	0.2267			3.85E-05
	all λ and q are varying 10 -3703.355 7427.005 15.84 ≈0.005 0.0645	0.0693	0.1736		0.0394		1.55E-05	1.33E-07	7.90E-07	0.00155	9.81E-06	2.73E-06
	all µ and q are varying 10 -3738.133 7496.560 85.40	0		0.1258		0.1067 0.1030	0.0000	1.87E-05	7.14E-08	2.63E-07	0.00043	2.73E-06	0.00044
	All rates are free	12 -3693.372 7411.163	0	0.99	0.0649	0.0378	0.3196 0.0398 0.0000	0.2244	1.09E-05	2.43E-07	6.45E-07	0.00224	1.32E-06	0.00027

TABLE 4 .

 4 Supports for gut microbiota as drivers of diversification

	Models	NP	logL	AICc	∆AIC AICω	λ1	λ2	λ3	µ1	µ2	µ3	q12	q13	q21	q23	q31	q32
	Null model	3	-3840.575 7687.181 274.01	0		0.1959			0.1764				3.03E-05
	all λ are varying	5	-3710.247 7430.573 17.41	0	0.1166 0.1270 0.2157		0.0941				1.72E-05
	all µ are varying	5	-3731.518 7473.115 59.95	0		0.1642		0.1456	0.1373 0.0447			1.21E-05
	all q are varying	8	-3835.418 7687.028 273.86	0		0.1957			0.1763		2.40E-05 6.75E-07 4.36E-06 0.00046 1.88E-06	3.10E-06
	all λ and µ are varying	7	-3699.510 7413.168	0	0.703 0.1070 0.1894 0.3022 0.0832	0.1652 0.1979			2.39E-05
	all λ and q are varying 10	-3703.777 7427.847 14.68	0	0.1147 0.1254 0.2157		0.0921		1.98E-05 2.74E-07 3.45E-06 0.00025 3.21E-06	2.81E-06
	all µ and q are varying 10	-3735.700 7491.694 78.53	0		0.1683		0.1500	0.1304 0.0555 2.01E-05 7.33E-07 3.67E-06 0.01504 2.54E-06	0.00220
	All rates are free	12	-3695.237 7414.892	1.72	0.297 0.1037 0.1730 0.3033 0.0795	0.1472 0.1975 3.11E-05 4.02E-07 4.87E-06 0.00029 4.52E-06	2.54E-06

ACKNOWLEDGMENTS

We thank the associate editor Brian Wiegmann and two anonymous referees for the excellent comments that improved the study. We also thank Giovanni Fagua, Gaël Kergoat, Hélène Morlon, André Nel and Felix Sperling for proofreading and commenting on this manuscript.

FUNDING

Financial support was provided by a Marie Curie International Outgoing Fellowship (project 627684 BIOMME) to FLC.

Cockroaches