Unraveling the evolutionary history of the snakefly family Inocelliidae (Insecta: Raphidioptera) through integrative phylogenetics

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 biogeographical 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.r. and Parainocellia stat.r. et emend.n. are relegated to subgeneric status within Inocellia, whereas a newly erected subgenus of Inocellia, Epinocellia subgen.n., accommodates the former Parainocellia burmana (U. Aspöck and H. Aspöck, 1968) plus a new species Inocellia (Epinocellia) weii sp.n. Further, the Inocellia crassicornis group constitutes the nominotypical subgenus Inocellia stat.n., 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 biogeographical 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.


Introduction
Raphidioptera (snakeflies) are the smallest order of the holometabolous insects, comprising c. 250 extant species in 33 genera worldwide Oswald, 2022). Currently, the order is placed as the sister group of Megaloptera + Neuroptera in a number of recent phylogenetic studies based on different types of data ranging from morphology to genomes (Asp€ ock et al., 2001;Misof et al., 2014;Vasilikopoulos et al., 2020;Wang et al., 2017;Winterton et al., 2018). Morphologically, adult snakeflies are characterized by a prognathous head, a slender elongate prothorax, the presence of a pigmented pterostigma and a long ovipositor. The larvae can be recognized by their elongate, flattened body, a prognathous head with chewing-mandibulate mouthparts and a soft 10-segmented abdomen Asp€ ock et al., 1991). Biologically, extant snakeflies are distinguished by an obligatory period with reduced (cooler) temperature for regular ontogenesis and their arboreal, predatory mode of life (Asp€ ock, 1998(Asp€ ock, , 2002. Although the taxonomy of extant snakeflies is well-founded owing to extensive studies over the past decades (Asp€ ock et al., 1991Liu et al., 2009aLiu et al., , 2009bLiu et al., , 2010aLiu et al., , 2010bLiu, Asp€ ock, Zhan, et al., 2012;Liu, Asp€ ock, Zhang, et al., 2012;Liu et al., 2013;Liu, Asp€ ock, et al., 2014;Liu & Hajong, 2015;Liu et al., 2018;Shen et al., 2019Shen et al., , 2021, several key questions on the evolution of snakeflies remain unresolved. One regards the largely undetermined phylogenetic relationships at the level of genus and species within the order. Because of this, little is known about the timescale of divergence and the biogeographical history which led to a modern range confined to the Northern Hemisphere (Asp€ ock, 2002;Haring et al., 2011). Early snakefly evolution is another key issue that remains elusive as a result of the vague relationships between the diverse fossil snakeflies and their modern counterparts (Bechly & Wolf-Schwenninger, 2011;Liu, Ren, et al., 2014;Makarkin & Archibald, 2014;Willmann, 1994).
Inocelliidae is one of the two extant families of Raphidioptera, with a modern fauna composed of seven genera and 44 species (Oswald, 2022) (File S1). Adult inocelliids are characterized by absence of ocelli, whereas Raphidiidae, the other extant snakefly family, have three distinct ocelli. Compared with Raphidiidae, Inocelliidae also can be distinguished by the absence of a veinlet within the pterostigma, the absence of a proximal veinlet in the pterostigma, shortened, funnel-like tarsomere 2 and shell-like male gonocoxite 9 (Asp€ ock et al., 1991). Nevertheless, many diagnostic characters of Inocelliidae are shared by some Mesozoic fossil snakeflies, and such similarity obscures the systematic position of this family in Raphidioptera (Liu et al., 2016;Lu et al., 2020;Lyu et al., 2020). As far as known, the larvae of all Inocelliidae are corticolous, which means that these insects require trees with a suitable bark for their development. Thus, in contrast to many Raphidiidae with terricolous larvae, Inocelliidae do not occur in arid areas. The development of Inocelliidae lasts at least two, however, often three (or even more) years.
The extant genera of Inocelliidae known before this work comprised Amurinocellia H. Asp€ ock and U. Asp€ ock, from East Asia (three species), Fibla Nav as from Mediterranean Europe and North Africa (four species), Indianoinocellia U. Asp€ ock and H. Asp€ ock from western North and Central America (two species), Inocellia Schneider which is widespread in Eurasia (26 species), Negha Nav as from western North America (three species), Parainocellia H. Asp€ ock & U. Asp€ ock from Mediterranean Europe and Southeast Asia (four species), and Sininocellia Yang from East Asia (two species) (Fig. 1). The majority of species (31 species in four genera) occur in the eastern Palaearctic and Oriental regions with Inocellia crassicornis (Schummel) distributed also in the western Palaearctic Liu et al., 2009a, 2010a, 2010b, Liu, Asp€ ock, Zhan, et al., 2012Liu, Asp€ ock, et al., 2014;Liu et al., 2018;Shen et al., 2019Shen et al., , 2021, whereas the remaining 13 species are endemic to either the western Palaearctic or western Nearctic regions (Asp€ ock et al., 1991). Based on fossil records, three extinct genera of inocelliids are known from the Eocene of Europe and North America, namely Electrinocellia Engel, Succinofibla U. Asp€ ock & H. Asp€ ock and Paraksenocellia Makarkin, Archibald & Jepson (Asp€ ock & Asp€ ock, 2004;Carpenter, 1958;Engel, 1998;Makarkin et al., 2019). The fossil record thus documents over 50 Myr of the evolutionary history of Inocelliidae. Nonetheless, molecular dating estimated a Cretaceous or even Jurassic origin of Inocelliidae (Vasilikopoulos et al., 2020;Wang et al., 2017;Winterton et al., 2018), although no definite inocelliid fossil has been reported from the Mesozoic.
However, their phylogenetic hypothesis was not recovered in our cladistic analysis. Subsequently, two phylogenetic analyses on the higher phylogeny of Neuropterida have peripherally provided new evidence to understand the intergeneric phylogeny of Inocelliidae (Vasilikopoulos et al., 2020;Winterton et al., 2018); both studies yielded the clade Negha + (Fibla + (Inocellia + Parainocellia)). Unfortunately, these results for the phylogeny within Inocelliidae are inconsistent, and none of the studies was based on a comprehensive sampling of taxa.
Here we present the first phylogenetic analysis of Inocelliidae using an integrative approach based on a comprehensive sampling and a combined dataset of morphological characters and mitochondrial genome sequences. We also estimate a timescale of cladogenesis within Inocelliidae based on molecular dating methods, and by combining the ancestral area reconstruction, we were able to trace the biogeographical history of this family. Our study provides a phylogenetic framework of Inocelliidae, which serves as the basis for a proposed revised classification. Moreover, the spatio-temporal pattern of divergence within Inocelliidae presently recovered sheds light on the evolutionary history of Raphidioptera.

Taxon sampling
All 44 described inocelliid species (File S1) of all seven extant genera were sampled in our phylogenetic analysis based on morphology.
A partial set in these taxa (six genera and 18 species) (Table S1), for which specimens were available for the molecular study, was sampled in the mitochondrial phylogenomic analysis: Amurinocellia H. Asp€ ock & U. Asp€ ock (two species), Fibla Nav as (two species), Inocellia Schneider (nine species), Negha Nav as (one species), Parainocellia H. Asp€ ock & U. Asp€ ock (two species) and Sininocellia Yang (two species). Only Indianoinocellia U. Asp€ ock & H. Asp€ ock was not included in the molecular phylogenetic analysis due to lack of material. In addition, a new species, Inocellia (Epinocellia) weii sp.n., following the revised classification presented here, was included in both the morphological and molecular analyses.

Morphological characters
Morphological characters used in the phylogenetic analysis comprise a total of 48 adult characters (see Appendix 1). Of these, 41 are coded as binary and seven as multistate characters. Unknown characters were coded as "?", whereas inapplicable characters were coded as "-". The character matrix was prepared in Mesquite v.3.61 (Maddison & Maddison, 2019). The selection of characters and assessment of character states was based on a comprehensive comparative morphological study. There are nine nongenital characters, and the remaining 39 characters were obtained from male and female genitalia. Illustrations are provided for most characters (Figs 7-9, Figs S1-S6), with some figures modified and adapted for the present study from the literature (Asp€ ock & Asp€ ock, 1968;Asp€ ock et al., 1982Asp€ ock et al., , 1991Asp€ ock et al., , 2011Liu et al., 2009aLiu et al., , 2009bLiu et al., , 2010aLiu et al., , 2010bLiu, Asp€ ock, Zhan, et al., 2012;Liu, Asp€ ock, Zhang, et al., 2012;Liu et al., 2013;Liu, Asp€ ock, et al., 2014;Liu & Hajong, 2015;Liu et al., 2018;Shen et al., 2019Shen et al., , 2021. Habitus photographs of adults were taken with Nikon D800 and Nikon D850 digital cameras with Nikon Micro Nikkor 105 mm lens. Photographs of genitalia were taken with a Nikon D850 on a Leica DM 2000 optical microscope. Illustrations of genitalia were handdrawn using a Nikon SMZ745 stereo microscope. Genital preparations were made by heating the apex of the abdomen in a 10% aqueous solution of potassium hydroxide (KOH) at 120°C on a hot plate for about 5-6 min. After rinsing the KOH with distilled water, the apex of the abdomen was transferred to glycerin for examination. Terminology of genital sclerites follows Asp€ ock and Asp€ ock (2008), whereas the terminology of wing venation generally follows

DNA extraction and mitochondrial genome sequencing
Total genomic DNA was extracted from the muscle tissue of meso-or metathorax or hind legs using TIANamp Genomic DNA Kit (Tiangen, China). The sequencing library was generated using Truseq Nano DNA HT Sample Prep Kit (Illumina, USA) with the size of 350 bp. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using Novaseq PE Cluster Kit v. 2.5 (Illumina). The DNA libraries were sequenced with Illumina Hiseq platform and generated 150-bp paired-end reads.

Genome assembly and annotation
The contigs were assembled in two methods: de novo assembly was performed with IDBA-UD v.1.1.1 (Peng et al., 2012). The assemblies were constructed using 200 for the setting of minimum size of contig, and an initial k-mer size of 20, an iteration size of 10 and a maximum k-mer size of 90. Reference-guided assembly was performed with Geneious v.11.1.4 (Kearse et al., 2012) (https://www. geneious.com) for further manual examination. The published mitogenome of Negha inflata (Hagen) (GenBank accession number: KT425086) was used as a reference in our study. The preliminary mitogenome annotations were conducted using the MITOS webserver, under default settings and the invertebrate genetic code for mitochondria (Bernt et al., 2013). The complete mitogenomes are deposited in GenBank (Table S1). For the mitochondrial proteincoding genes, we first removed the stop codon of each sequence, and then aligned nucleotide sequences based on their corresponding amino acid translations using the MAFFT algorithm implemented in TranslatorX (Abascal et al., 2010). For the ribosomal (r) RNA genes, each of them was aligned using MAFFT v.7.123 (Katoh & Standley, 2013) under the iterative refinement method incorporating the most accurate local pairwise alignment information (E-INS-i). Poorly aligned sites were removed with Gblocks v.0.91 (Talavera & Castresana, 2007). All alignments then were checked and corrected manually in MEGA v.7.0 (Kumar et al., 2016) and Mesquite v.3.61 (Maddison & Maddison, 2019). Finally, all fasta alignments were concatenated using FASconCAT_v. 1.0 (K€ uck & Meusemann, 2010) to construct a full dataset of protein-coding genes (PCGs) and rRNA (Files S7, S10, S13, S16).

Phylogenetic analyses
Three types of data were analysed: the morphology-only (see above), the molecule-only and the morphology + molecule data. The molecule-only data further comprised four different datasets: PCG123 (nucleotide data of 13 PCGs), PCG123rRNA (nucleotide data of the 13 PCGs and rRNA genes), PCG12 (nucleotide data of 13 PCGs without the third position of PCGs) and PCG_AA (amino acid sequences of the 13 PCGs). The morphology + molecule dataset was composed of the morphological data from 48 characters of 44 species and the PCG_AA dataset from 22 species. We used the PCG_AA results to estimate divergence times and ancestral areas because the tree topologies were congruent with previous results from larger genomic datasets (Vasilikopoulos et al., 2020;Winterton et al., 2018) (also see Results and Discussion below).
The analysis of the morphology-only dataset was performed using the maximum parsimony (MP) method in TNT v.1.5 (Goloboff & Catalano, 2016) with a traditional search (random seed = 1, replicates = 100, Tree-bisection-reconnection (TBR), 10 trees saved per replication). Bootstrap values were calculated with the function implemented in TNT (TBR from existing trees, retain trees suboptimal by 100 steps). All characters were treated as unordered and with equal weight. Only unambiguous character changes were mapped on the strict consensus tree. Sininocellia chikun (Liu, H. Asp€ ock, Zhan, & U. Asp€ ock) was omitted from the MP analysis as a consequence of insufficient data.
For the molecule-only dataset, nucleotide substitution saturation was tested for each gene partition using Xia's method (Xia et al., 2003) implemented in DAMBE v.7.0 (Xia, 2018). Phylogenetic trees were reconstructed under maximum-likelihood (ML), Bayesian inference (BI) and Maximum parsimony (MP). We used PartitionFinder v.2.1.1 (Lanfear et al., 2016) to determine the best partitioning schemes for the datasets under the Bayesian Information Criterion (BIC). ML analysis was performed in IQ-TREE (Nguyen et al., 2015;Trifinopoulos et al., 2016), with various data partition schemes and best-fitting models determined by PartitionFinder. Detailed information on the partitions and the best models selected (GTR-GAMMA model for nucleotide data and LG + GAMMA model for amino acid data) are summarized in Table S2. BI analysis was performed in MrBayes v.3.2.6 (Ronquist et al., 2012) implemented on XSEDE (Extreme Science and Engineering Discovery Environment) through the CIPRES Science Gateway (Miller et al., 2010) (http:// www.phylo.org/) portal, using default parameters, random starting trees and two independent runs. The GTR model for nucleotide and the WAG model for amino acid were chosen separately. The number of substitution types was set to 6. Markov Chain Monte Carlo (MCMC) chains were run for 2 000 000 generations and trees were sampled every 1000 generations. Analyses were terminated when the standard deviation of clade frequencies fell below 0.01, which indicated that stationarity had been reached. The first 25% of the obtained tree files were regarded as burn-in and the remaining trees were used to construct the majority-rule consensus tree, and then the posterior probabilities (PP) calculated. MP analysis was performed in TNT 1.5 (Goloboff & Catalano, 2016), with the same parameters as above in the analysis of the morphological data.
The morphology + molecule dataset was analysed using BI in MrBayes v.3.2.6 (Ronquist et al., 2012) and MP in TNT v.1.5 (Goloboff & Catalano, 2016). The amino acid and morphological data were separately partitioned by genes and by the number of character states (see Appendix 1). The morphological data were partitioned into three groups of characters, respectively with 41 characters (each with two states), six characters (each with three states) and one character (with four states). We performed a heuristic search using the 'search = user' option in PartitionFinder 2 (Lanfear et al., 2016) to select the best-fitting model for amino acid data, and using BI criteria for model selection, each partition was treated as a separate data block in PartitionFinder. The WAG model was selected for the amino acid data. All morphological characters were treated as unordered, the nst command was set to mixed, a distinct Lewis Mkv substitution model (Lewis, 2001) was applied to each partition separately, combined with gamma distributed rates (+G) with a shared shape parameter for each partition to account for variation in substitution rates of different morphological characters. Each analysis consisted of four MCMC chains run simultaneously for 10 million generations and trees were sampled every 1000 generations. We ensured that the average standard deviation of split frequencies fell below 0.01, which indicated that stationarity had been reached. The first 25% of the obtained tree files were regarded as burn-in and the remaining trees were used to construct the majority-rule consensus tree, and then the PP were calculated. For the MP analysis as well as the calculation of the bootstrap values, the parameters were set the same as the above analyses using TNT.

Divergence time estimation
Estimation of divergence times was performed in BEAST v.2.6.0 via the CIPRES web portal (Miller et al., 2010). We utilized the PCG_AA dataset and the topology resulting from the ML analysis. We defined partitions and site models that were identical to those used for the above phylogenetic analysis in BEAUti. The Yule or a birth death model of speciation was applied as tree prior, with timetree and clock model linked across partitions, estimated base frequencies, gamma model for site heterogeneity and four gamma categories. A lognormal relaxed clock and the uniform, exponential or lognormal prior distributions were employed in the analyses (And ujar et al., 2016;Chen et al., 2020;Toussaint & Short, 2017. To test the influence of different parameters on the divergence time estimation, six different parameter settings were chosen: (1) Yule model and exponential prior distributions; (2) birth death model and exponential prior distributions; (3) Yule model and lognormal prior distributions; (4) birth death model and lognormal prior distributions; (5) Yule model and uniform prior distributions; (6) birth death model and uniform prior distributions. We used two calibrations in the molecular dating for the divergences within Inocelliidae: (1) 132.1 Ma set as the minimum age for the divergence between Inocelliidae and Raphidiidae estimated by Vasilikopoulos et al. (2020). The ancestor of crown Raphidioptera was estimated to have lived at the beginning of the Cretaceous period (132. . Thus, under the exponential prior distribution we enforced the minimum age as 132.1 Ma and the soft maximum bound as 235 Ma (BEAUti settings: Mean = 28.4, Offset = 129.8). We also used uniform (132. 1-136.8) and lognormal (Mean = 20, stdev = 0.8, Offset = 128.3) prior distributions; (2) Fibla carpenteri Engel, 1998 from the Eocene Baltic amber (Engel, 1998) for the split between Fibla and Inocellia. Baltic amber is generally dated during the middle Eocene (Lutetian: 44.1 AE 1.1 Ma) with a minimum age of 43 Ma (Ritzkowski, 1997). As a result, under the exponential prior distribution we enforced the minimum age as 43 Ma and the soft maximum bound as 125 Ma (BEAUti settings: Mean = 24, Offset = 42.4). We also used a uniform prior distribution (43-45.2) and a lognormal prior distribution (Mean = 20, stdev = 0.8, Offset = 40).
We ran two independent MCMCs with 300 million generations each, sampling parameters and trees every 1000 generations. We used Tracer v.1.7  to check the convergence of MCMCs and whether the effective sample sizes exceeded 200. We resampled the tree files of both runs with a frequency of 10 000 and combined them in LogCombiner (BEAST package), of which the first 25% were removed as burn-in. We summarized the subsampled trees in a maximum clade credibility tree with median heights as node heights using TreeAnnotator (BEAST package). The median heights and 95% highest probability density (95% HPD) values were displayed in FigTree v. 1.4.4 (Rambaut, 2018) (Table S4).

Ancestral area reconstruction
We inferred ancestral areas for the major lineages of Inocelliidae using the BioGeoBEARS package (Matzke, 2013) implemented in RASP (Reconstruct Ancestral State in Phylogenies) v.4.0, build 20200103 (Yu et al., 2020). To account for phylogenetic uncertainty, 6002 post-burn-in trees yielded from BEAST for the above divergence time estimation were included with removal of outgroup taxa. Under the condition of species formation parameter (J) including "jump diffusion" or "founder effect", respectively, DIVALIKE + J was selected because it was estimated as the best model by Akaike Information Criterion (AICc)_wt value (Table S5)
In addition to the paraphyletic Inocellia, the I. fulvostigmata group was not recovered as monophyletic. The phylogenetic positions of the species belonging to this group were poorly resolved. Some species of the I. fulvostigmata group, such as Inocellia cheni Liu Asp€ ock, formed a monophylum based on the complex of male gonocoxites + gonapophyses + gonostyli 10 with a conspicuously long median projection (23 : 1), complex of gonocoxites + gonapophyses + gonostyli 11 with a pair of small denticulate processes on ventral margin (28 : 1) and endophallus with two pairs of tufts (36 : 1). The I. crassicornis group was corroborated to be monophyletic and supported by the male gonocoxite 9 feebly extending beyond ectoproct (12 : 0) and the complex of male gonocoxites + gonapophyses + gonostyli 11 distinctly produced posterodorsad (26 : 1). The species of Parainocellia were segregated into two groups. The three species from Europe formed a monophylum supported by the male gonocoxite 9 slightly extending beyond ectoproct (12 : 0) and the female sternum 7 without a posterior incision (41 : 0). Parainocellia burmana from Indochina was assigned as sister group to Inocellia (Epinocellia) weii sp.n. from Guangxi, southern China based on the presence of a conspicuously long median projection on the complex of male gonocoxites + gonapophyses + gonostyli 10 (23 : 1) and the endophallus dorsally with a scabrous region (34 : 1). The monophyly of Amurinocellia is supported by numerous apomorphic characters, two of which are unique: gonapophyses 9 longitudinally crossed (21 : 1) and robust lobe on tergum 8 (44 : 1).

Phylogenetic analysis based on mitogenomic data
For the molecular phylogenetic analysis, four datasets were assembled: (1) the PCG123 dataset (11 145 bp), (2) the PCG123rRNA dataset (13 280 bp), (3) the PCG12 dataset (7430 bp) and (4) the PCG_AA dataset (3715 bp). No nucleotide substitution saturation was detected in the gene partitions (Iss < Iss.c, p = 0.000) (Table S3). Twelve phylogenetic trees were obtained from the ML, BI and MP inferences of the four datasets (Fig. 3, Figs S7-S17). Monophyly of Inocelliidae was recovered in all results. In addition, Negha and Sininocellia were consistently recovered as the sister group to the remaining in all four datasets under ML and BI tree. The paraphyletic Inocellia was recovered in all results owing to Amurinocellia and Parainocellia being nested within the former genus. The major topological difference among the results regarded the position of Fibla. This genus was recovered as the sister group to the clade of the remaining inocelliids based on the PCG123 and PCG12 datasets under ML and BI analyses. In the results from the PCG123rRNA dataset under ML and BI analyses, Fibla was recovered as sister group to Negha + Sininocellia. The PCGAA dataset under ML, BI and MP analyses, however, assigned Fibla as the sister group to the clade comprising most inocelliids except Negha + Sininocellia, which also is supported by the MP analysis of the PCG123, PCG123rRNA and PCG12 datasets. Additionally, the I. crassicornis group was recovered as sister group to Amurinocellia + European Parainocellia based on the PCG_AA dataset under ML and BI analyses, but as sister group to Amurinocellia by the MP analysis of the same dataset as well as by the ML and BI analyses of the PCG123, PCG123rRNA and PCG12 datasets.

Phylogenetic analysis based on morphological plus mitogenomic data
The results from the combined dataset recovered a phylogenetic scheme largely consistent with that from the PCG_AA dataset (Figs S18 and S19). Based on the combined dataset, the MP analysis as performed through TNT yielded 30 MPTs (length = 667, consistency index = 0.55, retention index = 0.68). Indianoinocellia, which is not currently represented in the molecular data, was recovered as the sister group to Negha and these two genera were grouped with Sininocellia under the BI analysis, Sininocellia is consistently sister group to the remaining genera under the MP analysis. Within the paraphyletic Inocellia there are several monophyletic lineages consisting of species of the I. fulvostigmata group, such as the monophylum comprising I. cheni, I. Asp€ ock also formed a monophylum, which was recovered only in this combined dataset. However, the relationships among these lineages were not resolved.

Divergence time estimation
There were only slight differences in estimates between analyses using different models and parameters (see Fig. S20). According to the minimum value of stepping-stone sampling marginal likelihood estimate (SS MLE), we considered the divergence times estimated under the setting of birth death model and exponential prior distributions as the best result (see Table S8). The early divergence within Inocelliidae, referring to the split between Negha + Sininocellia and the remaining inocelliids, was dated during the Eocene at 49.

Ancestral area reconstruction
The ancestral area reconstruction for the major inocelliid clades based on BioGeoBEARS under the DIVALIKE + J model ( Fig. 4; Table S7) suggested a wide distribution of the ancestral inocelliid from western North America (A), Oriental Asia (B), western Palaearctic (C) and eastern Palaearctic (D). Within Inocelliidae, a vicariance event along with an extinction from Oriental Asia (B) was recovered referring to the divergence between Negha + Sininocellia and the remaining inocelliids. Thus, the ancestral range for Negha + Sininocellia was recovered to be in western North America (A), whereas that for the ancestor of the remaining inocelliids was found to be in western Palaearctic (C). Concerning the divergence between Negha and Sininocellia, a dispersal event from western North America (A) to Oriental Asia (B), followed by a vicariance event (A?AB?A|B), was found. Likewise, the ancestral range of Fibla + the Inocellia clade was estimated to be in western Palaearctic, while a prior dispersal from western Palaearctic (C) to Oriental Asia (B) and a subsequent vicariance (C?BC?B|C) was recovered for the split between Fibla and the Inocellia clade. Oriental Asia (B) was reconstructed as the ancestral range of the Inocellia clade. Notably, subsequent dispersals from Oriental Asia (B) to Europe (C), or from Oriental Asia (B) via eastern Palaearctic (D) to Europe (C), were recovered alongside the divergence between Amurinocellia and European Parainocellia as well as the speciation within the I. crassicornis group.

Intergeneric relationships within Inocelliidae
Our phylogenetic analysis based on both morphological and molecular data represents the first critical investigation on the intergeneric relationships of Inocelliidae. The results from the various analyses, whether morphological, mitochondrial genomic, or combined, contradict the previous hypothesis proposed by Asp€ ock et al. (1991), Fibla + (Negha + (Indianoinocellia + (Inocellia + Parainocellia))), in the following two aspects.
First, according to the results from the PCG_AA dataset under ML and BI analysis, and all four molecular datasets under MP analysis, Fibla was assigned as sister group to the Inocellia clade, whereas Negha and Sininocellia constituted the sister group of Fibla + the Inocellia clade. It should be further noted that Fibla was recovered either as the sister group of the remaining inocelliids based on the PCG123 and PCG12 datasets under ML and BI analysis, or as the sister group of Negha + Sininocellia based on the PCG123rRNA dataset under ML and BI analysis, but both results received relatively lower support values than the result from the PCG_AA dataset. In addition, previous studies based on larger genome-scale datasets recovered Negha + (Fibla + (Inocellia + Parainocellia)) (Vasilikopoulos et al., 2020; Winterton et al., 2018), which is consistent with our result based on the PCG_AA dataset. Morphologically, the serratulum (a peculiar sclerite of the male genitalia in many inocelliids) may provide supportive evidence to our result. Given the position of this sclerite, close and ventral to the gonocoxites + gonapophyses + gonostyli 10 (formerly addressed as parameres), it seems originally to be a part of the complex of segment 10, putatively interpreted as the gonapophyses 10. Because the paired gonapophyses 10 (herein interpreted as the paired serratulum, which is present only in Negha) is unquestionably a plesiomorphic state, the unpaired and scabrous serratulum, which is shared by Fibla and many Inocellia species, serves as a potential synapomorphy of Fibla + the Inocellia clade.
Second, the two western Nearctic genera, Indianoinocellia and Negha, clustered together as sisters based on the absence of MA stem in hind wing (6 : 1) and the long ectoproct (40 : 1). Nevertheless, the former character state also was assigned as an autapomorphy to the Inocellia clade. The co-occurrence of these two genera in western North America also may be indicative of their close relationship. An unexpected find in the combined analysis of morphological and molecular data was the assignment of the Asian endemic genus Sininocellia as the sister group of Indianoinocellia + Negha. However, we were unable to find supportive morphological arguments. Morphologically, Sininocellia shares some apomorphic states with the Inocellia clade, such as the proximal shift of the male gonostylus 9 on the inner side of gonocoxite 9 (16 : 1) and the presence of thick bristles on lateral portions of endophallus (33 : 1). However, these shared apomorphic characters can equally be considered as convergent derivations in each of these two lineages.
The following three items are noteworthy concerning the present results that are inconsistent to the traditional morphology-based classification. First, a paraphyletic Inocellia was consistently recovered with the inclusion of Amurinocellia and Parainocellia grouped with a partial set of Inocellia species in all results from the various datasets. Previously, the hook-like male gonapophysis 9 (pseudostyli) were considered as a key character to verify the generic status of Amurinocellia and Parainocellia (Asp€ ock et al., 1991;Liu et al., 2009a). However, this character state is herein interpreted to be plesiomorphic by comparison with the raphidiid outgroups as well as the inocelliids with foliate male gonapophyses 9. Thus, despite the inclusion of this character state in our morphology-based analysis, it may not suggest the independent generic status of Amurinocellia and Parainocellia. It might be retained in some hetergeneous inocelliids, whereas the foliate male gonapophyses 9 might have become derived in the other genera and species. Furthermore, P. burmana and a related new species from Oriental Asia, despite having the hook-like male gonapophyses 9 (20 : 0), greatly differ from Amurinocellia and the European Parainocellia in other parts of male and female genitalia, and were not placed in the latter two groups. Therefore, even if the hook-like male gonapophyses 9 are indeed apomorphic, they might have evolved convergently in heterogeneous inocelliids. The absence of a serratulum in Amurinocellia, Parainocellia and the Inocellia crassicornis group also may support the monophylum including these three lineages. Given the above interpretation of the serratulum as the male gonapophyses 10, the complete loss of this sclerite may be the most derived condition. However, as a consequence of the unknown condition of the homologous sclerite of serratulum in Raphidiidae and the contradictory previous interpretation of the presence of serratulum as the apomorphic state , we did not code this character in our analysis, pending further examination of more snakefly taxa.
Second, the monophyly of the Inocellia fulvostigmata group, which is diagnosed by a single character--the male gonocoxite 9 much longer than wide (12 : 2) (Asp€ ock et al., 1991;Liu et al., 2010b)--is challenged by our results. The prolonged male gonocoxite 9 is also present in some species of Fibla and Negha, which represent relatively earlier divergent lineages in the family. Thus, the prolonged male gonocoxite 9 may be plesiomorphic. The phylogenetic positions of the species previously placed in the I. fulvostigmata group were poorly resolved using the present datasets and should await further investigation. Third, F. pasiphae as the single representative of the subgenus Reisserella H. Asp€ ock & U. Asp€ ock surprisingly appeared as nested within the subgenus Fibla in the present result. This species differs markedly from the other three species of Fibla as a consequence of striking differences of morphological characters, such as the antenna considerably longer than the forewing and the increasing number of the forewing discoidal cells (cells between the branches of MP vein). The subgeneric assignment of this species awaits further verification.
According to the present phylogenetic inference on the intergeneric relationships within Inocelliidae, the monophyly of each of the following genera Fibla, Indianoinocellia, Negha and Sininocellia was corroborated. Within the Inocellia clade, the monophyly of Amurinocellia and of the European Parainocellia was verified, but P. burmana together with a presently described new species represent a distant but related lineage. Because these three lineages all possess a distinct hook-like male gonapophyses 9, they are treated here as subgenera of Inocellia (i.e. Amurinocellia stat.r. and Parainocellia stat.r. et emend.n.), and particularly, P. burmana and the new species are placed in a new subgenus, namely Epinocellia subgen.n. The monophyletic I. crassicornis species group is now the sole constituent of the subgenus Inocellia stat.n. (see Appendix 2 and File S1).

Spatial and temporal diversification of Inocelliidae
Extant snakefly families are probably relicts that survived the Cretaceous-Palaeogene (K-Pg) mass extinction owing to their adaption to a cold climate (Asp€ ock, 1998), with the divergence of Inocelliidae and Raphidiidae taking place during the Early Cretaceous or as early as the middle Jurassic based on evidence from molecular dating (Vasilikopoulos et al., 2020;Wang et al., 2017;Winterton et al., 2018). This implies that during the Mesozoic certain stem-group inocelliids or extinct lineages had a closer relationship to Inocelliidae than to Raphidiidae. Such a hypothesis was supported by a recent finding of diverse species of the extinct family Baissopteridae, which possess male genital characters similar to Inocelliidae (Lu et al., 2020). Nevertheless, no definite inocelliid has been found from the Mesozoic, and the relationships between the diverse fossils and extant snakeflies remain elusive (Lu et al., 2020;Lyu et al., 2020). Therefore, at present, the cladogenesis and its spatio-temporal mode can only be traced back to the Eocene at a time when inocelliids, according to fossil evidence, explicitly occurred (Asp€ ock & Asp€ ock, 2004; Makarkin et al., 2019).
Based on the present divergence time estimation, diversification of the extant inocelliid lineages initiated during the early Eocene. This estimation corresponds to that for the divergence of Negha and Fibla + Inocellia based on transcriptome data (Vasilikopoulos et al., 2020). Moreover, the estimate is compatible with the recent finding of the hitherto oldest fossil Inocelliidae, which is from the early Eocene (Ypresian) in Okanagan Highlands, Canada (Makarkin et al., 2019). The fossil record of Eocene Inocelliidae from Europe and western North America (Carpenter, 1958;Engel, 1998Engel, , 2002Asp€ ock & Asp€ ock, 2004;Makarkin & Archibald, 2014;Makarkin et al., 2019) corroborates the occurrence of ancestral inocelliids, as herein inferred, in these two regions. A notable scenario refers to the early Eocene extinction of ancestral inocelliids in East Asia where currently a particularly rich and diverse set of modern species of the family is present. This event may have been caused by a broad arid/steppe band extending from East China to Central Asia which formed during the early Eocene Zhang et al., 2012) and made the habitat in this region barely suitable for the arboreal lifemode of inocelliids.
Two dispersal events to Oriental Asia were recovered that occurred in the mid to late Eocene, one from western North America and the other from the western Palaearctic. The colonization of ancestors of Sininocellia from western North America to East Asia was probably routed over the Bering land bridge or was at least facilitated by the land bridge, which is considered an important channel for Holarctic biotic exchange (Gladenkov et al., 2002;Jiang, Klaus, et al., 2019;Sanmart ın et al., 2001;Wen et al., 2016). However, it would be problematic to postulate a dispersal of the Inocellia clade from Europe to East Asia during the middle Eocene. The Turgai Sea, which persisted from the middle Jurassic until the Oligocene (Sanmart ın et al., 2001), probably acted as a major barrier for the faunal exchange between the western and eastern sides of the Palaearctic. As such, it is unlikely that during the Eocene, the Inocellia clade could have originated in Europe and subsequently colonized East Asia. But, notably, the Eocene fossil record of Fibla ['Fibla' exusta (Cockrell & Custer, 1925)] from western North America (Makarkin & Archibald, 2014) may provide an alternative explanation, whereby the ancestor of the (i) female genital segments, lateral view; (j) bursa copulatrix, lateral view. Abbreviations: ab, atrium bursae; e, ectoproct; ep, endophallus; gp9, gonapophysis 9; gx9, gonocoxite 9; gx10, complex of fused gonocoxites, gonapophyses and gonostyli 10; gx11, fused gonocoxites 11; gr, glandula receptaculi; gst9, gonostylus 9; hi, hypandrium internum; r, receptaculum seminis; sb, sacculus bursae; sg, fused gonocoxites 8; S, sternum; T, tergum. Scale bars = 0.5 mm. Fig. 7. Illustrations of male genital sclerites of inocelliids. Arrows indicate character state used in the phylogenetic analysis of the morphologyonly dataset as indicated in the text. Colour scheme is presented at the bottom. Genital segments of male are symbolized by colours: light blue stands for gonocoxite 9, medium blue for gonostylus 9, dark blue for gonapophyses 9, red for complex of fused gonocoxites, gonapophyses, gonostyli 10 (fused parameres), light green for fused gonocoxites 11 (gonarcus), orange for endophallus and purple for serratulum. Scale bars = 0.5 mm.
Inocellia clade would have dispersed from western North America, where its sister group (Fibla) occurred, to East Asia via the Bering land bridge, taking the same route contemporaneously to the dispersal leading to the origin of Sininocellia in western North America. Nevertheless, the generic affiliation of the single Fibla fossil from North America is provisional because the specimen is known only by a piece of hind wing (Makarkin & Archibald, 2014). Some of its characters, especially the presence of the long MA stem, also are shared by Sininocellia and the Baltic amber Succinofibla Asp€ ock & Asp€ ock, 2004. Therefore, the present view on the origin of the Inocellia clade requires further investigation. Fig. 8. Illustrations of male genital sclerites of snakeflies. Arrows indicate character state used in the phylogenetic analysis of the morphologyonly dataset as indicated in the text. Colour scheme is presented at the bottom. Genital segments of male are symbolized by colours: light blue, gonocoxite 9; medium blue, gonostylus 9; dark blue, gonapophyses 9; red, complex of fused gonocoxites, gonapophyses and gonostyli 10 (fused parameres); light green, fused gonocoxites 11 (gonarcus); orange, endophallus; and purple, serratulum. Scale bars = 0.5 mm.
After the Eocene divergence of recent inocelliid genera, speciation within the genera was dated to be mainly during the Oligocene and Miocene (Fig. 4). Taking the Inocellia clade (the bulk of the inocelliid diversifications) as an example, its speciation was initiated during the late Oligocene and early . The latest speciation, represented by the split between I. (A.) sinica and I. (A.) sinensis, took place during the late . The development of the East Asian monsoon and the uplift of the Tibetan Plateau within the Neogene brought momentous changes in climate, landform and biota (Cao et al., 2012;Ding & Liao, 2019;Li et al., 2016;Quan et al., 2021). As the predominant inocelliid group from East Asia, the speciation and distribution of the Inocellia clade is likely to have been influenced by the dramatic topographic and climatic effects during the Miocene. Based on the present distribution, most species of Inocellia occur in a vast area from south of the Qinling Mountains to northern Indochina. This area belongs to the Oriental region but its northern part is Fig. 9. Illustrations of female genital sclerites of snakeflies. Arrows indicate character state used in the phylogenetic analysis of the morphologyonly dataset as indicated in the text. Colour scheme is presented at the bottom. Genital segments of female are symbolized by colours: light yellow, fused gonocoxites 8; dark yellow, gonapophyses 8; light blue, gonocoxite 9; and grey, atrium bursae. Scale bars = 1.0 mm. considered a transitional zone between the Palaearctic and Oriental regions (Zhang et al., 2015;Schmitt, 2020). The summer monsoon is a defining factor in this area, which harbours a rich diversity of Inocellia, suggesting that the ecological preference of Inocellia is for humid habitats and not arid. Since the early Miocene, the arid belt, which once extended to southeast China during the Palaeogene, retreated northwestward, which led to humidification in southwest and southeast China (Guo et al., 2008). This Miocene climate change probably triggered an increase of suitable habitats for Inocellia in East Asia. Synchronously, the Tibetan uplift resulting from the collision between India and Eurasia was accelerated, and the large-scale orogeny associated with the growing Himalayas (Spicer et al., 2021;Su et al., 2021;Xu et al., 2021), such as the formation of the Hengduan Mountains, might have played an important role in providing geographical isolation for allopatric speciation. The high species diversity and endemism of Inocellia along the Himalayas as well as in the Hengduan Mountains (Liu, Asp€ ock, Zhang, et al., 2012, Liu et al., 2018Shen et al., 2019) stresses the importance of this topographic factor as a driving force for speciation.
Aside from this broad impact on the diversification of Inocellia from East Asia, the long-distance dispersal from East Asia to the western Palaearctic of some inocelliids is another notable issue recovered in this study. First, the distribution of the subgenus Parainocellia is restricted to regions north of the Mediterranean, whereas the parent genus, Inocellia, has an East Asian-Tethyan disjunct distribution, in common with a wide range of animals and plants (Jiang, Hipp, et al., 2019). The middle Miocene dispersal of ancestral Parainocellia from East Asia to Mediterranean Europe may have occurred via the corridor in Tibet-Himalaya and landmasses south of the Paratethys (represented mainly by the Alpine-Carpathian-Dinarides mountain chain). This dispersal route is considered important for biotic exchange within Eurasia during the Oligocene into the early Neogene based on an empirical study on the biogeographical history of European holly oaks, which share the East Asian-Tethyan disjunct distribution (Jiang, Hipp, et al., 2019). Second, the I. crassicornis group, which has many species in eastern Palaearctic (e.g. Inocellia biprocessus Liu, H. Asp€ ock, Yang & U. Asp€ ock, Inocellia fujiana Yang and Inocellia japonica Okamoto) might have diversified during the expansion of its distribution northward from Oriental China into eastern Palaearctic since the late Miocene. The wide distribution of I. crassicornis in both eastern and western Palaearctic is likely to have been shaped by a Pliocene dispersal originating from the eastern Palaearctic. This migration differs from the putative dispersal route for the ancestral Parainocellia and probably followed a pathway along the coniferous belt of the northern Palaearctic (Hines, 2008).
Our study primarily documented the cladogenesis of Inocelliidae during the timespan from the early Eocene to Pliocene. However, the recent distribution pattern of inocelliids is apparently not only a consequence of various Tertiary factors but also the Quaternary glaciations. Asp€ ock et al. (1991) and some of their later studies (Asp€ ock, 1998;Asp€ ock et al., 2011) attempted to analyse the biogeography of Raphidioptera by considering species to be faunal elements of a certain refugial centre following the concept of de Lattin (1967). This concept is plausible in particular for the glacial periods (and preferably the last glacial period) and concerns the Palaearctic region (Asp€ ock et al., 2011), but also the Nearctic region. However, the diversification center of Inocelliidae appears to be located in Oriental China, although the Oriental species are distributed at higher elevations and were once considered as the elements of a transition zone between the Palaearctic and Oriental regions . Radiation of Inocellia at the species level in Oriental Asia during the Miocene, as inferred in this study, clearly demonstrates that at least a number of East Asian inocelliids represent the Oriental elements. Nonetheless, because the sampling used in our molecular phylogenetic and biogeographical analyses did not include all species of the family, certain potential scenarios shaping the modern distribution of Inocelliidae in relation to the Quaternary glaciations cannot be tested precisely at present.

Conclusions
The present study is the first to elucidate the evolutionary history of the extant snakefly family Inocelliidae through phylogenetic analysis, molecular dating and ancestral range reconstruction based on a combined dataset from both morphological characters and mitogenomes. Under the phylogenetic framework herein inferred, the paraphyly of the species-rich genus Inocellia, as well as its relationship with Amurinocellia and Parainocellia, was recovered, which resulted in a new classification of Inocelliidae. The diversification within the family is traced back at least to the early Eocene. Both the Eocene and the Miocene are crucial phases for the intergeneric and interspecific divergences along the timeline of the inocelliid cladogenesis. The diverse inocelliid fauna from East Asia might have a western North American origin via dispersal across Beringia during the early Tertiary. Thereafter, Oriental China and adjacent areas to the south hosted the radiation of inocelliids with an increase in geographical isolation and suitable coniferous habitats resulting from the Tibetan uplift and development of East Asian summer monsoon. Moreover, from this inocelliid diversification centre, long-distance dispersals via the Tibet-Himalayan corridor or eastern Palaearctic to western Palaearctic would explain the European distribution of some younger lineages in the Inocellia clade. Taking the present work as a starting point, a totalevidence phylogenetic analysis combining all extant and fossil taxa as well as larger genomic-scale data would be essential for further understanding the evolutionary history of Inocelliidae, and a phylogeographic study on certain widespread species together with new biological and ecological data from additional species will be helpful to investigate the Quaternary glacial impact on inocelliid evolution.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.     Liu et al. (2010b); (j) modified from Liu and Hajong (2015); (m) modified from Liu, Asp€ ock, Zhang, et al. (2012); (n) modified from Liu et al. (2009b).   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. Fig. S19. Maximum parsimony tree based on the combined analysis of morphological and molecular data. 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. 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. Table S5. Likelihood scores and model comparison for ancestral area reconstruction using BioGeoBEARS. Table S6. Major biogeographical events of Inocelliidae estimated with BioGeoBEARS. Table S7. Summary of ancestral area reconstructions based on DIVALIKE + J model. Table S8. Beast median divergence times for major nodes in Fig. 4, with 95% credibility intervals and marginal likelihoods.
File S1. Checklist of extant species of Inocelliidae. File S2. Specimens examined for the present study. File S3. Nexus file of MP tree based on morphological data.
File S4. TNT file of MP tree based on morphological data.
File S5. Nexus file of BI tree based on dataset PCG123 (nucleotide data of 13 PCGs).
File S6. Phylip file of ML tree based on dataset PCG123 (nucleotide data of 13 PCGs).
File S7. Fasta file of MP tree based on dataset PCG123 (nucleotide data of 13 PCGs).
File S8. Nexus file of BI tree based on dataset PCG123rRNA (nucleotide data of the 13 PCGs and rRNA genes).
File S9. Phylip file of ML tree based on dataset PCG123rRNA (nucleotide data of the 13 PCGs and rRNA genes).
File S10. Fasta file of MP tree based on dataset PCG123rRNA (nucleotide data of the 13 PCGs and rRNA genes).
File S11. Nexus file of BI tree based on dataset PCG12 (nucleotide data of the 13 PCGs without the third position of PCGs).
File S12. Phylip file of ML tree based on dataset PCG12 (nucleotide data of the 13 PCGs without the third position of PCGs).
File S13. Fasta file of MP tree based on dataset PCG12 (nucleotide data of the 13 PCGs without the third position of PCGs).
File S14. Nexus file of BI tree based on dataset PCG_AA (amino acid sequences of the 13 PCGs).
File S15. Phylip file of ML tree based on dataset PCG_AA (amino acid sequences of the 13 PCGs).
File S16. Fasta file of MP tree based on dataset PCG_AA (amino acid sequences of the 13 PCGs).
File S17. Nexus file containing the molecular and morphological alignments.
File S18. TNT file containing the molecular and morphological alignments.