<figure><img src="image_1.png" /></figure>

HAL
open science

# A Jurassic Gondwanan origin with a Mesozoic vicariance and Cenozoic dispersals explains the biogeographic history of crown webspinners (Insecta: Embioptera)

Corentin Jouault, Fabien L Condamine, Frédéric Legendre

<figure><img src="image_2.png" /></figure>

## To cite this version:

Corentin Jouault, Fabien L Condamine, Frédéric Legendre. A Jurassic Gondwanan origin with a Mesozoic vicariance and Cenozoic dispersals explains the biogeographic history of crown webspinners (Insecta: Embioptera). Evolutionary Journal of the Linnean Society, 2025, 4, 10.1093/evolinskiean/kzaf011. mnhn-05283283

HAL Id: mnhn-05283283

https://mnhn.hal.science/mnhn-05283283v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Original Article

# A Jurassic Gondwanan origin with a Mesozoic vicariance and Cenozoic dispersals explains the biogeographic history of crown webspinners (Insecta: Embioptera)

Corentin Jouault<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,*, Fabien L. Condamine<fnref n="3" />,+, Frédéric Legendre<fnref n="2" />,†

## ABSTRACT

Embioptera, commonly called webspinners, represent a small and understudied order of insects whose evolutionary history is poorly understood. To date, neither the timing of their origin, diversification, nor biogeographic history has been thoroughly investigated. In this study, we establish the ground-work for addressing these gaps by estimating divergence times for most extant webspinner lineages and reconstructing their biogeographic history using molecular and morphological datasets. Regardless of the sampling fractions (representing either the current or estimated diversity of webspinners), the clock models, and the inclusion or exclusion of morphological data, our divergence time estimates indicated that webspinners originated at around 280 Mya in the Permian, whereas the crown Embioptera appeared later, around 170 Mya in the Middle Jurassic. Our sensitivity analyses indicated that the divergence time estimates were robust to known biases in molecular dating analyses. With this new time-calibrated phylogenetic framework, we investigated the ancestral geographic origin of the order and reconstructed the historical biogeography of webspinners. Our analyses confirmed the expected Gondwanan origin of the crown group evident in the distribution of extant families, and a first radiation in Australia (Australembiidae). Our results suggested that the distribution of extant Embioptera results from a combination of vicariance, due to the break-up of Gondwana, and dispersal events following continental collisions during the Cenozoic. Through our analyses, we present a first comprehensive scenario for the diversification and origin of Embioptera, which provides testable hypotheses for further studies.

**Keywords:** Bayesian inference; biogeography; divergence times; node calibration; sampling fraction

## INTRODUCTION

Within the insect superorder Polyneoptera, the Embioptera are called webspinners because of their ability to spin silk, which they use to create domiciles (Grimaldi and Engel 2005, Misof et al. 2014, Wipfler et al. 2019; Fig. 1). Today, the order Embioptera is depauperate in terms of species number, as only about 400 valid extant species are known (Szumik et al. 2008, Maehr 2000). However, their extant diversity is underestimated, and taxonomists assume that the order diversity may be at least twice and possibly up to four times the current number of known species (Ross 1991). Female embiopterans are wingless and often found on tree or rock surfaces, under rocks, or in leaf litter, sometimes with their

offspring. They show varying degrees of sociality, sometimes living alone and sometimes in groups of females (Edgerly 1997, Morinaga et al. 2023). Males are often winged but rarely collected, which not only complicates our understanding of the embiopteran subclades delineation, but also prevents a clear estimate of the diversity within the order (Ross 1991, Szumik et al. 2008, Miller et al. 2012). Most of the characters used to distinguish the species or the higher clades are found in the male head and terminalia (Szumik et al. 2008, Miller et al. 2012).

Molecular dating studies suggest that Embioptera diverged from a common ancestor shared with the Phasmatodea (stick insects) during the Jurassic (Misof et al. 2014), Triassic

<footnote n="1">¹Oxford University Museum of Natural History, University of Oxford, Parks Road, Oxford OX1 3PW, United Kingdom</footnote>

<footnote n="2">²Institut de Systématique, Évolution, Biodiversité (ISYEB), MNHN, CNRS, SU, EPHE-PSL, UA, CP50, Paris F-75005, France</footnote>

<footnote n="3">³Institut des Sciences de l'Évolution de Montpellier (UMR 5554), Université de Montpellier, CNRS, Montpellier 34095, France</footnote>

<footnote n="4">⁴Géosciences Rennes (UMR 6118), Université de Rennes, CNRS, Rennes F-35000, France</footnote>

<footnote>*Corresponding author. Oxford University Museum of Natural History, University of Oxford, Parks Road, Oxford OX1 3PW, United Kingdom.</footnote>

<footnote>E-mail: corentin.jouault@oum.ox.ac.uk</footnote>

<footnote>†These authors jointly supervised this work.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Diversity of extant and extinct webspinners. A, extant Oligotomidae (© CSIRO). B, fossil nymph in Mid-Cretaceous Kachin amber (courtesy of Michele Baldi). C, Sorellembiidae in Mid-Cretaceous Kachin amber (courtesy of Enrico Bonino). D, adult on tree bark (© CSIRO).</figcaption></figure>

(Tong *et al.* 2015, Tihelka *et al.* 2020), or Permian (Montagna *et al.* 2019). Interestingly, their crown group is estimated to have arisen much later during the Cretaceous (Simon *et al.* 2019, Tihelka *et al.* 2020), suggesting that an important extinct diversity of stem lineages may have existed between the origin of the order and their crown group. However, no study has focused on estimating the time of divergence of the order and its constituent lineages, with most previous estimates derived from time-calibrated phylogenies reconstructed to investigate their relationships with Phasmatodea and the origin of the latter order (e.g. Simon *et al.* 2019, Tihelka *et al.* 2020). Similarly, the impact of uncertainty in the ‘true’ number of Embioptera species on divergence time estimates has not been explored, although it is directly involved in dating analyses under the parameter defining the sampling probability of terminal lineages (e.g. Zhang *et al.* 2016, Luo *et al.* 2023). Recent simulations have suggested that increasing taxon-sampling density does not necessarily improve divergence time estimates in tip-dating analyses (Luo *et al.* 2023). On the other hand, increasing taxon sampling in node dating analyses has often improved topological inference (e.g. Heath *et al.* 2008, Soares and Schrago 2012) and the accuracy of branch-length estimates by reducing the node-density effect (e.g. Hugall and Lee 2007). However, empirical testing of these results in node-dating analyses remains limited for insects (Jouault *et al.* 2025).

The fossil record of the Embioptera, including over 30 species and recent discoveries from the Permian, offers valuable insights into the timing of the origin of the order and its extant lineages (The Paleobiology Database, https://paleobiodb.org). These fossils challenge previous assumptions about the age of the order and its crown group, prompting a closer look at the timing of the origin of Embioptera and its constituent lineages (e.g. Engel *et al.* 2011, 2016, Cui *et al.* 2020). They also underscore the necessity to re-evaluate the classification of numerous taxa known from the Permian or Triassic, considering their wing venation and

preserved body characters. For all these reasons, a comprehensive study of the phylogenetic relationships and dating estimates in Embioptera appears crucial. A robust topology would provide a backbone to study, among other things, the biogeographic history of webspinners.

The current geographic distribution of extant Embioptera clades reflects a complex evolutionary history, with some families being widespread and others restricted in distribution (Szumik *et al.* 2008). Currently, 18 families of Embioptera are recognized as valid (Hopkins 1900), though additional ‘hidden’ lineages and undocumented diversity are likely present within the order (Miller *et al.* 2012). Anisembiidae are confined to the New World (North + South America) (Szumik *et al.* 2008, fig. 10), and Notoligotomidae and Australembiidae are exclusively present in Australia, whereas Embiidae or Oligotomidae have a nearly global distribution (Szumik *et al.* 2008, figs 7, 13, 17). Previous phylogenetic reconstructions suggested Australembiidae as the earliest diverged embiopteran family, with the subsequent clade primarily comprising Neotropical lineages (Anisembiidae, Notoligotomidae, Archembiidae, and Andesembiidae) (Miller *et al.* 2012). Whether these distributions and relationships support a Gondwanan origin for the order remains to be investigated. To date, no study has thoroughly explored the early historical biogeography of the clade, leaving gaps in our understanding of the ancestral range of the extant webspinner clades and the order itself.

To provide new insights into the evolutionary history of webspinners, we propose the first time-calibrated phylogeny of the order, allowing the timing of the origin of stem and crown Embioptera to be investigated. To account for uncertainty in the ‘true’ number of Embioptera species, the impact of model choice, and the inclusion of a morphological partition on divergence time estimates, we test multiple combinations in a Bayesian framework. Lastly, we use the newly reconstructed time-calibrated phylogeny and the known distribution of current lineages to
explore the biogeographic evolution of webspinners and test for
a Gondwanan origin.

# MATERIAL AND METHODS

## Molecular and morphological data and matrix assembly

The molecular dataset was modified from Miller et al. (2012) and is composed of six outgroup taxa distributed in the Dictyoptera, Grylloblattodea, Mantophasmatodea, Phasmatodea, and 81 Embioptera species composing the ingroup. DNA sequences were retrieved from five markers: 16S rRNA (580 bp), 18S rRNA (2351 bp), 28S rRNA (2808 bp), cytochrome oxidase subunit I (coxI; 1282 bp), and histone H3 (328 bp). After cross-referencing our dataset with NCBI (https://www.ncbi.nlm.nih.gov), we did not find substantial new molecular data that would warrant inclusion. For instance, adding a species represented by only a single gene (e.g. coxI) would likely weaken our dataset and potentially introduce branching artifacts, complicating the dating analysis. The dataset of Miller et al. (2012) remains the most comprehensive and provides the best balance between quality and taxonomic coverage. We excluded the taxon EB142_Embiidae_sp from our analyses due to its undetermined nature and its ambiguous and unresolved position in prior analyses (Miller et al. 2012). It has been suggested that this specimen may belong to an unrecognized family or new clade (Miller et al. 2012). Because we could not examine the specimen directly to check it further, we chose to omit it from the analysis. The sequences were aligned using Muscle 3.8 (Edgar 2004). The concatenated matrix was generated in Aliview 1.28 (Larsson 2014). The morphological data were also retrieved from Miller et al. (2012) and encompass 95 unordered, unweighted characters (sometimes multistates). Two matrices were established using Mesquite 3.61 (Maddison and Maddison 2019): (1) with molecular data only; and (2) with a combination of molecular and morphological data. In both, gaps and characters coded as inapplicable were treated as missing data and coded ‘-’, while unknown characters were coded ‘?’. The final concatenated molecular matrix comprises 7349 characters and 87 taxa. The final concatenated matrix encompassing molecular and morphological data comprises 7444 characters and 87 taxa (Figshare data repository (Jouault 2025a)).

## Maximum likelihood phylogenetic analyses

We performed maximum likelihood analyses using IQ-TREE 2.2.2.6 (Minh et al. 2020) with a dataset partitioned following the results of the PartitionFinder analysis. For each partition, best-fitting substitution models (-m MFP) were selected by ModelFinder (Kalyaanamoorthy et al. 2017). ModelFinder computes the log-likelihoods of an initial parsimony tree for many different models, and the Akaike information criterion (AIC), corrected Akaike information criterion (AICc), and the Bayesian information criterion (BIC). Then, ModelFinder chooses the model that minimizes the BIC score (Subset 1: GTR+F+I+G4; Subset 2: GTR+F+I+R4; Subset 3: TIM2+F+R3; Subset 4: SYM+I+R4). Phylogenetic analyses were carried out with the thorough and slower nearest neighbour interchange (NNI) search (-allnni). The perturbation strength was defined as 0.2 as recommended for datasets containing short sequences (-pers 0.2), and the number of unsuccessful iterations to stop was set to 500 (-nstop 500).

The final analysis was repeated 10 times to test the robustness and stability of the results and to avoid the tree search getting stuck in a local optimum. Support for nodes was evaluated with 1000 ultra-fast bootstrap (UFBS) approximations, considering UFBS values ≥ 95 as strong node support (Hoang et al. 2018). The strong support and monophyly of families are used to justify constraints implemented for Bayesian inferences to optimize computational time. Files are available in the Figshare data repository (Jouault 2025a).

## Bayesian phylogenetic and dating inferences

To estimate divergence times and infer their credibility intervals, we performed Bayesian relaxed-clock analyses using MrBayes 3.2.7a (Huelsenbeck and Ronquist 2001, Ronquist et al. 2012b, Ronquist and Huelsenbeck 2003). As MrBayes does not include as many evolutionary models as IQ-TREE, we used Partition-Finder 2 (Lanfear et al. 2017) to identify partitions that could be combined (settings: branch lengths = unlinked, models = mrbayes, model_selection = BIC, search = greedy), and the best-fitting evolutionary model for all partitions (GTR+I+G). Four partitions were suggested by PartitionFinder 2 as follows: (16S), (18S, 28S), (third position of coxI) and (H3-1, H3-2, H3-3, coxI-1, coxI-2). Morphological characters were grouped in one partition. The model parameters for the molecular partitions were set as nst = 6, rates = invgamma, and covarion = no (i.e. GTR+I+G), with all substitution model parameters unlinked across partitions. Morphological characters were analysed with a one-parameter Markov model (Lewis 2001) with a gamma rate variation across characters (coding = variable; rates = gamma). All analyses were computed using two different clock models: (1) an independent gamma relaxed clock model (Lepage et al. 2007) [prset clockratepr = log-norm(-7,0.6), prset clockvarpr = igr, prset igrvarpr = exp(10)] (Table 1A); and (2) an autocorrelated relaxed lognormal clock model (Thorne and Kishino 2002) [prset clockratepr = log-norm(-7,0.6), prset clockvarpr = tk02, prset tk02varpr = exp(10)] (Table 1B). Since substitution model selection impacts divergence time estimates and topology inference (Jouault et al. 2023), we performed an additional series (i.e. not using the GTR+I+G) of analyses without substitution models selected a priori, for the molecular partitions, but using the reversible jump Markov chain Monte Carlo (rjMCMC) model. The rjMCMC model, a form of model averaging across different dimensions of parameter space, explores all possible time-reversible substitution models, thereby incorporating uncertainty in model selection. Through rjMCMC, MrBayes integrates all substitution models (six rates and 203 possible time-reversible models), with the Markov chain sampling a nucleotide substitution model in proportion to its marginal likelihood, spending more time in the best likelihood (Huelsenbeck et al. 2004).

Similarly, we tested for the impact of the sampling fraction (SF) parameter, because the number of known and putative species frequently differ in insect lineages, including Embioptera (Aberlenc et al. 2020). We used two different sampling fractions: the first one was set to 0.21 and accounted for the current diversity of webspinners (about 400 species; Miller et al. 2012), and the second was set to 0.10375 and accounted for the estimated, yet undiscovered, diversity of the order (about 800 species; Ross 1991). The prior probability distribution on branch lengths was set to
<table><caption>**Table 1.** Time divergence estimates of the main clades in Embioptera, according to different data sets, substitution models, sampling fractions and in-group calibrations (analyses in bold are illustrated in Fig. 3).</caption><thead><tr><th colspan="13">IGR clock model (Part.A)</th></tr><tr><th>Data type</th><th>Models</th><th>Sampling fraction</th><th>Calibrations</th><th>Number of calibrations</th><th>Crown-Embioptera</th><th>Andesemblidae + archenblidae</th><th>Anisembidae</th><th>Archembidae</th><th>Australembidae</th><th>Clothodidae</th><th>Embidae</th><th>Notolectomidae</th><th>Oligotomidae</th><th>Pitocerembidae</th><th>Scelambidae</th><th>Teratembidae</th></tr></thead><tbody><tr><td rowspan="10">Molecular</td><td>GTR+I+G</td><td>≈ 0.1</td><td>Stem</td><td>3</td><td>156.53<br>[131.5-186.76]</td><td>101.88<br>[72.1-132.75]</td><td>96.6<br>[70.93-123.97]</td><td>70.45<br>[41.05-101.54]</td><td>103.88<br>[70.52-138.4]</td><td>44.97<br>[17.7-73.66]</td><td>91.02<br>[13.23-53.21]</td><td>30.44<br>[68.12-115.8]</td><td>86.39<br>[66.26-107.4]</td><td>43.83<br>[25.33-65.56]</td><td>95.04<br>[75.54-117.45]</td><td>59.29<br>[35.99-84.34]</td></tr><tr><td>GTR+I+G</td><td>≈ 0.21</td><td>Stem</td><td>3</td><td>156.25<br>[131.73-185.18]</td><td>102.6<br>[71.53-132.87]</td><td>97.46<br>[72.45-125.3]</td><td>71.8<br>[42.86-102.35]</td><td>102.74<br>[70.33-138.27]</td><td>46.52<br>[18.81-75.05]</td><td>90.84<br>[68-115.24]</td><td>30.76<br>[12.99-53.3]</td><td>86.08<br>[67.19-107.46]</td><td>43.66<br>[25.24-65.91]</td><td>94.9<br>[75.29-117.31]</td><td>59.65<br>[36.69-84.45]</td></tr><tr><td>RJMCMC</td><td>≈ 0.1</td><td>Stem</td><td>3</td><td>156.09<br>[129.78-183.83]</td><td>100.96<br>[71.57-133.6]</td><td>96.75<br>[71.81-124.36]</td><td>70.47<br>[43.23-99.45]</td><td>102.55<br>[71.76-140.75]</td><td>46.09<br>[17.77-74.83]</td><td>91<br>[67.74-115.75]</td><td>29.63<br>[13.26-52.48]</td><td>85.97<br>[65.89-108.8]</td><td>43.61<br>[25.39-65.45]</td><td>95.2<br>[76.48-118.28]</td><td>59.54<br>[36.6-84.61]</td></tr><tr><td>RJMCMC</td><td>≈ 0.21</td><td>Stem</td><td>3</td><td>155.95<br>[130.18-185.42]</td><td>102.38<br>[71.02-134.2]</td><td>97.54<br>[72.09-124.2]</td><td>71.53<br>[42.35-101.77]</td><td>102.9<br>[70.48-140.08]</td><td>45.94<br>[18.68-75.43]</td><td>91.29<br>[68.88-114.87]</td><td>30.71<br>[12.9-55.1]</td><td>86.14<br>[66.62-107.26]</td><td>43.76<br>[25.28-65.39]</td><td>95.01<br>[75.32-118.22]</td><td>58.72<br>[36.68-83.46]</td></tr><tr><td>Molecular</td><td>GTR+I+G</td><td>≈ 0.1</td><td>Crown</td><td>5</td><td>170.51<br>[142.08-200.23]</td><td>109<br>[76.59-142.8]</td><td>104.04<br>[74.97-133.34]</td><td>75.46<br>[43.67-108.78]</td><td>112.0<br>[74.9-152.61]</td><td>49.21<br>[17.95-80.02]</td><td>100.37<br>[75.7-126.77]</td><td>32.61<br>[14.64-58.65]</td><td>103.62<br>[98-117.96]</td><td>47.77<br>[27.57-71.68]</td><td>107.31<br>[98-124.88]</td><td>64.29<br>[39.06-91.67]</td></tr><tr><td rowspan="10">Molecular + Morphological</td><td>GTR+I+G</td><td>≈ 0.21</td><td>Crown</td><td>5</td><td>169.9<br>[144.05-200.05]</td><td>110.88<br>[76.52-143.86]</td><td>106.8<br>[77.45-135.23]</td><td>77.33<br>[43.29-110.88]</td><td>113.17<br>[73.79-152.05]</td><td>50.67<br>[20.63-83.2]</td><td>99.4<br>[72.86-125.8]</td><td>33.79<br>[14.47-59.84]</td><td>103.43<br>[98-117.7]</td><td>47.33<br>[27.45-70.53]</td><td>106.42<br>[98-123.66]</td><td>65.29<br>[41.04-92.45]</td></tr><tr><td>RJMCMC</td><td>≈ 0.1</td><td>Crown</td><td>5</td><td>169.93<br>[142.8-199.47]</td><td>108.6<br>[75.19-142.97]</td><td>103.8<br>[76.73-136.6]</td><td>75.09<br>[42.91-109.47]</td><td>113.13<br>[74.74-151.91]</td><td>48.31<br>[18.73-79.4]</td><td>98.43<br>[72.64-123.13]</td><td>32.47<br>[14.77-58.39]</td><td>103.42<br>[98-118.47]</td><td>46.73<br>[27.11-69.77]</td><td>106.5<br>[98-123.08]</td><td>65.33<br>[40.02-93.58]</td></tr><tr><td>RJMCMC</td><td>≈ 0.21</td><td>Crown</td><td>5</td><td>170.64<br>[143.93-200.21]</td><td>111.72<br>[80.22-143.89]</td><td>106.78<br>[78.33-135.14]</td><td>77.34<br>[47.05-110.91]</td><td>113.45<br>[67.97-150.82]</td><td>50.97<br>[21.14-84.24]</td><td>98.43<br>[72.92-124.79]</td><td>33.98<br>[14.38-59.84]</td><td>103.39<br>[98-118.78]</td><td>48.26<br>[28.79-70.3]</td><td>106.66<br>[40.27-91.94]</td><td>64.96<br>[40.67-91.94]</td></tr><tr><td>GTR+I+G/Mkv+G</td><td>≈ 0.1</td><td>Stem</td><td>3</td><td>159.95<br>[135.39-189.16]</td><td>NA</td><td>100.07<br>[75.86-129.03]</td><td>74.61<br>[46.13-103.56]</td><td>100.2<br>[68.52-135.46]</td><td>45.67<br>[18.41-73.81]</td><td>91.47<br>[67.08-115.36]</td><td>27.69<br>[12.32-48.55]</td><td>84.27<br>[64.46-105.09]</td><td>41.98<br>[25.18-62.67]</td><td>95.04<br>[74.72-116.94]</td><td>58.78<br>[37.03-84.87]</td></tr><tr><td>GTR+I+G/Mkv+G</td><td>≈ 0.21</td><td>Stem</td><td>3</td><td>158.69<br>[133.49-189.37]</td><td>NA</td><td>100.53<br>[75.05-126.43]</td><td>75.39<br>[45.88-103.81]</td><td>99.18<br>[67.78-133.57]</td><td>46.97<br>[21.12-75.53]</td><td>90.67<br>[67.19-115.66]</td><td>27.63<br>[12.32-49.38]</td><td>84.26<br>[65.3-104.88]</td><td>41.39<br>[23.9-61.83]</td><td>94.46<br>[74.33-116.97]</td><td>58.84<br>[35.65-82.3]</td></tr><tr><td>RJMCMC/Mkv+G</td><td>≈ 0.1</td><td>Stem</td><td>3</td><td>159.49<br>[134.78-190.1]</td><td>NA</td><td>100.64<br>[75.7-130.67]</td><td>74.99<br>[46.15-103.3]</td><td>99.74<br>[66.59-133.13]</td><td>46.35<br>[20.85-74.6]</td><td>90.84<br>[67.08-116.06]</td><td>27.64<br>[12.68-49.26]</td><td>84.64<br>[64.92-105.4]</td><td>42.46<br>[25.24-63.18]</td><td>95.16<br>[74.11-117.03]</td><td>59.27<br>[36.85-83.75]</td></tr><tr><td>RJMCMC/Mkv+G</td><td>≈ 0.21</td><td>Stem</td><td>3</td><td>160.2<br>[135.53-190]</td><td>NA</td><td>100.82<br>[75.19-126.92]</td><td>75.16<br>[46.47-104.74]</td><td>100.15<br>[69.65-135.27]</td><td>46.99<br>[20.92-76.4]</td><td>91.57<br>[68.65-114.99]</td><td>28.01<br>[12.41-48.97]</td><td>85.51<br>[66.09-106.17]</td><td>41.99<br>[24.62-62.88]</td><td>95.62<br>[74.6-116.54]</td><td>59.45<br>[37.08-83.1]</td></tr><tr><td rowspan="10">Molecular + Morphological</td><td>GTR+I+G/Mkv+G</td><td>≈ 0.1</td><td>Crown</td><td>5</td><td>174.55<br>[147.75-206.41]</td><td>NA</td><td>108.14<br>[80.28-138.19]</td><td>79.81<br>[47.87-113.3]</td><td>109.24<br>[73.19-148.66]</td><td>50.17<br>[20.13-82.02]</td><td>100.7<br>[74.88-126.02]</td><td>30.34<br>[14.26-53.79]</td><td>103.14<br>[98-117.29]</td><td>46.69<br>[27.29-69.78]</td><td>100.7<br>[98-125.37]</td><td>65.46<br>[39.76-91.92]</td></tr><tr><td>GTR+I+G/Mkv+G</td><td>≈ 0.1</td><td>Crown</td><td>5</td><td>174.82<br>[148.02-204.76]</td><td>NA</td><td>108.27<br>[80.73-137.08]</td><td>80.27<br>[50.6-114.32]</td><td>109.64<br>[73.75-146.35]</td><td>49.84<br>[19.29-82.5]</td><td>99.61<br>[74.38-125.88]</td><td>30.19<br>[13.77-53.08]</td><td>102.96<br>[98-117.12]</td><td>46.13<br>[27.05-69.33]</td><td>107.28<br>[98-125.01]</td><td>66<br>[40.91-92.21]</td></tr><tr><td>RJMCMC/Mkv+G</td><td>≈ 0.1</td><td>Crown</td><td>5</td><td>175.46<br>[147.59-207.29]</td><td>NA</td><td>109.05<br>[79.69-140.02]</td><td>80.03<br>[49.64-112.09]</td><td>109.95<br>[75.1-150.15]</td><td>51.39<br>[20.68-84.55]</td><td>99.88<br>[73.66-126.4]</td><td>31.18<br>[14.7-54.76]</td><td>103.47<br>[98-118.14]</td><td>46.3<br>[26.1-68.68]</td><td>107.66<br>[98-126.17]</td><td>66.02<br>[40.38-93.5]</td></tr><tr><td>RJMCMC/Mkv+G</td><td>≈ 0.21</td><td>Crown</td><td>5</td><td>174.11<br>[148.13-205.97]</td><td>NA</td><td>109.33<br>[81.39-139.5]</td><td>80.08<br>[47.63-112.7]</td><td>109<br>[75.56-147.47]</td><td>50.68<br>[21.11-81.49]</td><td>100.45<br>[74.81-126.46]</td><td>30.59<br>[13.85-54.47]</td><td>103.17<br>[98-117.76]</td><td>46.02<br>[26.13-69.46]</td><td>107.24<br>[98-125.4]</td><td>65.23<br>[40.2-92.24]</td></tr><tr><td>GTR+I+G</td><td>≈ 0.1</td><td>Variable</td><td>7</td><td>170.83<br>[142.99-199.73]</td><td>111.74<br>[77.82-144.44]</td><td>103.82<br>[77.19-131.27]</td><td>76.38<br>[44.69-109.9]</td><td>112.35<br>[76.42-151.8]</td><td>48.46<br>[19.52-80]</td><td>99.49<br>[75.26-128.28]</td><td>32.68<br>[14.97-59.09]</td><td>103.35<br>[98-117.63]</td><td>47.91<br>[26.8-71.19]</td><td>107.3<br>[98-125.41]</td><td>65.83<br>[40.08-95.11]</td></tr><tr><td>GTR+I+G</td><td>≈ 0.21</td><td>Variable</td><td>7</td><td>170.97<br>[145.16-202.98]</td><td>113.13<br>[80.82-145.27]</td><td>103.47<br>[77-132.07]</td><td>78.02<br>[49.78-110.88]</td><td>113.62<br>[79.33-151.97]</td><td>51.29<br>[19.89-83.89]</td><td>101.46<br>[73.53-127.74]</td><td>34.09<br>[14.57-59.64]</td><td>103.43<br>[98-117.16]</td><td>48.05<br>[27.66-71.71]</td><td>107.19<br>[98-125.13]</td><td>65.32<br>[40.81-92.05]</td></tr><tr><td>RJMCMC</td><td>≈ 0.1</td><td>Variable</td><td>7</td><td>171.79<br>[144.95-201.3]</td><td>112.25<br>[77.98-145.69]</td><td>103.19<br>[75.31-131.64]</td><td>77.04<br>[42.73-110.76]</td><td>113.81<br>[73.58-152.2]</td><td>49.41<br>[19.69-83.97]</td><td>99.99<br>[70.49-126.62]</td><td>33.77<br>[15.52-60.31]</td><td>103.5<br>[98-118.58]</td><td>48.39<br>[27.54-72.36]</td><td>107.45<br>[98-125.73]</td><td>66.27<br>[41.47-94.3]</td></tr><tr><td>RJMCMC</td><td>≈ 0.21</td><td>Variable</td><td>7</td><td>170.57<br>[142.65-200.75]</td><td>112.72<br>[80.22-146.66]</td><td>103.9<br>[78.26-130.02]</td><td>77.78<br>[45.32-111.71]</td><td>112.27<br>[75.21-152.04]</td><td>100.03<br>[73.09-125.7]</td><td>33.94<br>[14.03-60.41]</td><td>103.66<br>[98-118.34]</td><td>48.17<br>[28.31-70.95]</td><td>107.13<br>[98-125.37]</td><td>66.72<br>[40.36-94.46]</td></tr></tbody></table>

(Continued)
<table><caption>Table 1. (Continued)</caption><thead><tr><th colspan="13">IGR clock model (Part A)</th></tr><tr><th>Data type</th><th>Models</th><th>Sampling fraction</th><th>Calibrations</th><th>Number of calibrations</th><th>Crown-Embioptera</th><th>Andensembiida + Embioptera</th><th>Anensembiida</th><th>Archemiidae</th><th>Australembiidae</th><th>Embiidae</th><th>Clothodidae</th><th>Embiidae</th><th>Notoligotomidae</th><th>Oligotomidae</th><th>Ptilocerembiidae</th><th>Scelembiidae</th><th>Teratembiidae</th></tr></thead><tbody><tr><td rowspan="6">Molecular + Morphological</td><td>GTR+I+G/Mkv+G ≈ 0.1</td><td>Variable</td><td>7</td><td>177.54<br/>[152.06-208.69]</td><td>NA</td><td>84.61<br/>[84.7-141.44]</td><td>11.69<br/>[57.01-115.52]</td><td>11.69<br/>[75.24-150.14]</td><td>49.41<br/>[19.22-81.46]</td><td>101.17<br/>[75.57-127.98]</td><td>30.87<br/>[14.96-53.9]</td><td>103.62<br/>[98-118.45]</td><td>46.78<br/>[27.26-70.51]</td><td>108.14<br/>[98-126.8]</td><td>66.5<br/>[40.88-94.32]</td><td></td><td></td></tr><tr><td>GTR+I+G/Mkv+G ≈ 0.21</td><td>Variable</td><td>7</td><td>176.12<br/>[149.77-207.05]</td><td>NA</td><td>112.31<br/>[85.17-142.02]</td><td>85.36<br/>[56.06-114.97]</td><td>110.48<br/>[76.16-150.06]</td><td>51.16<br/>[22.29-83.28]</td><td>100.14<br/>[73.92-125.13]</td><td>31.9<br/>[14.52-55.19]</td><td>103.44<br/>[98-117.69]</td><td>47.3<br/>[27.58-72.28]</td><td>107.61<br/>[98-127.01]</td><td>66<br/>[39.64-93.42]</td><td></td><td></td></tr><tr><td>RJMCMC/Mkv+G ≈ 0.1</td><td>Variable</td><td>7</td><td>177.25<br/>[151.18-207.16]</td><td>NA</td><td>110.19<br/>[80.68-139.99]</td><td>83.92<br/>[54.49-114.56]</td><td>109.93<br/>[74.76-148.26]</td><td>50.87<br/>[19.34-82.01]</td><td>100.44<br/>[73.96-127.17]</td><td>30.36<br/>[15-54.82]</td><td>103.34<br/>[98-117.43]</td><td>46.32<br/>[25.91-68.41]</td><td>65.95<br/>[98-126.6]</td><td>65.95<br/>[40.69-92.87]</td><td></td><td></td></tr><tr><td>RJMCMC/Mkv+G ≈ 0.21</td><td>Variable</td><td>7</td><td>176.47<br/>[151.47-207.53]</td><td>NA</td><td>112.16<br/>[81.92-141.14]</td><td>83.73<br/>[53.56-114.25]</td><td>110.6<br/>[76.01-148.29]</td><td>51.64<br/>[22.08-83.39]</td><td>100.66<br/>[74.08-126.14]</td><td>31.08<br/>[14.6-53.13]</td><td>103.21<br/>[98-117.5]</td><td>45.82<br/>[26.67-68.78]</td><td>108.38<br/>[98.01-126.73]</td><td>65.5<br/>[40.78-94.2]</td><td></td><td></td></tr><tr><td colspan="13">TK02 clock model (Part B)</td></tr><tr><th>Data type</th><th>Models</th><th>Sampling fraction</th><th>Calibrations</th><th>Number of calibrations</th><th>Crown-Embioptera</th><th>Andensembiida + Embioptera</th><th>Anensembiidae</th><th>Archemiidae</th><th>Australembiidae</th><th>Embiidae</th><th>Clothodidae</th><th>Embiidae</th><th>Notoligotomidae</th><th>Oligotomidae</th><th>Ptilocerembiidae</th><th>Scelembiidae</th><th>Teratembiidae</th></tr><tr><td rowspan="10">Molecular + Morphological</td><td>GTR+I+G ≈ 0.1</td><td>Stem</td><td>3</td><td>172.63<br/>[147.55-199.42]</td><td>118.33<br/>[88.17-150.24]</td><td>85.16<br/>[63.53-108.95]</td><td>88.9<br/>[60.58-119.46]</td><td>121.72<br/>[87.91-169.21]</td><td>67.37<br/>[42.71-93.92]</td><td>106.42<br/>[83.82-131.14]</td><td>33.65<br/>[19.85-51.57]</td><td>123.85<br/>[92.69-156.8]</td><td>80.97<br/>[57.11-105.75]</td><td>129.05<br/>[105.52-154.26]</td><td>107.89<br/>[78.24-138.73]</td><td></td><td></td></tr><tr><td>GTR+I+G ≈ 0.21</td><td>Stem</td><td>3</td><td>170.91<br/>[145.01-197.35]</td><td>121.45<br/>[92.72-152.67]</td><td>87.27<br/>[65.69-110.68]</td><td>91.54<br/>[64.17-122.61]</td><td>132.86<br/>[90.95-172.54]</td><td>66.11<br/>[42.88-90.85]</td><td>104.77<br/>[20.96-53.92]</td><td>35.08<br/>[10.98-29.29]</td><td>117.32<br/>[87.59-150.4]</td><td>78.06<br/>[54.44-102.06]</td><td>127.11<br/>[103.8-152.33]</td><td>101.58<br/>[73.56-132.43]</td><td></td><td></td></tr><tr><td>RJMCMC ≈ 0.1</td><td>Stem</td><td>3</td><td>173.34<br/>[147.65-198.92]</td><td>116.44<br/>[85.86-143.5]</td><td>84.82<br/>[62.48-105.62]</td><td>87.44<br/>[58.89-114.17]</td><td>116.07<br/>[83.2-105.62]</td><td>67.14<br/>[43.21-92.28]</td><td>107.67<br/>[84.72-130.84]</td><td>33.41<br/>[20.42-50.36]</td><td>129.22<br/>[97.1-156.22]</td><td>85.55<br/>[59-107.69]</td><td>130.69<br/>[105.81-154.29]</td><td>112.63<br/>[83.69-140]</td><td></td><td></td></tr><tr><td>RJMCMC ≈ 0.21</td><td>Stem</td><td>3</td><td>170.02<br/>[144.99-194.53]</td><td>116.27<br/>[87.83-143.31]</td><td>84.56<br/>[63.06-106.27]</td><td>87.31<br/>[61.27-114.03]</td><td>117.19<br/>[85.3-163.46]</td><td>64.98<br/>[42.14-89.39]</td><td>104.38<br/>[81.2-126.5]</td><td>33.51<br/>[20.22-50.58]</td><td>122.88<br/>[92.07-150.91]</td><td>79.83<br/>[56.95-103.76]</td><td>126.95<br/>[103.17-150.18]</td><td>106.55<br/>[77.33-132.83]</td><td></td><td></td></tr><tr><td>GTR+I+G ≈ 0.1</td><td>Crown</td><td>5</td><td>172.61<br/>[148.81-196.13]</td><td>118.07<br/>[85.54-147.24]</td><td>84.92<br/>[64.42-107.92]</td><td>88.64<br/>[59.37-115.83]</td><td>119.73<br/>[88.76-167.23]</td><td>67.35<br/>[44.2-91.74]</td><td>106.84<br/>[86.11-129.38]</td><td>33.52<br/>[20.29-50.32]</td><td>125.25<br/>[99.1-150.51]</td><td>81.49<br/>[59.62-104.41]</td><td>129.65<br/>[108.02-151.86]</td><td>108.51<br/>[83.31-135.68]</td><td></td><td></td></tr><tr><td>GTR+I+G ≈ 0.21</td><td>Crown</td><td>5</td><td>172.57<br/>[150.45-197.21]</td><td>117.72<br/>[88.87-147.25]</td><td>85.17<br/>[63.89-107.27]</td><td>88.36<br/>[60.11-115.53]</td><td>117.57<br/>[87.77-163.21]</td><td>66.54<br/>[43.47-91.89]</td><td>106.34<br/>[19.86-50.67]</td><td>33.66<br/>[98.05-150.73]</td><td>126.01<br/>[58.96-105.49]</td><td>81.69<br/>[58.96-105.49]</td><td>129.04<br/>[107.25-151.32]</td><td>109.59<br/>[82.24-136.14]</td><td></td><td></td></tr><tr><td>RJMCMC ≈ 0.1</td><td>Crown</td><td>5</td><td>172.78<br/>[150.2-197.35]</td><td>118.62<br/>[90.86-147.48]</td><td>85.78<br/>[64.47-107.8]</td><td>88.82<br/>[62.45-116.49]</td><td>120.36<br/>[89.55-164.9]</td><td>67.25<br/>[43.23-90.81]</td><td>106.78<br/>[85.61-128.81]</td><td>34.19<br/>[20.65-51.97]</td><td>125.62<br/>[98-150.15]</td><td>81.97<br/>[59.14-105.69]</td><td>129.55<br/>[108.22-136.21]</td><td>109.03<br/>[83.81-136.21]</td><td></td><td></td></tr><tr><td>RJMCMC ≈ 0.21</td><td>Crown</td><td>5</td><td>171.36<br/>[148.32-195.04]</td><td>115.4<br/>[85-146.31]</td><td>83.81<br/>[63.14-107.72]</td><td>86.49<br/>[68.4-119.38]</td><td>114.63<br/>[84.43-160.31]</td><td>66.05<br/>[43.74-91.74]</td><td>105.51<br/>[83.71-127.17]</td><td>33.28<br/>[20.18-50.63]</td><td>124.69<br/>[99.47-149.29]</td><td>80.8<br/>[58.89-103.97]</td><td>128.16<br/>[107.07-151.53]</td><td>108.6<br/>[84.65-134.56]</td><td></td><td></td></tr><tr><td>GTR+I+G/Mkv+G ≈ 0.1</td><td>Stem</td><td>3</td><td>167.1<br/>[142.59-192.13]</td><td>NA</td><td>82.79<br/>[61.4-105.07]</td><td>92.74<br/>[63.12-108.16]</td><td>131.4<br/>[81.8-167.06]</td><td>61.17<br/>[40.9-84.41]</td><td>98.42<br/>[77.64-120.89]</td><td>29.62<br/>[18.48-44.23]</td><td>104.98<br/>[79.04-135.49]</td><td>69.78<br/>[48.69-91.93]</td><td>119.73<br/>[97.37-141.72]</td><td>91.42<br/>[68.21-122.18]</td><td></td><td></td></tr><tr><td>GTR+I+G/Mkv+G ≈ 0.21</td><td>Stem</td><td>3</td><td>167.35<br/>[143.21-191.6]</td><td>NA</td><td>86.3<br/>[66.12-108.16]</td><td>97.02<br/>[73.12-122.07]</td><td>140.85<br/>[111.05-168.74]</td><td>61.14<br/>[40.05-82.87]</td><td>97.34<br/>[76.79-117.54]</td><td>31.02<br/>[18.94-45.53]</td><td>99.27<br/>[79.44-119.7]</td><td>66.99<br/>[48.04-86.79]</td><td>118.09<br/>[97.65-139.83]</td><td>85.99<br/>[67.25-105.77]</td><td></td><td></td></tr><tr><td>RJMCMC/Mkv+G ≈ 0.1</td><td>Stem</td><td>3</td><td>167.2<br/>[142.11-191.67]</td><td>NA</td><td>83.8<br/>[63.09-105.82]</td><td>94.69<br/>[69.67-121.24]</td><td>138.31<br/>[90.82-170.09]</td><td>61.08<br/>[38.46-82.79]</td><td>97.82<br/>[77.43-119.2]</td><td>30.47<br/>[18.71-45.72]</td><td>101.14<br/>[78.98-130.75]</td><td>68.15<br/>[49.58-89.82]</td><td>118.64<br/>[97.73-141.22]</td><td>87.87<br/>[66.46-114.71]</td><td></td><td></td></tr><tr><td>RJMCMC/Mkv+G ≈ 0.21</td><td>Stem</td><td>3</td><td>166.78<br/>[141.55-191.36]</td><td>NA</td><td>85.73<br/>[64.92-108.6]</td><td>96.29<br/>[72.06-121.36]</td><td>140.35<br/>[110.61-168.78]</td><td>60.21<br/>[39.58-82.84]</td><td>96.49<br/>[76.6-118.68]</td><td>31.13<br/>[19.27-46.51]</td><td>98.42<br/>[78.72-120.72]</td><td>66.2<br/>[48.31-85.79]</td><td>117.18<br/>[96.2-139.8]</td><td>85.04<br/>[65.82-105.91]</td><td></td><td></td></tr></tbody></table>

(Continued)
<table><caption>Table 1. (Continued)</caption><caption>TK02 clock model (Part B)</caption><thead><tr><th>Data type</th><th>Models</th><th>Sampling fraction</th><th>Calibrations</th><th>Number of Crown-Embionta</th><th>Andeensembli + Archeembli</th><th>Anisembli</th><th>Archeembli</th><th>Australembli</th><th>Clothodli</th><th>Embiide</th><th>Notoligomide</th><th>Oligomide</th><th>Pilosecreembli</th><th>Sclerembli</th><th>Teremlihne</th></tr></thead><tbody><tr><td rowspan="7">Molecular +<br>Morphological</td><td>GTR+I+G/Mk+G ≈ 0.1</td><td>Crown</td><td>5</td><td>173.61<br>[153.87-193.9]</td><td>NA</td><td>89.4<br>[69.22-110.19]</td><td>100.24<br>[75.64-124.29]</td><td>146.26<br>[103.14-174.46]</td><td>65.08<br>[43.83-86.97]</td><td>102.46<br>[85.3-122.01]</td><td>31.93<br>[20.09-48.16]</td><td>106.6<br>[98-124.45]</td><td>72.61<br>[54.97-90.72]</td><td>124.14<br>[107.89-142.81]</td><td>92.72<br>[76.28-111.86]</td></tr><tr><td>GTR+I+G/Mk+G ≈ 0.21</td><td>Crown</td><td>5</td><td>173.95<br>[153.85-194.91]</td><td>NA</td><td>89.86<br>[69.93-110.24]</td><td>101.2<br>[76.21-125.53]</td><td>146.6<br>[100.51-176.76]</td><td>65.12<br>[44.23-87.68]</td><td>102.86<br>[85.18-122.91]</td><td>32.33<br>[19.92-47.87]</td><td>106.58<br>[98-125.4]</td><td>72.01<br>[54.44-91.38]</td><td>124.48<br>[107.17-143.11]</td><td>92.6<br>[76.9-113.18]</td></tr><tr><td>RJMCMC/Mk+G ≈ 0.1</td><td>Crown</td><td>5</td><td>174.28<br>[154.36-194.85]</td><td>NA</td><td>89.85<br>[70.24-110.84]</td><td>101.45<br>[76.38-124.41]</td><td>148.07<br>[123.32-174.04]</td><td>66.5<br>[45.48-88.05]</td><td>103.88<br>[86.22-123.41]</td><td>32.39<br>[20.03-48.13]</td><td>107.02<br>[98-123.38]</td><td>73.08<br>[55.64-91.95]</td><td>125.25<br>[109.1-144.67]</td><td>93.03<br>[77.74-111.53]</td></tr><tr><td>RJMCMC/Mk+G ≈ 0.21</td><td>Crown</td><td>5</td><td>172.8<br>[153.38-194.4]</td><td>NA</td><td>88.53<br>[68.56-110.31]</td><td>100.03<br>[74.99-123.21]</td><td>144.79<br>[89.46-171.74]</td><td>63.56<br>[43.39-85.03]</td><td>101.64<br>[83.57-120.47]</td><td>32.03<br>[20.07-47.42]</td><td>106.12<br>[98-124.82]</td><td>71.11<br>[53.23-89.55]</td><td>123.23<br>[107.07-141.68]</td><td>92.11<br>[76.31-113.25]</td></tr><tr><td>GTR+I+G ≈ 0.1</td><td>Variable</td><td>7</td><td>174.63<br>[151.74-197.63]</td><td>119.09<br>[93.63-148.12]</td><td>86.31<br>[68.71-107.53]</td><td>89.32<br>[64.11-117.84]</td><td>120.04<br>[91.35-163.92]</td><td>68.63<br>[45.51-94.24]</td><td>108.87<br>[88.23-130.41]</td><td>34.66<br>[21.68-51.43]</td><td>128.77<br>[101.7-154.49]</td><td>84.46<br>[62.8-107.44]</td><td>131.88<br>[111.79-154.07]</td><td>112.13<br>[86.8-138.21]</td></tr><tr><td>GTR+I+G ≈ 0.21</td><td>Variable</td><td>7</td><td>172.58<br>[151.1-196.18]</td><td>121.2<br>[94.69-149.04]</td><td>87.16<br>[67.77-108.55]</td><td>90.17<br>[63.47-117.82]</td><td>127.18<br>[92.16-169.61]</td><td>67.35<br>[44.53-91.99]</td><td>106.23<br>[85.45-128.79]</td><td>35.35<br>[21.54-52.67]</td><td>121.02<br>[98-148.8]</td><td>80.28<br>[57.69-103.79]</td><td>128.99<br>[108.63-151.57]</td><td>105.27<br>[80.65-134.4]</td></tr><tr><td>RJMCMC ≈ 0.1</td><td>Variable</td><td>7</td><td>172.77<br>[149.83-195.57]</td><td>117.17<br>[92.37-146.79]</td><td>84.89<br>[66.84-104.77]</td><td>87.2<br>[60.8-115.72]</td><td>117.93<br>[90.14-163.21]</td><td>67.18<br>[44.04-91.08]</td><td>106.87<br>[86.19-128.32]</td><td>34.12<br>[21.31-51.11]</td><td>125.66<br>[98.96-150.87]</td><td>82.18<br>[60.54-107.17]</td><td>129.74<br>[109.24-151.25]</td><td>109.69<br>[83.23-135.55]</td></tr><tr><td>RJMCMC ≈ 0.21</td><td>Variable</td><td>7</td><td>173.35<br>[150.01-196.37]</td><td>119.75<br>[93.78-148.04]</td><td>86.41<br>[67.11-108.62]</td><td>89.88<br>[63.49-117.61]</td><td>121.55<br>[90.28-167.66]</td><td>67.34<br>[43.81-91.91]</td><td>106.91<br>[85.98-129.18]</td><td>35.02<br>[20.86-52.52]</td><td>124.23<br>[98-150.06]</td><td>81.42<br>[59.23-105.88]</td><td>129.99<br>[109.15-152.38]</td><td>107.89<br>[82.25-135.7]</td></tr><tr><td rowspan="7">Molecular +<br>Morphological</td><td>GTR+I+G/Mk+G ≈ 0.1</td><td>Variable</td><td>7</td><td>174.34<br>[155.33-195.31]</td><td>NA</td><td>89.44<br>[70.65-110.61]</td><td>102.04<br>[78.8-125.47]</td><td>147.77<br>[122.2-173.44]</td><td>65.35<br>[44.07-87.04]</td><td>102.96<br>[85.28-122.32]</td><td>32.32<br>[19.94-48]</td><td>106.63<br>[98-122.91]</td><td>72.25<br>[54.9-91.35]</td><td>124.7<br>[108.3-143.18]</td><td>92.72<br>[77.26-110.57]</td></tr><tr><td>GTR+I+G/Mk+G ≈ 0.21</td><td>Variable</td><td>7</td><td>173.93<br>[156.19-193.8]</td><td>NA</td><td>90.47<br>[71.46-110.38]</td><td>102.23<br>[79.95-123.97]</td><td>147.21<br>[123.32-171.57]</td><td>64.77<br>[44.46-85.61]</td><td>102.19<br>[86.07-121.7]</td><td>32.72<br>[20.34-47.62]</td><td>105.61<br>[98-119.73]</td><td>71.38<br>[54.28-88.82]</td><td>123.84<br>[108.14-147.77]</td><td>91.45<br>[76.63-107.45]</td></tr><tr><td>RJMCMC/Mk+G ≈ 0.1</td><td>Variable</td><td>7</td><td>174.53<br>[155.84-195.54]</td><td>NA</td><td>90.66<br>[70.18-112.24]</td><td>102.85<br>[79.67-125.7]</td><td>148.13<br>[123.32-173.27]</td><td>65.95<br>[44.98-88.56]</td><td>103.72<br>[86.76-123.07]</td><td>32.97<br>[20.9-48.84]</td><td>106.99<br>[98-123.13]</td><td>72.62<br>[54.76-91.32]</td><td>125.17<br>[109.61-143.95]</td><td>92.99<br>[77.29-110.93]</td></tr><tr><td>RJMCMC/Mk+G ≈ 0.21</td><td>Variable</td><td>7</td><td>173.7<br>[154.66-194.73]</td><td>NA</td><td>90.15<br>[70.64-111.75]</td><td>101.24<br>[77.3-123.64]</td><td>145.7<br>[94.4-172.35]</td><td>64.71<br>[44.53-87.64]</td><td>102.92<br>[84.82-121.82]</td><td>32.23<br>[20.24-47.14]</td><td>107.2<br>[98-126.7]</td><td>72.31<br>[54.04-90.8]</td><td>124.46<br>[108.01-142.92]</td><td>93.03<br>[76.23-115.87]</td></tr></tbody></table>
clock: birthdeath. Sampling strategy of taxa was set to diversity (prset samplestrat = diversity) wherein taxa are sampled to maximize diversity. An exponential prior and a beta prior were used for the net speciation rate and the relative extinction rate, respectively [prset speciationpr = exp(10); prset extinctionpr = beta(1,1)].

To test for the impact of parameters and calibrations on divergence time estimates, we performed 48 analyses. We assigned calibrations with uniform priors to outgroup taxa: (1) Dictyoptera [minimum age: 245 Mya, maximum age: 315 Mya; corresponding to several occurrences of crown-Blattodea from the Triassic from Vosges (France) and elsewhere (Cariglino et al. 2020), and to the upper limit of the age estimates for the appearance of Dictyoptera following Legendre et al. (2015) and to several fossils from the Carboniferous]; (2) Phasmatodea [minimum age: 165 Mya, maximum age: 251.9 Mya; corresponding to the oldest unambiguous crown-Phasmatodea and to an age derived from Tihelka et al. (2020) for the appearance of crown-Phasmatodea]. We assigned uniform priors to the root of the tree [minimum age: 315 Mya, maximum age: 380 Mya; corresponding to several fossils from the Carboniferous, and to the age estimated for Polyneoptera in comprehensive studies (Misof et al. 2014, Tong et al. 2015)].

For the ingroup, we used different calibration approaches to explore the sensitivity of our dating estimates. First, we calibrated five nodes with uniform priors (Supporting information, Table S1): (1) crown Embiidae with Galloembia raholai (Oise amber; De Franceschi and De Ploëg 2003, Nel et al. 2004, Falières et al. 2021); (2) crown Scelembiidae with several species in Kachin amber (Engel and Grimaldi 2006, Anisyutkin and Perkovsky 2022, Liu et al. 2024) with plesiomorphic wing venations but body characters resembling species of crown Scelembiidae; (3) crown Teratembiidae with a fossil species of the genus Oligembia from the Miocene Dominican amber (Szumik 1994, Iturralde-Vinent and MacPhee 1996, Iturralde-Vinent 2001, Penney 2010); (4) crown Oligotomidae with Litoclostes delicatus (from Kachin amber, c. 98 Mya: Cruickshank and Ko 2003, Shi et al. 2012, Yu et al. 2019) presenting a generalized oligotomid morphology (Engel et al. 2016); and (5) the divergence between Embioptera and Phasmatodea using an unpublished specimen of stem Embioptera from the Capitanian (c. 260 Mya; Yao et al. 2015, Zhang et al. 2019). Although not yet formally described, this fossil only slightly pre-dates two other stem-Embioptera species belonging to the genus Alexarasnia, which were described from the Wuchipingian (254.14–259.51 Mya) (Gorochov 2011, Aristov 2017). It is worth noting that this calibration, although slightly older, aligns well with the discovery of an Embioptera fossil from the Triassic of the Mount San Giorgio Lagerstätte (240 Mya; Montagna et al. 2019, fig. 1f). This fossil exhibits enlarged basal foretarsomeres, a trait found in the crown Embioptera; but information on its wing venation and most of its body parts is lacking. Therefore, we refrain from considering it as definitive evidence for the crown Embioptera and instead consider that it calibrates the same node as the Capitanian fossil (i.e. Embioptera + Phasmatodea). The use of the unpublished Capitanian fossil as calibration point will bolster the longevity of our findings. Subsequently, we examined the impact of treating the calibrations as ‘stem calibrations’ on our results (Supporting information, Table S1). Specifically, we calibrated, with uniform priors, the node (Oligotomidae + Teratembiidae) using Litoclostes delicatus, the node (Scelembiidae

+ Embiidae) using one of the Scelembiidae species from the same Kachin amber deposit, and the divergence between Embioptera and Phasmatodea using the Capitanian fossil. Lastly, we calibrated seven nodes with uniform priors, to assess the impact of increasing the number of ‘stem calibrations’ on our estimates (Supporting information, Table S1): the same five as in the first set of analysis and two additional calibrations: (1) the node C2 with a Clothodidae (Henoclothoda simplex Cui & Engel, 2020) from Kachin amber (Cui et al. 2020); and (2) the node Notoligotomidae + its sister group (i.e., either Archembiidae or Anisembiidae) with a Notoligotomidae from Kachin amber (Burmitembia venosa Cockerell, 1919).

Henoclothoda simplex is classified within the crown group Clothodidae, specifically in the subfamily Clothodinae (Cui et al. 2020), and it could potentially be used to calibrate the divergence between Clothodinae and its sister group. However, our knowledge of the family's internal relationships remains limited, and our sampling of Clothodidae includes only two species, representing a small fraction of its crown group. Given these constraints, we have chosen to use Henoclothoda simplex to calibrate the node preceding the crown Clothodidae (i.e. C2), as this represents the most cautious and methodologically sound approach with the available data. The species Burmitembia venosa is classified in its own notoligotomid subfamily, the Burmitembiinae (Cockerell 1919, Engel and Grimaldi 2006). Therefore, this fossil cannot be used to calibrate the crown group Notoligotominae, which encompasses only one extant genus Notoligotoma, and it can only be used here to calibrate the node preceding the crown Notoligotominae (Engel and Grimaldi 2006, Engel et al. 2016). The detail of the combination of priors and calibrations is summarized in Table 1 and justifications for the choice of fossils are available in Supporting information, Table S1. This approach provides a valuable framework for exploring the uncertainty in fossil taxon placement and its impact on divergence time estimates (Klopfstein 2021).

MrBayes analyses comprised two runs and four Markov chain Monte Carlo (MCMC) for 50 to 150 million generations (until convergence was reached). MCMC were sampled every 5000 generations, and a burn-in fraction of 0.25 was used. Convergence diagnostics were checked for each analysis with the average standard deviation of split frequencies (< 0.01), the potential reduction scale factors (PRSF; close to 1.0), and effective sample size (ESS; > 200) in Tracer 1.7.2 (Rambaut et al. 2018). All consensus trees were visualized and drawn using FigTree 1.4.4 (Rambaut 2009) and modified with Adobe Illustrator CC2019. Files are available in the Figshare data repository (Jouault 2025b).

## Ancestral range estimation

To estimate ancestral ranges, we determined five or six areas for the current distribution of extant Embioptera species following Wallace's zoogeographical regions by Holt et al. (2013): New World (divided into the Neotropics and Nearctic in the second analysis), Palearctic, Australian (including Oceanian), Oriental, and Afrotropical. The definition and choice of these regions are designed to allow comparisons with the results of similar analyses for the Phasmatodea (Simon et al. 2019). For invasive species, only the area of origin is considered in the analyses. We conducted
the analysis using the Dispersal-Extinction-Cladogenesis model (DEC) implemented in the R package BioGeoBEARS 1.1.3 (Matzke 2013, 2018) based on two dated phylogenetic trees reconstructed using the best fitting models, and considering the most probable placement of fossils, along with a sampling fraction that reflects our current understanding of Embioptera diversity (1: GTR+I+G, crown calibrations, and SF = 0.21; 2: GTR+I+G and Mkv+G, crown calibrations, and SF = 0.21). We allowed a maximum of two areas to compose a range. Additionally, we fitted the DEC+J model (Matzke 2014) but refrained from comparing it to DEC due to concerns regarding statistical validity and model choice in BioGeoBEARS (Ree and Sanmartín 2018). We opted for DEC over DEC+J as the latter often indicates null or exceedingly low extinction rates (Sanmartín and Meseguer 2016). Furthermore, the DEC+J model tends to favour direct dispersal over widespread ranges, a preference that may not be suitable for reconstructing the history of old groups (Rolland and Condamine 2019). Files are available in the Figshare data repository (Jouault 2024).

## RESULTS

Our undated tree inferred with IQ-TREE (Supporting information, Fig. S1) revealed rate heterogeneity across the phylogeny. Consequently, we considered the results obtained using the IGR clock model to be more appropriate than those from the autocorrelated TK02 model and therefore focussed primarily on discussing the IGR results. Divergence time estimates from the TK02 analyses are provided in Table 1B.

### Family-group relationships among Embioptera

Deep nodes and inter-familial relationships consistently received maximal or relatively good support (PP ≥ 0.90; UFBS ≥ 80; Figs 2–3; Supporting Information, Fig. S1) with few exceptions observed in the clade comprising Anisembiidae, Notoligotomidae, Archembiidae, and Andesembiidae. Notably, the support for the clade (Anisembiidae + Notoligotomidae) was weak using molecular data only with both Bayesian (PP = 0.69) and ML analyses

(UFBS = 59). Similarly, in Bayesian analyses incorporating molecular and morphological data, the clade ((Andesembiidae + (Notoligotomidae + Archembiidae)) showed modest support (PP = 0.73), and in ML analyses, the relationships within the clade ((Archembiidae + (Notoligotomidae + Anisembiidae)) lacked strong support (Fig. 3; Supporting Information, Fig. S1).

The monophyly of Embioptera families examined in our analyses was robustly supported, receiving maximal support in Bayesian analyses (PP = 1). While some families did not reach maximum support in ML analyses, they still had high support, with the Archembiidae (UFBS = 96), the Scelembiidae (UFBS = 94), the Oligotomidae (UFBS = 99), and to a lesser extent, the Australembiidae (UFBS = 70), receiving UFBS scores superior or equal to 70.

We consistently found the Australembiidae as the sister lineage to the rest of Embioptera (PP = 1; UFBS = 100), and this position was retrieved in all analyses. We found three main Embioptera lineages (Figs 2–3). The first one (C1) encompassed (Anisembiidae + Notoligotomidae + Archembiidae + Andesembiidae) and received good support in Bayesian (PP = 1) and ML analyses (UFBS = 96). The second clade (C2) encompassed (Scelembiidae + Embiidae + Clothodidae) and received maximal support in Bayesian analyses (PP = 1) but was not supported in ML analyses (UFBS = 66). The third clade (C3) encompassed (Oligotomidae + Teratembiidae + Ptilocerembiidae + *Oedembia*) and was well supported by Bayesian inference (PP = 0.98, 0.9), but not supported by ML analyses (UFBS = 74) (Figs 2–3; Supporting Information, Fig. S1). It is worth mentioning that in some TK02 analyses, based on molecular data only, the clade C3 was divided into two clades (Ptilocerembiidae + *Oedembia*) and (Oligotomidae + Teratembiidae).

### Divergence time estimates for major embiopteran lineages

We estimated divergence times for and within Embioptera (Fig. 3). Using node-dating analyses and molecular data (GTR+I+G and a sampling fraction of ≈0.21), we estimated the divergence between Embioptera and Phasmatodea in the Cisuralian (Permian) around 283.66 Mya [95% HPD = 264.27, 298.88 Mya] and the crown Embioptera in the Middle Jurassic around 169.9 Mya

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of family-group relationships among Embioptera from each of the three methods used in this study (Bayesian topologies obtained with an IGR clock model).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Dated phylogenetic tree of webspinners based on Bayesian analyses with crown calibrations, GTR+I+G and Mkv+G models, and a sampling fraction of 0.21. A, molecular data only. B, molecular and morphological data combined. Support values (coloured dots) and ages with the 95% HPD represented by blue bars are indicated for each node.</figcaption></figure>
[95% HPD = 144.05, 200.05 Mya]. Similar ages were found when morphological data were integrated into the analyses (GTR+I+G and Mkv+G, and a sampling fraction of ≈0.21), with the divergence between Embioptera and Phasmatodea also estimated during the Cisuralian 284.34 Mya [95% HPD = 264.24, 298.89 Mya] and the crown Embioptera in the Middle Jurassic around 174.82 Mya [95% HPD = 148.02, 204.76 Mya].

Using molecular data only (GTR+I+G and a sampling fraction of ≈0.21), the origin of C1 was estimated to 141.27 Mya [95% HPD = 115.62, 169.56 Mya]. The origin of C2 was estimated to 135.98 Mya [95% HPD = 115.51, 159.8 Mya]. The origin of C3 was estimated to 136.02 Mya [95% HPD = 115.9, 159.72 Mya]. Likewise, using molecular data combined with morphological data (GTR+I+G and Mkv+G, and a sampling fraction of ≈0.21), the origin of C1 was estimated to 145.4 Mya [95% HPD = 117.02, 172.14 Mya]. The origin of C2 was estimated to 138.59 Mya [95% HPD = 116.69, 162.31 Mya]. The origin of C3 was estimated to 138.57 Mya [95% HPD = 117.14, 162.2 Mya]. The divergence time estimates for the crown Embioptera and webspinner families, for all the analyses and the different dating scenarios, are summarized in Table 1. They showed a marginal effect of the clock model choice, sampling fraction, the number of fossil calibrations, and the inclusion of a morphological partition; the effect of the stem vs. crown calibrations was somewhat more pronounced but remained insignificant at this time scale.

## Historical biogeography

All DEC analyses were consistent; minor variations affected only a few internal nodes. Our four analyses of ancestral range estimation (ARE), with the DEC model, at the root of the crown Embioptera were equivocally reconstructed to be New World (Nearctic + Neotropical) and Australia (Fig. 4). The ARE of crown Australi- an (Fig. 4; Supporting Information, Figs S2, S3, S4, S5). The AREs of the clades C1 and C2 were unequivocally reconstructed to be in the New World (Fig. 4; Supporting Information, Figs S2, S3, S4, S5). The ARE of the third embiopteran lineage (C3) was equivocally reconstructed to be in the New World and Oriental (Fig. 4; Supporting Information, Figs S2, S3, S4, S5).

## DISCUSSION

### Family-group relationships among Embioptera

We focus our discussion mostly on familial relationships, as Miller *et al.* (2012) have already greatly discussed the relationships within embiopteran families. Our analyses have provided a more resolved topology of Embioptera by identifying only monophyletic groups. We also obtained a similar topology for C2 and C3 in ML and Bayesian inferences (Fig. 2; Supporting Information, Fig. S1), whereas the topologies obtained by Miller *et al.* (2012) differed. These disparities, despite the use of similar datasets, can be attributed to the use of different models. Our methodology tends to eschew over-parameterization and presents an alternative approach in terms of partitioning, while offering improved compensation for site heterogeneity. For instance, in previous ML analyses, nine partitions—corresponding to the genes and their codon positions—were implemented (Miller *et al.* 2012), whereas our study employed only four after estimating the best-fit partitioning scheme. Additionally, we determined the best-fitting model for each partition, maximizing the compatibility between the data and the model (Kalyaanamoorthy *et al.* 2017), by incorporating parameters that account for rate heterogeneity among sites and variable discrete rate categories. In contrast, Miller *et al.* (2012) used a GTR-MIX model with a fixed number of discrete rate categories (i.e. four) for their ML analysis.

Our results strongly confirm the monophyly of the Australi-emiidae, positioning them as the most anciently diverging lineage within Embioptera (PP = 1, UFBS = 100; Fig. 3). This topology supports the results of Miller *et al.* (2012). The subsequent clade (C1) marks a new radiation of Embioptera and includes families that have often been considered closely related, such as Archemiidae and Notoligotomidae (Figs 2–3; Supporting Information, Fig. S1). Previous analyses relying on morphological data have indicated the proximity between these families (Szumik *et al.* 2008). When morphological data were combined with molecular data, the positioning of Anisembiidae and Andesembiidae varied, aligning with the inconsistent placement observed in previous analyses (e.g. Szumik 1996, Szumik *et al.* 2008, Miller *et al.* 2012). As our ML and Bayesian analyses, relying on either molecular data alone or a combination of molecular and morphological data, resulted in distinct topologies for C1 (Figs 2–3; Supporting Information, Fig. S1), selecting a definitive topology becomes impossible. We anticipate that the incorporation of new sequences and the broadening of available molecular data for the order will contribute to resolving this puzzle.

In clade C2, the close relationship between Embiidae and Scelembiidae was previously emphasized in morphological analyses, as representatives of Scelembiidae were initially considered part of Embiidae (Szumik 1996, fig. 9). However, the newfound association of Clothodidae with these families is a recent discovery (Miller *et al.* 2012) and contradicts their previous placement as the sister lineage to the rest of Embioptera, as indicated by morphological analyses (Szumik 1996, Szumik *et al.* 2008, Miller *et al.* 2012). Due to their seemingly generalized morphology in the head, wings, and male genitalia, members of Clothodidae have often been considered the sister group to the remaining Embioptera (Davis 1940, Ross 1970, 1987, Szumik 1996, Szumik *et al.* 2008). However, the placement of Clothodidae as a sister to the remaining Embioptera taxa has relied on authoritative assumptions about the polarity of certain characters without thorough testing (discussed in Miller *et al.* 2012).

The clade C3 includes three families (Ptiloceremiidae, Oligotomidae, and Teratembiidae) (Figs 2–3; Supporting Information, Fig. S1). Oligotomidae and Teratembiidae are consistently grouped together by morphological and molecular data (Szumik 1996, Szumik *et al.* 2008). More recent Bayesian analyses also support the grouping of Oligotomidae and Teratembiidae (Miller *et al.* 2012). However, morphological-based analyses have sometimes suggested that Ptiloceremiidae are closely related to either Notoligotomidae (Szumik 1996) or Archemiidae (Szumik *et al.* 2008). These discrepancies are maybe the results of plesiomorphic characters. The genus *Oedembia*,
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Ancestral range estimates of Embioptera as inferred with the DEC model. The nodal pie charts show the relative probabilities of the geographic ranges according to the in-figure colour code. Teratembiidae (Terat.), Ptilocerembiidae (Ptiloc.), Clothodidae (C.), Andesembiidae (A), Anisembiidae (Anisembiid.), Archembiidae (Arc), Quaternary (Q). Palaeomaps used with permission © 2020 Colorado Plateau Geosystems Inc.</figcaption></figure>

currently classified within the Embiidae, is the only genus of this
family found outside the clade and is resolved in a clade with
*Ptilocerembia* (Fig. 3). This topology is consistent with previous
analyses, in which the authors considered the possibility of
expanding the definition of Ptilocerembiidae to include *Oedem-
bia* (Miller *et al.* 2012). However, the formal transfer of this
genus to Ptilocerembiidae was deferred due to the absence of

unequivocal synapomorphies for this clade and the limited sam-
pling of other *Oedembia* and closely related taxa (Miller *et al.
2012*). While we concur with this earlier study, the maximum
support observed for this topology (UFBS = 100) in our analy-
ses led us to remove *Oedembia* from Embiidae, thus rendering
Embiidae monophyletic in our analyses, and designating the
genus as having uncertain familial affinities.
## Origin and early diversification of Embioptera

Our time-calibrated analyses indicate that the divergence between the Embioptera and the Phasmatodea occurred during the Artinskian (283.5–290.1 Mya), approximately 284 Mya (Fig. 3; Table 1). This result aligns with recent time-calibrated phylogenies of Phasmatodea, derived from genomic datasets, indicating their divergence from Embioptera between 260 Mya and 305 Mya (Tihelka et al. 2020). Notably, our analyses consistently support this result, with the youngest median age estimate at around 281 Mya and the oldest at approximately 284.5 Mya (Table 1). The radiation of stem Embioptera likely occurred during the Permian, preceding the extinction events of that period (Jouault et al. 2022a). However, the scarcity of their fossil record during the Permian and Triassic hampers an in-depth exploration of their diversity dynamics during these periods. It can be postulated that the concealed subterranean lifestyle of webspinners facilitated their survival during these periods.

Our divergence time analyses, employing fossils as crown calibrations, yielded consistent age estimates for the crown Embioptera, suggesting an origin in the Middle Jurassic or the latest Early Jurassic, approximately between 169 Mya and 175 Mya (Table 1). These results imply that ages derived from Phasmatodea phylogenies produced estimates that were too recent to capture the true timing of the origin of the crown Embioptera. Interestingly, the use of fossils as stem calibrations leads to slightly younger ages for the crown Embioptera, suggesting an origin in the early Late Jurassic, between 156 and 160 Mya (Table 1). This result is not necessarily surprising because when using stem calibrations, we reduced the number of calibrations in the in-group from four to two. Regarding the ages of embiopteran families, there is an approximate difference of ≈10 Mya between analyses employing stem calibrations and those employing crown calibrations, which confirms that our calibrations are extremely conservative (Table 1).

Interestingly, our results suggest a radiation of the crown Embioptera during the Cretaceous, with seven of the eleven families included in the analyses originating during the Cretaceous. The Teratembidae would have originated during the Palaeocene, the Clothodidae during the Early Eocene, the Ptilocerembiidae during the Middle Eocene, and the Notoligotomidae during the Oligocene (Fig. 3). This diversification of the crown Embioptera coincides with the Angiosperm Terrestrial Revolution (ATR), an event spanning 100 to 50 Mya (Benton et al. 2022). The ATR is characterized by the rapid diversification of the angiosperms during the Cretaceous, which likely participated in the decline of conifers (Condamine et al. 2020). The diversification of the flowering plants is assumed to lead to the rise of mutualistic relationships (pollination, herbivory) with insects, and is often mentioned as the trigger of the diversification of the ‘big five’ insect clades (e.g. McKenna et al. 2019, Kawahara et al. 2023, Jouault et al. 2024, Peris and Condamine 2024; but see Asar et al. 2022). While this event is acknowledged for its role in the diversification of certain insect groups, it might have also induced ecological changes leading to the decline of other lineages (e.g. Sinitshenkova 2003, Jouault et al. 2022b, Sroka et al. 2023, Jouault et al. 2024).

Considering the cryptic lifestyle of embiopterans, we propose that the diversification of angiosperms created novel ecological niches, food sources, and habitats, especially litter-type habitats,

thereby triggering the diversification of embiopterans, a scenario also proposed for ants and Blattodea (Evangelista et al. 2019, Jouault et al. 2024). Future studies should investigate this hypothesis, particularly by augmenting molecular data for the group and expanding the taxonomic sampling.

## Impact of morphological partition, priors, and model choices on time divergence estimates

It is noteworthy that in all our analyses, the inclusion of morphological data had little impact on the age estimates for the crown Embioptera and webspinner families (Fig. 3; Table 1). A previous study (Barba-Montoya et al. 2021) proposed that node calibrations, particularly when integrating molecular and morphological data, were more effective for nodes with numerous taxa—a methodology we adopted here. Similarly, they demonstrated a high correlation between time divergence estimates inferred using molecular data and those inferred from discrete morphological data from extant taxa (Barba-Montoya et al. 2021). Our results corroborate these empirical tests and suggest that their results can be generalized.

We also investigated how altering the sampling fraction (i.e. 0.21 vs. 0.1) affects time divergence estimates. It was demonstrated that the value chosen for the sampling fraction can have strong effects on dating estimates (Welch et al. 2005, Beaulieu et al. 2015). However, similarly to the inclusion of a morphological partition, reducing the sampling fraction by a factor of two had no or a reduced impact on time divergence estimates (Table 1). In most Bayesian analyses, taxon sampling is assumed to be random concerning phylogenetic relationships and is informed for the entire tree (referred to as global sampling fraction) (Bromham et al. 2018). While this assumption is violated in many datasets, especially those designed to include representatives of all higher taxa regardless of differences in species richness (e.g. one species *per* family), our dataset shows a relatively representative sampling of each family's diversity (e.g. Beaulieu et al. 2015, Bromham et al. 2018). Families with more genera and species, such as Oligotomidae (six genera) and Teratembidae (five genera), are represented by more tips in the topology compared to families with fewer genera and species, like Andesembiidae (two genera). We attribute the minimal fluctuations in time divergence estimates to this sampling strategy, underscoring the importance of maintaining this balance in future studies of the evolutionary history of Embioptera.

While tree topology exhibits low sensitivity to model misspecification, other parameters may be affected (Alfaro and Huelsenbeck 2006). Given that our study primarily focuses on time divergence estimates, we investigated how model selection could influence age estimates. We replicated all analyses conducted under the GTR+I+G substitution model using an RJMCMC model. Interestingly, the 16 analyses, which used comparable parameters, produced similar time divergence estimates (Table 1). These findings suggest that within the GTR family, model choice may have minimal impact on time divergence estimates in a node-dating context.

Future studies should evaluate how incorporating fossil taxa as tips in the topology influences divergence time estimates for Embioptera families, using tip-dating or total-evidence dating approaches. Additionally, improving the molecular coverage and increasing the number of species represented in the tree (aiming for a sampling
fraction of 1) will be crucial for gaining a better understanding of the timing and tempo of Embioptera diversification.

## Origin and early historical biogeography

Our AREs, derived either from molecular data-based or molecular and morphological data-based time-calibrated trees, support several vicariant events in the deep time and dispersal events during the Cenozoic. We found a New World (i.e. Neotropics) + Australian origin of webspinners in the Jurassic (Fig. 4; Supporting Information, Figs S2, S3, S4, S5), compatible with a Gondwanan origin for the crown Embioptera. It is also worth highlighting that our analyses with the DEC model inferred an ancient vicariance event between Australia and the New World (i.e. Neotropics).

Since the Jurassic, the crown Australembiidae likely inhabited eastern Gondwana, particularly what is now Australia (Fig. 4). We interpret the current distribution of these webspinners to be a consequence of a vicariant event, where an ancestral lineage thrived in the eastern Gondwana (comprising Australia and Antarctica) until the subsequent separation of the landmass during the Cretaceous and Early Cenozoic periods (Müller et al. 2000, McLoughlin 2001, Sanmartín and Ronquist 2004, Seton et al. 2012).

The C1 clade likely originated in Gondwana, encompassing lineages from the Neotropics, Nearctic, and Australia (Fig. 4; Supporting Information, Figs S2, S3, S4, S5). Our analysis of AREs points to a New World + Australian origin, for this clade, with the Notoligotomidae representing an Australian radiation (Fig. 4). Considering the divergence age between the Notoligotomidae and other lineages within C1 (i.e. Early Cretaceous), it is plausible that this family colonized Australia before the complete Gondwanan break-up, following a narrative akin to the Australembiidae and other insects (e.g. Almeida et al. 2012, Vicente et al. 2017, Letsch et al. 2021). The Archembiidae are today distributed in the Neotropics and the central eastern part of Africa. This distribution likely resulted from the split, during the break-up of Gondwana, of an ancestral lineage distributed in the central Gondwana. The present distribution of Anisembiidae likely resulted from a subsequent dispersal event from the Neotropics to the Nearctic during the Palaeogene (Rota et al. 2016).

The C2 clade has an unambiguous origin estimated in the New World, specifically in the Neotropics. All Scelembiidae and Clothodidae included in our analyses are thought to have originated from the Neotropics (Fig. 4; Supporting Information, Figs S2, S3, S4, S5). Within C2, the origin of Embiidae is ambiguous, possibly Palearctic or Afrotropical, although our various analyses tend to favour a Palearctic origin during the Early Cretaceous (Supporting Information, Figs S2, S3, S4, S5). While our analyses alone make the distribution challenging to interpret, the fossil record may provide additional insight into this pattern. An Eocene North American fossil Embiidae suggests their likely presence in this area before the Eocene (Ross 1984). It is plausible that crown Embiidae experienced declines in specific areas and underwent dispersal events along migration routes during the Palaeocene and Eocene periods, yet further data are required to comprehensively unravel their biogeographic history (e.g. Archibald et al. 2006, Archibald and Makarkin 2006, Brikiatis 2014).

Clade C3 has an ambiguous origin, estimated to be in the New World + Oriental regions (Fig. 4; Supporting Information,

Figs S2, S3, S4, S5). Our AREs unequivocally indicate a New World (i.e. Neotropical) origin for the Teratembiidae, aligning with previous studies and suggesting that Nearctic Teratembiidae may have followed a dispersal road comparable to Anisembiidae (Szumik et al. 2008). Their presence on the western coast of Africa lets us assume that they were ancestrally present in central Gondwana and later split during the Gondwana break-up (i.e. vicariance).

The Ptilocerembiidae + *Oedembia* and Oligotomidae trace their origin in the Oriental region during the Early Cretaceous, emerging from an origin in the New World + Oriental regions. Considering the Neotropical origin of the Teratembiidae and the putative Gondwanan origin of all the clades preceding the Ptilocerembiidae + *Oedembia* (Fig. 4; Supporting Information, Figs S2, S3, S4, S5), it is plausible that the Ptilocerembiidae + *Oedembia* and the Oligotomidae originated from insular India in the Cretaceous and Palaeogene before further diversification ('out-of-India' hypothesis) (e.g. Bossuyt and Milinkovitch 2001, Karanth 2006). This scenario aligns not only with the timing of the AREs but also with the fossil record (Engel et al. 2016). Following the collision of India and Asia, the Ptilocerembiidae likely expanded in the Oriental region, and some Oligotomidae dispersed to Australia, then to the Nearctic via the Beringian land bridge, and subsequently to the Palearctic through the De Geer or Thulean routes (Condamine et al. 2013, fig. 1).

## CONCLUSION

Our study presents a first time-calibrated phylogeny for the order Embioptera, 'harmonizing' previously conflicting topologies. Our analyses converge on an estimated Upper or Middle Jurassic origin for the crown Embioptera, with a split from the Phasmatodea during the Permian. These results are consistent with the fossil record, which includes extinct families and putative stem-group representatives (such as Sinembiidae and Alexarasniidae), whose systematic affinities have yet to be investigated within a phylogenetic framework. To assess the robustness of our age estimates, we conducted many dating analyses employing two different clock models, as well as different substitution models, sampling fractions, and fossil calibrations, and the inclusion or exclusion of morphological data. Our findings indicate that incorporating morphological data and transitioning across different clock and substitution models minimally affected age estimates for the crown Embioptera and webspinner families. Similarly, our analyses demonstrate that halving the sampling fraction had either no impact or a marginal effect on time divergence estimates. By adopting a comprehensive approach, we have overcome many of the known biases associated with dating analyses.

However, our study comes with limitations. A major constraint lies in our sampling, which includes only 11 of the 18 currently recognized embiopteran families. Although we built upon the most extensive molecular dataset available, with a good balance between gene coverage and completeness, our analyses remain affected by incomplete lineage sorting, which likely contributes to some unresolved relationships in the phylogeny. Embioptera is a particularly challenging group to collect, with few specialists, complicating efforts to generate a fully comprehensive molecular
framework. Therefore, we advocate for the sampling and sequencing of new taxa and new molecular data, especially for the missing families. In the genomic era, it is worth noting that no whole genome has yet been produced for Embioptera (GenBank, last assessed 29 April 2025). All these efforts would provide a foundation for resolving uncertain nodes in the Embioptera tree of life, for deciphering their diversification, and for bringing attention to this understudied group of insects.

To unravel the biogeographic history of the order, we conducted ancestral range estimations, using topologies derived from both molecular and morphological data, as well as molecular data alone. Our results suggest a Gondwanan origin for the crown Embioptera, indicating that the distribution of extant Embioptera families stems from numerous vicariant events in ancient times, succeeded by dispersal events during the Cenozoic era. This is particularly evident when considering the Burmese amber biota, which provided many of the calibrations used here. Located near present-day South Africa during the Jurassic and most of the Cretaceous, this biota belongs to the Burmese Terrane, which later migrated and collided with Laurasia. This complex geological history likely created vicariance events, followed by episodes of dispersal once the biota came into contact with the continental landmasses. Integrating fossil species as tips in the phylogeny will be important to advance our understanding of the biogeographic history of the order, refine divergence time estimates, and clarify the boundaries of webspinner families. Our study provides a temporal framework for these insects and lays the foundation for future research on their biogeographic history, both of which are crucial for understanding the timing and diversification patterns of Embioptera evolution.

## ACKNOWLEDGEMENTS

We thank three anonymous reviewers and the associate editor for their constructive comments. We thank Enrico Bonino and Michele Baldi for letting us use their photographs of fossil Embioptera trapped in Kachin amber, and Prof. André Nel (MNHN) for discussions on the early version of this paper.

## AUTHOR CONTRIBUTIONS

Corentin Jouault: Conceptualization (equal); Investigation (lead); Writing—original draft (lead); Writing—review and editing (lead); Visualization (lead); Validation (lead); Methodology (lead); Data curation (equal); Software (lead); Formal analysis (lead); Funding acquisition (supporting); Project administration (lead); Resources (equal); Supervision (equal); Fabien L. Condamine: Conceptualization (lead); Investigation (equal); Writing—original draft (equal); Writing—review and editing (equal); Visualization (supporting); Validation (equal); Methodology (supporting); Data curation (supporting); Software (equal); Formal analysis (supporting); Funding acquisition (lead); Project administration (lead); Resources (equal); Supervision (lead); and Frédéric Legendre: Conceptualization (lead); Investigation (equal); Writing—original draft (equal); Writing—review and editing (equal); Visualization (supporting); Validation (equal);

Methodology (supporting); Data curation (supporting); Software (equal); Formal analysis (supporting); Funding acquisition (lead); Project administration (lead); Resources (equal); Supervision (lead).

## SUPPLEMENTARY DATA

Supplementary data is available at *Evolutionary Journal of the Linnean Society* online.

Figure S1. Maximum-likelihood phylogeny of the Embioptera. Topology obtained with IQ-TREE, with ultrafast bootstrap values written at nodes.

Figure S2. Ancestral range estimates for Embioptera inferred using all models implemented in BioGeoBEARS, based on a molecular phylogeny and a five-area biogeographic scheme.

Figure S3. Ancestral range estimates for Embioptera inferred using all models implemented in BioGeoBEARS, based on a phylogeny built with morphological and molecular data and a five-area biogeographic scheme.

Figure S4. Ancestral range estimates for Embioptera inferred using all models implemented in BioGeoBEARS, based on a molecular phylogeny and a six-area biogeographic scheme.

Figure S5. Ancestral range estimates for Embioptera inferred using all models implemented in BioGeoBEARS, based on a phylogeny built with morphological and molecular data and a six-area biogeographic scheme.

Table S1. Calibration strategies, fossil selection, and justifications for their phylogenetic placement.

## CONFLICT OF INTEREST

The authors declare no conflicts of interest.

## FUNDING

This project benefitted from the ‘Investissements d’Avenir’ programme managed by the Agence Nationale de la Recherche (CEBA, ref. ANR-10-LABX-25-01).

## DATA AVAILABILITY

The data that support the findings of this study are openly available. Sequences used in the analyses are available in GenBank (https://www.ncbi.nlm.nih.gov/genbank/; Miller *et al.* 2012). Alignments used for tree reconstruction analyses and all input files are provided in a Figshare data repositories (Jouault 2024, Jouault 2025a, b).

## REFERENCES

Aberlenc H-P, Albouy V, Barthélémy D, Beaucournu J-C, Blandin P, Cliquennois N *et al.* *Les Insectes du Monde. Biodiversité. Classification. Clés de détermination des familles*. Versailles, Montpellier & Plaissan: Quae & Museo éditions, 2020, Tome 1, pp. 1192; Tome 2, pp. 656.

Alfaro ME, Huelsenbeck JP. Comparative performance of Bayesian and AIC-based measures of phylogenetic model uncertainty. *Systematic Biology* 2006;55:89–96. https://doi.org/10.1080/10635150500433565
Almeida EAB, Pie MR, Brady SG *et al*. Biogeography and diversification of colletid bees (Hymenoptera: Colletidae): emerging patterns from the southern end of the world. *Journal of Biogeography* 2012;39:526–44. https://doi.org/10.1111/j.1365-2699.2011.02624.x

Anisyutkin LN, Perkovsky EE. New data on embiids (Insecta: Embiodea) from Mid-Cretaceous Burmese amber, with description of new genus and two new species. *Cretaceous Research* 2022;134:105149. https://doi.org/10.1016/j.cretres.2022.105149

Archibald SB, Cover SP, Moreau CS. Bulldog ants of the Eocene Okanagan Highlands and history of the subfamily (Hymenoptera: Formicidae: Myrmeciinae). *Annals of the Entomological Society of America* 2006;99:487–523. https://doi.org/10.1603/0013-8746(2006)99[487:BAOTEO]2.0.CO;2

Archibald SB, Makarkin VN. Tertiary giant lacewings (Neurotera: Polystoechotidae): revision and description of new taxa from Western North America and Denmark. *Journal of Systematic Palaeontology* 2006;4:119–55. https://doi.org/10.1017/S1477201906001817

Aristov DS. New insects (Insecta: Eoblattida, Embiida) from the Permian of Russia and the Triassic of Kyrgyzstan, with observations on the origin of webspinners. *Paleontological Journal* 2017;51:161–70. https://doi.org/10.1134/S0031030117020046

Asar Y, Ho SYW, Sauquet H. Early diversifications of angiosperms and their insect pollinators: were they unlinked? *Trends in Plant Science* 2022;27:858–69. https://doi.org/10.1016/j.tplants.2022.04.004

Barba-Montoya J, Tao Q, Kumar S. Molecular and morphological clocks for estimating evolutionary divergence times. *BMC Ecology and Evolution* 2021;21:83. https://doi.org/10.1186/s12862-021-01798-6

Beaulieu JM, O'Meara BC, Crane P *et al*. Heterogeneous rates of molecular evolution and diversification could explain the Triassic age estimate for angiosperms. *Systematic Biology* 2015;64:869–78. https://doi.org/10.1093/sysbio/syv027

Benton MJ, Wilf P, Sauquet H. The Angiosperm Terrestrial Revolution and the origins of modern biodiversity. *New Phytologist* 2022;233:2017–35. https://doi.org/10.1111/nph.17822

Bossuyt F, Milinkovitch MC. Amphibians as indicators of Early Tertiary "out-of-India" dispersal of vertebrates. *Science* 2001;292:93–5. https://doi.org/10.1126/science.1058875

Brikiatis L. The De Geer, Thulean and Beringia routes: key concepts for understanding Early Cenozoic biogeography. *Journal of Biogeography* 2014;41:1036–54. https://doi.org/10.1111/jbi.12310

Bromham L, Duchêne S, Hua X *et al*. Bayesian molecular dating: opening up the black box. *Biological Reviews* 2018;93:1165–91. https://doi.org/10.1111/brv.12390

Cariglino B, Lara MB, Zavattieri AM. Earliest record of fossil insect oothecae confirms the presence of crown-dictyopteran taxa in the Late Triassic. *Systematic Entomology* 2020;45:935–47. https://doi.org/10.1111/syen.12442

Cockerell TDA. Two interesting insects in Burmese amber. *The Entomologist* 1919;52:193–5.

Condamine FL, Silvestro D, Koppelhus EB *et al*. The rise of angiosperms pushed conifers to decline during global cooling. *Proceedings of the National Academy of Sciences of the United States of America* 2020; 117:28867–75. https://doi.org/10.1073/pnas.200557111

Condamine FL, Sperling FAH, Kergoat GJ. Global biogeographical pattern of swallowtail diversification demonstrates alternative colonization routes in the Northern and Southern Hemisphere. *Journal of Biogeography* 2013;40:9–23. https://doi.org/10.1111/j.1365-2699.2012.02787.x

Cui Y, Chen ZT, Engel MS. New species of webspinners (Insecta: Embiodea) from Mid-Cretaceous amber of northern Myanmar. *Cretaceous Research* 2020;113:104457. https://doi.org/10.1016/j.cretres.2020.104457

Cruickshank RD, Ko K. Geology of an amber locality in the Hukawng Valley, northern Myanmar. *Journal of Asian Earth Sciences* 2003;21:441–55. https://doi.org/10.1016/S1367-9120(02)00044-5

Davis C. Family classification of the order Embioptera. *Annals of the Entomological Society of America* 1940;33:677–82.

De Franceschi D, De Ploëg G. Origine de l'ambre des faciès sparnaciens (Éocène inférieur) du Bassin de Paris: le bois de l'arbre producteur. *Geodiversitas* 2003;25:633–47.

Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. *Nucleic Acids Research* 2004;32:1792–7. https://doi.org/10.1093/nar/gkh340

Edgerly JS, Life beneath silk walls: a review of the primitively social Embiidae. In: Choe JC, Crespi B (eds), *The Evolution of Social Behavior in Insects and Arachnids*. Cambridge, UK: Cambridge University Press, 1997, 14–25.

Engel MS, Grimaldi DA. The earliest webspinners (Insecta: Embiodea). *American Museum Novitates* 2006;3514:1–15.

Engel MS, Grimaldi DA, Singh H *et al*. Webspinners in Early Eocene amber from western India (Insecta, Embiodea). *ZooKeys* 2011;148:197–208. https://doi.org/10.3897/zookeys.148.1712

Engel MS, Huang DY, Breitkreuz LCV *et al*. Two new species of Mid-Cretaceous webspinners in amber from northern Myanmar (Embiodea: Clothodidae, Oligotomidae). *Cretaceous Research* 2016;58:118–24. https://doi.org/10.1016/j.cretres.2015.10.007

Evangelista DA, Wipfler B, Béthoux O *et al*. An integrative phylogenomic approach illuminates the evolutionary history of cockroaches and termites (Blattodea). *Proceedings of the Royal Society B: Biological Sciences* 2019;286:20182076. https://doi.org/10.1098/rspb.2018.2076

Falières E, Engel MS, Nel A. Earliest occurrence of Embiidae: a new genus from Earliest Eocene Oise amber (Insecta: Embiodea). *Comptes Rendus Palévol* 2021;20:799–805. https://doi.org/10.5852/cr-palevol2021v20a38

Gorochov AV. A new, enigmatic family for new genus and species of Polyneoptera from the Upper Permian of Russia. *ZooKeys* 2011;130:131–6. https://doi.org/10.3897/zookeys.130.1487

Grimaldi DA, Engel MS. *Evolution of the Insects*. Cambridge, England, United-Kingdom: Cambridge University Press, 2005.

Heath TA, Hedtke SM, Hillis DM. Taxon sampling and the accuracy of phylogenetic analyses. *Journal of Systematics and Evolution* 2008;46:239–57. https://doi.org/10.3724/SP.J.1002.2008.08016

Holt BG, Lessard J-P, Borregaard MK *et al*. An update of Wallace's zoo geographic regions of the world. *Science* 2013;339:74–8. https://doi.org/10.1126/science.1228282

Hopkins HEL. Embioptera species file, 1900. Retrieved on 25 March 2025 from: https://embioptera.speciesfile.org/otus/915589/overview

Huelsenbeck JP, Larget B, Alfaro ME. Bayesian phylogenetic model selection using reversible jump Markov chain Monte Carlo. *Molecular Biology and Evolution* 2004;21:1123–33. https://doi.org/10.1093/molbev/msh123

Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogeny. *Bioinformatics* 2001;17:754–5. https://doi.org/10.1093/bioinformatics/17.8.754

Hugall AF, Lee MSY. The likelihood node density effect and consequences for evolutionary studies of molecular rates. *Evolution* 2007;61:2293–307. https://doi.org/10.1111/j.1558-5646.2007.00188.x

Hoang DT, Chernomor O, von Haeseler A *et al*. UFBoot2: improving the ultrafast bootstrap approximation. *Molecular Biology and Evolution* 2018;35:518–22. https://doi.org/10.1093/molbev/msx281

Iturralde-Vinent MA. Geology of the amber-bearing deposits of the Greater Antilles. *Caribbean Journal of Science* 2001;37:141–67.

Iturralde-Vinent MA, MacPhee RDE. Age and paleogeographical origin of Dominican amber. *Science* 1996;273:1850–2. https://doi.org/10.1126/science.273.5283.1850

Jouault C. Supporting Files for Likelihood analyses. figshare. Dataset. 2025a. https://doi.org/10.6084/m9 figshare.29037446.v1

Jouault C. Supporting Files for Dating analyses. figshare. Dataset. 2025b. https://doi.org/10.6084/m9 figshare.29037437.v1

Jouault C, Condamine FL, Legendre F, Perrichot V. The Angiosperm Terrestrial Revolution buffered ants against extinction. *Proceedings of the National Academy of Sciences of the United States of America* 2024;121:e2317795121. https://doi.org/10.1073/pnas.2317795121

Jouault C, Coutret B, Konhauser KO *et al*. New odonatans (Odonata: Gomphaeschnidae; Synlestidae) from the Paleocene Paskapoo Formation: systematic and biogeographical implications. *Journal of Systematic Palaeontology* 2023;21:2261457. https://doi.org/10.1080/14772019.2023.2261457
Jouault C, Nel A, Legendre F *et al.* Estimating the drivers of diversification of stoneflies through time and the limits of their fossil record. *Insect Systematics and Diversity* 2022b;6:1–14. https://doi.org/10.1093/isd/ ixac017

Jouault C, Nel A, Perrichot V *et al.* Multiple drivers and lineage-specific insect extinctions during the Permo–Triassic. *Nature Communications* 2022a;13:7512. https://doi.org/10.1038/s41467-022-35284-4

Jouault C, Oyama N, Álvarez-Parra S *et al.* The radiation of Hymenoptera illuminated by Bayesian inferences from the fossil record. *Current Biology* 2025;35:2164–74.e4. https://doi.org/10.1016/j.cub.2025.03.002

Kalyaanamoorthy S, Minh BQ, Wong TKF *et al.* ModelFinder: fast model selection for accurate phylogenetic estimates. *Nature Methods* 2017;14:587–9. https://doi.org/10.1038/nmeth.4285

Karanth KP. Out-of-India Gondwanan origin of some tropical Asian biota. *Current Science* 2006;90:789–92.

Kawahara AY, Storer C, Carvalho APS *et al.* A global phylogeny of butterflies reveals their evolutionary history, ancestral hosts and biogeographic origins. *Nature Ecology Evolution* 2023;7:903–13. https://doi.org/ 10.1038/s41559-023-02041-9

Klopstein S. The age of insects and the revival of the minimum age tree. *Austral Entomology* 2021;60:138–46. https://doi.org/10.1111/aen. 12478

Lanfear R, Frandsen PB, Wright AM *et al.* PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. *Molecular Biology and Evolution* 2017;34:772–3. https://doi.org/10.1093/molbev/msw260

Larsson A. AliView: a fast and lightweight alignment viewer and editor for large data sets. *Bioinformatics* 2014;30:3276–8. https://doi.org/ 10.1093/bioinformatics/btu531

Legendre F, Nel A, Svenson GJ, Robillard T, Pellens R, Grandcolas P. Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence. *PLoS One* 2015;10:e0130127. https://doi.org/10.1371/journal. pone.0130127

Lepage T, Bryant D, Philippe H *et al.* A general comparison of relaxed molecular clock models. *Molecular Biology and Evolution* 2007;24:2669–80. https://doi.org/10.1093/molbev/msm193

Letsch H, Simon S, Frandsen PB *et al.* Combining molecular datasets with strongly heterogeneous taxon coverage enlightens the peculiar biogeographic history of stoneflies (Insecta: Plecoptera). *Systematic Entomology* 2021;46:952–67. https://doi.org/10.1111/syen.12505

Lewis P. A likelihood approach to estimating phylogeny from discrete morphological character data. *Systematic Biology* 2001;50:913–25. https:// doi.org/10.1080/106351501753462876

Liu S, Peng Z, Shi C *et al.* New genus and species of webspinners (Insecta: Embioptera) from the Mid-Cretaceous of Myanmar with a catalog of fossil members. *Insects* 2024;15:636. https://doi.org/10.3390/insects 15090636

Luo A, Zhang C, Zhou Q-S *et al.* Impacts of taxon-sampling schemes on Bayesian tip dating under the fossilized birth-death process. *Systematic Biology* 2023;72:781–801. https://doi.org/10.1093/sysbio/syad011

Maddison WP, Maddison DR. Mesquite: a modular system for evolution-ary analysis. Version 3.61, 2019. http://www.mesquiteproject.org.

Maehr MD. SF embioptera: embioptera species file (version 5.0, Jun 2018). In: Roskov Y, Ower G, Orrell T, Nicolson D, Bailly N, Kirk PM, Bourgoin T, DeWalt RE, Decock W, Nieukerken E, Zarucchi J, Penev L (eds), *Species 2000 & ITIS Catalogue of Life, 2019 Annual Checklist*. Leiden, The Netherlands: Naturalis, Species 2000. Digital resource at www. catalogueoflife.org/annual-checklist/2019. ISSN 2405-884X, 2019.

Matzke NJ. Probabilistic historical biogeography: new models for found-er-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. *Frontiers of Biogeography* 2013;5:242–8. http://escholarship.org/uc/item/44j7n141

Matzke NJ. Model selection in historical biogeography reveals that found-er-event speciation is a crucial process in island clades. *Systematic Biol- ogy* 2014;63:951–70. https://doi.org/10.1093/sysbio/syu056

Matzke NJ. BioGeoBEARS: BioGeography with Bayesian (and likelihood) Evolutionary Analysis with R Scripts. version 1.1.3, 2018. https:// github.com/nmatzke/BioGeoBEARS

McKenna DD, Shin S, Ahrens D *et al.* The evolution and genomic basis of beetle diversity. *Proceedings of the National Academy of Sciences, USA* 2019;116:24729–37. https://doi.org/10.1073/pnas

McLoughlin S. The breakup history of Gondwana and its impact on pre-Ce- nozoic floristic provincialism. *Australian Journal of Botany* 2001;49:271– 300. https://doi.org/10.1071/BT00023

Miller KB, Hayashi C, Whiting MF *et al.* The phylogeny and classification of Embioptera (Insecta). *Systematic Entomology* 2012;37:550–70. https://doi.org/10.1111/j.1365-3113.2012.00628.x

Minh BQ, Schmidt HA, Chernomor O *et al.* IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. *Molecular Biology and Evolution* 2020;37:1530–4. https://doi.org/10.1093/ molbev/msaa015

Misof B, Liu S, Meusemann K *et al.* Phylogenomics resolves the timing and pattern of insect evolution. *Science* 2014;346:763–7. https://doi.org/ 10.1126/science.1257570

Montagna M, Tong KJ, Magoga G *et al.* Recalibration of the insect evolu-tionary time scale using Monte San Giorgio fossils suggests survival of key lineages through the End-Permian Extinction. *Proceedings of the Royal Society B* 2019;286:20191854. http://dx.doi.org/10.1098/rspb. 2019.1854

Morinaga G, Soghigian J, Edgerly JS. Macroecology and potential drivers of diversity in webspinner maternal care (order Embioptera). *Insect Systematics and Diversity* 2023;7:11. https://doi.org/10.1093/isd/ixac031

Müller RD, Gaina C, Tikku A, Mihut D, Cande S, Stock JM. Mesozoic/ Cenozoic tectonic events around Australia. In: Richards M, Gordon R (eds), *The History and Dynamics of Global Plate Motions*. American Geophysical Union Monograph, Vol. 121. Hoboken, New Jersey, U.S.: John Wiley & Sons, 2000, 161–188.

Nel A, de Ploëg G, Millet J *et al.* The French ambers: a general conspectus and the Lowermost Eocene amber deposit of Le Quesnoy in the Paris Basin. *Geologica Acta* 2004;2:3–8.

Peris D, Condamine FL. The angiosperm radiation played a dual role in the diversification of insects and insect pollinators. *Nature Communications* 2024;15:552. https://doi.org/10.1038/s41467-024-44784-4

Ree RH, Sanmartín I. Conceptual and statistical problems with the DEC+ J model of founder-event speciation and its comparison with DEC via model selection. *Journal of Biogeography* 2018;45:741–9. https://doi. org/10.1111/jbi.13173

Rambaut A. FigTree v1.4.4. 2009. http://tree.bio.ed.ac.uk/software/ figtree/ (10 January 2024, date last accessed).

Rambaut A, Drummond AJ, Xie D *et al.* Posterior summarisation in Bayes-ian phylogenetics using Tracer 1.7. *Systematic Biology* 2018;67:901–4. https://doi.org/10.1093/sysbio/syy032

Rolland J, Condamine FL. The contribution of temperature and continental fragmentation to amphibian diversification. *Journal of Biogeography* 2019;46:1857–73. https://doi.org/10.1111/jbi.13592

Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. *Bioinformatics* 2003;19:1572–4. https://doi. org/10.1093/bioinformatics/btg180

Ronquist F, Teslenko M, van der Mark P *et al.* MRBAYES 3.2: efficient Bayesian phylogenetic inference and model selection across a large model space. *Systematic Biology* 2012b;61:539–42. https://doi. org/10.1093/sysbio/sys029

Ross ES. Biosystematics of the Embioptera. *Annual Review of Entomology* 1970;15:157–72.

Ross ES. A synopsis of the Embiidina of the United States. *Proceedings of the Entomological Society of Washington* 1984;86:82–93.

Ross ES. Studies in the insect order Embiidina: a revision of the family Clotho-didae. *Proceedings of the California Academy of Sciences* 1987;45:9–34.

Ross ES. Embioptera—Embiidina (Embiids, webspinners, foot-spinners). In: Naumann ID, Carne PB, Lawrence JF *et al.* (eds), *The Insects of Australia, Chapter 26*. Melbourne: Melbourne University Press, 1991, 405–409.
Rota J, Peña C, Miller SE. The importance of long-distance dispersal and establishment events in small insects: historical biogeography of metalmark moths (Lepidoptera, Choreutidae). *Journal of Biogeography* 2016;43:1254–65. https://doi.org/10.1111/jbi.12721

Sanmartín I, Meseguer AS. Extinction in phylogenetics and biogeography: from timetrees to patterns of biotic assemblage. *Frontiers in Genetics* 2016;7:35. https://doi.org/10.3389/fgene.2016.00035

Sanmartín I, Ronquist F. Southern Hemisphere biogeography inferred by event-based models: plant versus animal patterns. *Systematic Biology* 2004;53:216–43. https://doi.org/10.1080/10635150490423430

Seton M, Müller RD, Zahirovic S *et al*. Global continental and ocean basin reconstructions since 200 Ma. *Earth-Science Reviews* 2012; **113**:212–70.
https://doi.org/10.1016/j.earscirev.2012.03.002

Shi GH, Grimaldi DA, Harlow GE *et al*. Age constraint on Burmese amber based on U-Pb dating of zircons. *Cretaceous Research* 2012; **37**:155–63.
https://doi.org/10.1016/j.cretres.2012.03.014

Simon S, Letsch H, Bank S *et al*. Old World and New World Phasmatodea: phylogenomics resolve the evolutionary history of stick and leaf insects. *Frontiers in Ecology and Evolution* 2019;7:345. https://doi.org/10.3389/fevo.2019.00345

Sinitshenkova ND. Main ecological events in aquatic insects history. *Acta Zoologica Cracoviensia* 2003; **46**:381–92.

Soares AER, Schrago CG. The influence of taxon sampling and tree shape on molecular dating: an empirical example from mammalian mitochondrial genomes. *Bioinformatics and Biology Insights* 2012; **6**:129–43.

Sroka P, Godunko RJ, Prokop J. Fluctuation in the diversity of mayflies (Insecta, Ephemerida) as documented in the fossil record. *Scientific Reports* 2023; **13**:16052. https://doi.org/10.1038/s41598-023-42571-7

Szumik CA. *Oligembia vetusta*, a new fossil teratembiid (Embioptera) from Dominican amber. *Journal of the New York Entomological Society* 1994; **102**:67–73.

Szumik CA. The higher classification of the order Embioptera: a cladistic analysis. *Cladistics* 1996; **12**:41–64.

Szumik CA, Edgerly JS, Hayashi C. Phylogeny of embiopterans (Insecta). *Cladistics* 2008; **24**:993–1005. https://doi.org/10.1111/j.1096-0031.2008.00228.x

Thorne JL, Kishino H. Divergence time and evolutionary rate estimation with multilocus data. *Systematic Biology* 2002; **51**:689–702. https://doi.org/10.1080/10635150290102456

Tihelka E, Cai C, Giacomelli M *et al*. Integrated phylogenomic and fossil evidence of stick and leaf insects (Phasmatodea) reveal a Permian–Triassic co-origination with insectivores. *Royal Society Open Science* 2020;7:201689. http://dx.doi.org/10.1098/rsos.201689

Tong JK, Duchêne S, Ho SYW *et al*. Comment on “Phylogenomics resolves the timing and pattern of insect evolution”. *Science* 2015; **349**:487. https://doi.org/10.1126/science.aaa5460

Vicente N, Kergoat GJ, Dong J *et al*. In and out of the Neotropics: historical biogeography of Eneopterinae crickets. *Journal of Biogeography* 2017; **44**:2199–210. https://doi.org/10.1111/jbi.13026

Welch JJ, Fontanillas E, Bromham L. Molecular dates for the “Cambrian explosion”: the influence of prior assumptions. *Systematic Biology* 2005; **54**:672–8. https://doi.org/10.1080/10635150590947212

Wipfler B, Letsch H, Frandsen PB *et al*. Evolutionary history of Polyneoptera and its implications for our understanding of early winged insects. *Proceedings of the National Academy of Sciences of the United States of America* 2019; **116**:3024–9. https://doi.org/10.1073/pnas.1817794116

Yao X, Zhou Y-Q, Hinnov L.A. Astronomical forcing of a Middle Permian chert sequence in Chaohu, South China. *Earth and Planetary Science* 2015; **422**:206–21. https://doi.org/10.1016/j.epsl.2015.04.017

Yu TT, Thomson U, Mu L *et al*. An ammonite trapped in Burmese amber. *Proceedings of the National Academy of Sciences* 2019; **116**:11345–50. https://doi.org/10.1073/pnas.1821292116

Zhang B-L, Yao S-P, Hu W-X *et al*. Development of a high-productivity and anoxic–euxinic condition during the Late Guadalupian in the Lower Yangtze region: implications for the Mid-Capitanian extinction event. *Palaeogeography, Palaeoclimatology, Palaeoecology* 2019; **531**:108630–16. https://doi.org/10.1016/j.palaeo.2018.01.021

Zhang C, Stadler T, Klopfstein S *et al*. Total-evidence dating under the fossilized birth–death process. *Systematic Biology* 2016; **65**:228–49. https://doi.org/10.1093/sysbio/syv080