Inocelliidae is one of the two extant families of the holometabolan order Raphidioptera (snakeflies), with the modern fauna represented by seven genera and 44 species. The evolutionary history of the family is little known. Here we present the first phylogenetic and biogeographic analyses based on a worldwide sampling of taxa and datasets combined with morphological characters and mitochondrial genomes, aiming to investigate the intergeneric phylogeny and historical biogeography of Inocelliidae. The phylogenetic inference from the combined analysis of morphological and molecular data recovered the sister-group relationship between a clade of (Negha + Indianoinocellia) + Sininocellia and a clade of Fibla + the Inocellia clade (interiorly nested by Amurinocellia and Parainocellia). Amurinocellia stat. rev. and Parainocellia stat. rev. et emend. nov. are relegated to subgeneric status within Inocellia, while a newly erected subgenus of Inocellia, Epinocellia subgen. nov., accommodates the former Parainocellia burmana (U. Aspöck and H. Aspöck, 1968) plus a new species Inocellia (Epinocellia) weii sp. nov. Further, the Inocellia crassicornis group constitutes the nominate subgenus Inocellia stat. nov., but the Inocellia fulvostigmata group is paraphyletic. Diversification within Inocelliidae is distinguished by an Eocene divergence leading to extant genera and a Miocene radiation of species. A biogeographic scenario depicts how the diverse inocelliid fauna from East Asia could have originated from western North America via dispersal across the Beringia during the early Tertiary, and how the Miocene ancestors of Inocellia could have accomplished long-distance dispersals via the Tibet‐Himalayan corridor or eastern Palaearctic to western Palaearctic. Our results shed new light specifically on the evolution of Inocelliidae and, in general, the Raphidioptera.

Data Matrix

File S1. Checklist of extant species of Inocelliidae.

File S2. Specimens examined for the present study.

File S3. Nexus file of Maximum parsimony tree based on morphological data.

File S4. TNT file of Maximum parsimony tree based on morphological data.

File S5. Nexus file of Bayesian inference tree based on dataset PCG123 (nucleotide data of 13 protein-coding genes).

File S6. Phylip file of Maximum likelihood tree based on dataset PCG123 (nucleotide data of 13 protein-coding genes).

File S7. Fasta file of Maximum parsimony tree based on dataset PCG123 (nucleotide data of 13 protein-coding genes).

File S8. Nexus file of Bayesian inference tree based on dataset PCG123rRNA (nucleotide data of the 13 protein-coding genes and rRNA genes).

File S9. Phylip file of Maximum likelihood tree based on dataset PCG123rRNA (nucleotide data of the 13 protein-coding genes and rRNA genes).

File S10. Fasta file of Maximum parsimony tree based on dataset PCG123rRNA (nucleotide data of the 13 protein-coding genes and rRNA genes).

File S11. Nexus file of Bayesian inference tree based on dataset PCG12 (nucleotide data of the 13 protein-coding genes without the third position of PCGs)

File S12. Phylip file of Maximum likelihood tree based on dataset PCG12 (nucleotide data of the 13 protein-coding genes without the third position of PCGs)

File S13. Fasta file of Maximum parsimony tree based on dataset PCG12 (nucleotide data of the 13 protein-coding genes without the third position of PCGs)

File S14. Nexus file of Bayesian inference tree based on dataset PCG_AA (amino acid sequences of the 13 protein-coding genes).

File S15. Phylip file of Maximum likelihood tree based on dataset PCG_AA (amino acid sequences of the 13 protein-coding genes).

File S16. Fasta file of Maximum parsimony tree based on dataset PCG_AA (amino acid sequences of the 13 protein-coding genes).

File S17. Nexus file containing the molecular and morphological alignments.

File S18. TNT file containing the molecular and morphological alignments.

Table S1. List of species sampled in the present phylogenetic analysis and GenBank accession number.

Table S2. The partitioning schemes and selected amino acid substitution models for phylogenetic analysis.

Table S3. Saturation test of the gene PCG and its different codon sites and gene rRNA.

Table S4. Divergence time estimates of Inocelliidae. Nodal code corresponds to that in Fig. 4. Time-scale units are in millions of years (Ma).

Table S5. Likelihood scores and model comparison for ancestral area reconstruction using BioGeoBEARS. Nodal code corresponds to that in Fig. 4.

Table S6. Major biogeographic events of Inocelliidae estimated with BioGeoBEARS. Nodal code corresponds to that in Fig. 4.

Table S7. Summary of ancestral area reconstructions based on DIVALIKE + J model. Node for nodes in Fig. 4.

Table S8. Beast median divergence times for major nodes in Fig. 4, with 95% credibility intervals and marginal likelihoods.

Fig. S1. Images of antennae of inocelliids. (a) Indianoinocellia mayana U. Aspöck, H. Aspöck, & H. Rausch; (b) Inocellia elegans Liu, H. Aspöck, Yang & U. Aspöck.

Fig. S2. Images of head and pronotum of snakeflies. (a) Amurinocellia sinensis (Navás), male; (b) Inocellia cheni Liu, H. Aspöck, Yang & U. Aspöck, male; (c) Inocellia elegans Liu, H. Aspöck, Yang & U. Aspöck, male; (d) Negha inflata (Hagen), male; (e) Sininocellia gigantos Yang, male; (f) Inocellia elegans Liu, H. Aspöck, Yang & U. Aspöck, female; (g) Amurinocellia sinensis (Navás), female; (h) Negha inflata (Hagen), female; (i) Agulla arizonica (Banks), female; (j) Sininocellia chikun Liu, H. Aspöck, Zhan, & U. Aspöck, female; (k) Fibla pasiphae (H. Aspöck & U. Aspöck), female; (l) Inocellia fulvostigmata fulvostigmata U. Aspöck & H. Aspöck, female; (m) Inocellia crassicornis (Schummel), female. Arrows indicate character state used in the phylogenetic analysis of the morphology-only dataset as indicated in the text.

Fig. S3. Images of wings of snakeflies. Arrows indicate character state used in the phylogenetic analysis of the morphology-only dataset as indicated in the text.

Fig. S4. Drawings of male fused gonocoxites 11 (gonarcus) of inocelliids. (a), (d), (e) and (g)–(i) modified from H. Aspöck et al. (1991); (b) modified from Shen et al. (2019); (f), (k) (l) and (o) modified from Liu et al. (2010a); (c) modified from Liu et al. (2010b); (j) modified from Liu & Hajong. (2015); (m) modified from Liu et al. (2012a); (n) modified from Liu et al. (2009b). Scale bar = 0.5 mm.

Fig. S5. Drawings of male fused gonocoxites 11 (gonarcus) of inocelliids. (a)–(n) modified from H. Aspöck et al. (1991); (o) modified from Liu et al. (2012b). Scale bar = 0.5 mm.

Fig. S6. Drawings of male fused gonocoxites 11 (gonarcus) of inocelliids. (e)–(h) modified from U. Aspöck et al. (2011); (b) modified from Liu et al. (2010a); (a) (j) (c) and (d) modified from Liu et al. (2010b); (i) modified from U. Aspöck et al. (2009); (k) modified from Liu et al. (2014a). Scale bar = 0.5 mm.

Fig. S7. Maximum likelihood tree based on PCG123 dataset.

Fig. S8. Bayesian inference tree based on PCG123 dataset.

Fig. S9. Maximum parsimony tree based on PCG123 dataset.

Fig. S10. Maximum likelihood tree based on PCG123rRNA dataset.

Fig. S11. Bayesian inference tree based on PCG123rRNA dataset.

Fig. S12. Maximum parsimony tree based on PCG123rRNA dataset.

Fig. S13. Maximum likelihood tree based on PCG12 dataset.

Fig. S14. Bayesian inference tree based on PCG12 dataset.

Fig. S15. Maximum parsimony tree based on PCG12 dataset.

Fig. S16. Bayesian inference tree based on PCG_AA dataset.

Fig. S17. Maximum parsimony tree based on PCG_AA dataset.

Fig. S18. Bayesian inference tree based on the combined analysis of morphological and molecular data. The genera in colored taxa are inocelliids, and the grey colored taxa are outgroups in Raphidiidae.

Fig. S19. Maximum parsimony tree based on the combined analysis of morphological and molecular data. The genera in colored taxa are inocelliids, and the grey colored taxa are outgroups in Raphidiidae.

Fig. S20. BEAST divergence times estimates for key nodes in inocelliid cladogenesis (node codes correspond to those in Fig. 4), under different priors Violin plots presenting mean ages (white dots) along with the posterior distribution (blue area) of the 95% credibility intervals (black bars) inferred in the different BEAST analyses. (a) node Ino4: divergence between Sininocellia and Negha, (b) node Ino8: crown Inocelliidae), (c) node Ino9: divergence between Fibla and Inocellia clade, (d) node Ino21: crown I. crassicornis group, (e) node Ino23: divergence between western Palaearctic and eastern Palaearctic I. crassicornis, (f) node Ino24: divergence between Parainocellia and Amurinocellia. YU, Yule Tree model; BD, birth-death Tree model; Exp, exponential prior distribution for fossil calibrations; LogN; lognormal prior distribution for fossil calibrations; Uni, uniform prior distribution for fossil calibrations.