Evolution, systematics and historical biogeography of Palparini and Palparidiini antlions (Neuroptera: Myrmeleontidae): Old origin and in situ diversification in Southern Africa

Palparine and palparidiine antlions constitute an emblematic clade of large and occasionally colourful insects that are only distributed in the western portion of the Eastern hemisphere, with about half of the known species diversity occurring exclusively in Southern Africa. Little is known about their evolutionary history, and the boundaries and relationships of most genera are still unresolved. In this study, we analyse a molecular dataset consisting of seven loci (five mitochondrial and two nuclear genes) for 144 antlion species and provide the first phylogenetic hypothesis for a representative sampling of Palparini and Palparidiini (62 Palparini species, representing 15 of the 17 known genera, and all three known Palparidiini species). In addition, we reconstruct their timing of diversification and historical biogeography. The resulting tree indicates that several extant palparine genera are polyphyletic or paraphyletic and provides interesting leads that ought to be helpful for future taxonomic revisions; it also enables us to re‐evaluate the taxonomic utility and relevancy of a number of morphological characters that were previously used to define some genera. Molecular dating analyses indicate that the most recent common ancestor of both groups originated about 92 million years ago (Ma) in the Late Cretaceous. Finally, the results of historical biogeography analyses provide strong support for an origin in Southern Africa, which further acted as both a cradle of diversification and a springboard for successive waves of northern dispersals.

Mervyn W. Mansell and Bruno Michel jointly supervised this work. Palparidiini are merged under Palparini (along with Maulini sensu Stange, 2004). This merging is supported by the fact that Palparidiini and Palparini are recovered as sister groups in all molecular and morphological investigations where they have been sampled Badano, Aspöck, Aspöck, & Haring, 2017;Machado et al., 2019). By contrast, Jones (2019) elevated the subfamily Palparinae to full family status. However, this decision was not substantiated, with only the family name Palparidae being mentioned in isolation. In Badano, Aspöck, Aspöck and Haring (2017), a new classification system is proposed for Palparinae, with two tribes (Palparini and Stilbopterygini) consisting of four (Dimarina, Echthromyrmecina, Palparidiina and Palparina) and two subtribes (Pseudimarina and Stilbopterygina), respectively; the corresponding changes are based on the analysis of a 29 species morphological dataset (with 51 characters) and the analysis of a 10 species molecular dataset (with three markers). All this taxonomic instability is further exacerbated by the fact that the recent phylogenomic studies using anchored hybrid enrichment data (Machado et al., 2019;Winterton et al., 2018) or transcriptomes (Vasilikopoulos et al., 2020) have yielded conflicting results regarding the monophyly of Myrmeleontidae, both with maximum statistical support. The latter is not unexpected given that phylogenomic analyses are sensitive to sampling and not impervious to systematic errors (see the review of Young & Gillung, 2020); the inference of conflicting topologies with an equally high support is also not unusual in phylogenomics analyses, where support for incorrect relationships is sometimes inflated by the sheer amount of data analysed (Brown & Thomson, 2017;Jeffroy, Brinkmann, Delsuc, & Philippe, 2006). Finally, it is worth underlining that the recent study of Badano et al. (2021), in which a part of the genomic dataset of Winterton et al. (2018) was combined with a morphological dataset, recovered the monophyly of Myrmeleontidae. Therefore, and until a firmer consensus based on additional investigations is available, we decided to exercise caution and adhere to the traditional classification system of Stange (2004).
Palparini antlions include some of the largest and most spectacular species in the order Neuroptera, and indeed, of all insects. They are characterized by their unusually large size and striking wings (wingspans ranging from 45 to 170 mm), and occasionally conspicuous colouration, which render them unmistakable (Figure 1). Palparini consists of 136 species assembled in 17 genera (Table S1) that are mostly distributed throughout the Afrotropical region and extending across the southern Palearctic to the Oriental region, but are absent from the Australasian region and the Western hemisphere (Mansell, 1990;Stange, 2004). Of the 98 species exclusively found in the Afrotropics, 44 are found in 11 genera occurring in Southern Africa (Botswana, Namibia, South Africa, southern Mozambique and Zimbabwe), with 35 species being endemic to the sub-region. Most of the endemicity, 28 species in nine genera, is centred in the Western and Northern Cape Provinces of South Africa and Namibia. Palparini are mostly found in dry or wet savannas, but few species can also be found in semi-arid deserts or woodlands (Mansell & Erasmus, 2002;Michel, 1999;Prost, 1995). Adults generally fly erratically and over short distances.
Palparini are generally nocturnal predators that are attracted by lights; however, this is not always the case and some species are diurnal, such as those in the genera Pamares Mansell and Pamexis Hagen (Mansell, 1990(Mansell, , 1992b. During daytime, they generally rest with their wings folded; their intricately patterned wings serve as an effective camouflage (e.g., see some of the pictures in Figure 1), usually matching the predominant vegetation of the environment in which they live (Mansell, 1999). The Palparini also have the distinction of having only two species of antlions that mimic lichens (Mansell & Ball, 2016; see the picture of Pamexis namaqua in Figure 1). Palparini larvae (Figure 2) are psammophilous and the protection afforded by their deep sand habitats is a key factor that potentially accounts for their large size, with the larvae of some species reaching a length of 35 mm (Mansell, 1999). Because of their size, Palparini larvae are able to subdue a large range of prey (including grasshoppers) and are an important component of the predator guild in the areas they inhabit (Mansell, 1999). Most live in deep sand and migrate to the surface in the late afternoon and evening when sand temperatures cool. All larvae whose ecology is known do not construct pitfall traps and just lay below the surface with their head and spread jaws exposed, waiting to be alerted by vibrations from approaching prey. The six eye facets set on a prominent tubercle play a role in directing the jaw strike towards their prey. Once secured by the tips and intermeshing teeth on the sickle-shaped mandibles, the larva moves backwards to subdue the prey and protect itself from retaliatory injury. Historically, most species of Palparini have been described in the genus Palpares Rambur, the type genus of the subfamily Palparinae and tribe Palparini. It has long been postulated that this genus is polyphyletic and several attempts have been made to clarify its taxonomy. Noticeably, Mansell (1992a) arranged this genus (and other palparine genera) into four divisions and a number of species groups, based on morphology.
Palparidiini consist of only three species belonging to the genus Palparidius Péringuey (Table S1). Adults are quite large in size (wingspans up to 90 mm); morphologically they share several characters with Dimarini and Palparini (Mansell, 1999;Stange, 2004) but differ by their highly modified male ectoproct. The larvae of two species (Palparidius capicola Péringuey and Palparidius concinnus Péringuey) are known (Stange, 2004;Tippett, 2022), and present characters that are similar to both Dimarini and Palparini (Stange, 2004). That said, the hypothesis that Palparidiini are more closely related to Palparini is clearly supported by the phylogenomic study of Machado et al. (2019) and by the cladistic study on larval characters of Badano et al. (2018).
Extant members of the tribe Palparidiini are endemic to arid and semiarid environments of Southern Africa. Palparidius larvae are psammophilous predators (Mansell, 1999) that likely live in deep sand; although little is known about their ecology, it can be postulated that both adults and larvae behave similarly to Palparini.
At the moment, no molecular studies have investigated in detail the evolutionary history of Palparini and Palparidiini at the generic and species levels. For instance, the phylogenetic study with the highest number of Palparini species to date (Michel et al., 2017) only encompasses 15 palparine species representing six genera. The current study is consequently a contribution to provide more data for molecular analyses investigating the evolutionary history of the clade consisting of Palparini and Palparidiini, for which a wide coverage of genera and species, especially from southern and western Africa and other disparate geographical regions, is now available. This opportunity arises from the extensive coverage of the group as a result of collections over many years in widely disparate areas, particularly in the Afrotropics. In addition, reconstructing the evolutionary history of this group offers the opportunity of exploring its biogeographic history, which presents a particular interest given the paucity of research focusing on the origin and diversification of Afrotropical insect groups.
Here, we provide the first comprehensive phylogenetic framework for Palparidiini and Palparini antlions. We further conduct molecular dating and historical biogeography analyses to deepen our understanding of the evolutionary history of the two clades, especially in relation to the diversification of Afrotropical lineages, by looking at their age, origin and colonization dynamic.

Taxon sampling, DNA extraction and sequencing
For this study, 121 specimens representing 55 Palparini species (with representatives of 14 of the 17 known genera) and all three Palparidiini species were sampled and sequenced successively, except two specimens (see Table S2). The material generally corresponds to specimens collected by the two senior authors, with most specimens col-  (Akoudjin & Michel, 2011;Mansell, 1990Mansell, , 1992aMansell, , 1992bMansell, , 1996Mansell, , 2004Mansell, 2018;Mansell & Ball, 2016). For each specimen, DNA was extracted from one hind leg using BioBasic EZ-10 96 Well Plate  Table S3). Initially, we used Sanger sequencing (see Michel et al., 2017 for details), but because DNA was often degraded, the majority of PCR failed (especially when targeting the 28S gene). Therefore, we relied on highthroughput sequencing (amplicon sequencing) to amplify gene fragments for the cox1, rrnL and rrnS genes. Amplicon libraries were constructed for these three genes following Galan et al. (2017); for the cox1, two overlapping fragments were also targeted following Shokralla et al. (2015). Compared with the settings of Galan et al. (2017), we made the following changes to lower the proportion of chimeric fragments: for the first PCR step, we changed the number of cycles to 40 and the extension period for the second PCR step, the extension duration was set to 120 s. The final library was paired-end sequenced on an Illumina MiSeq flowcell using a MiSeq Reagent Kit v2 (500 cycles) at the AGAP laboratory (Montpellier, France). Illumina reads were processed using the FROGS pipeline (http://frogs. toulouse.inra.fr/; Escudié et al., 2018) on the Genotoul Galaxy server using demultiplexing, pre-processing, clustering and chimera removal tools. The remaining contaminants were further detected using the BLAST tool (available at: https://blast.ncbi.nlm.nih.gov/Blast.cgi) and removed manually. The two overlapping fragments of cox1 were merged using Mesquite v3.70 (Maddison & Maddison, 2021). All mitochondrial and nuclear sequences were aligned using MAFFT v7 (Katoh & Standley, 2013) with default option settings and a gap opening penalty of 5.0, and further manually corrected using Mesquite. For all protein-coding genes (cob and cox1), coding frames and stop codons were also checked with Mesquite to detect potential pseudogenes.

Molecular datasets
Newly generated sequences were combined with extant data available on GenBank, including data (for 50 species) generated by our research group in the study of Michel et al. (2017) Table S2).

Phylogenetic analyses
Phylogenetic studies were conducted under maximum likelihood (ML) using IQ-TREE v2.1.3 (Minh et al., 2020). In a preliminary analysis, the specimen-level dataset was used to evaluate potential species paraphyly. We further analysed the species-level dataset, to generate a phylogenetic tree to be used as a reference tree for the dating and historical biogeography analyses. For both datasets, the same phylogenetic procedures were carried out. The concatenated dataset was partitioned into 13 partitions, with one partition established for each non-coding gene fragment (rrnL, rrnS, 18S and 28S) and three partitions (one for each codon position) defined for each coding gene fragment (cob, cox1 and cox3). The Bayesian information criterion (BIC) implemented in IQ-TREE was used to choose the best-fit substitution models and partition schemes (Table S4). Best-scoring trees were obtained using heuristic searches implementing 500 random-addition replicates, with the following settings: random-starting tree, thorough hill-climbing nearest neighbour interchange (NNI) search (Àallnni option), a perturbation strength set to 0.2 (Àpers 0.2 option), partition-resampling strategy (ÀÀsampling GENE option) and best partition scheme allowing the merging of partitions (Àm MFP + MERGE option).

Dating analyses
Divergence times were estimated using Bayesian relaxed clocks as implemented in BEAST v1.10.4 . Dating analyses relied on three vetted fossil constraints allowing to specify minimum age constraints (see Appendix S2 for additional details). A conservative maximum age of 201.3 Ma (corresponding to the Jurassic/Triassic boundary) was chosen for these three constraints since it is significantly older than the appearance of any myrmeleontoid in the fossil record Michel et al., 2017). These three constraints were enforced using uniform distributions; two distinct uncorrelated lognormal clocks were used for the mitochondrial and the nuclear genes, and the tree model was set to a birth-death speciation process (Gernhard, 2008). A fixed topology corresponding to the best-scoring tree from the ML analyses of the species-level dataset was used to

Historical Biogeography analyses
Historical biogeography analyses were carried out using parametric methods (Ree & Sanmartín, 2009) with RASP v4.2 (Yu, Blair, & He, 2020). We used the BioGeoBears (Matzke, 2014) implementation, which can run and compare both dispersal-vicariance analysis (DIVA; Ronquist, 1997), Dispersal-Extinction-Cladogenesis (DEC; Ree & Smith, 2008) and patterns, a maximum number of two areas was allowed (disjunct areas were also excluded). A time-stratified model was then defined ('reference model'), implementing a matrix of scaling factors for dispersal rates (DR) between areas (from 0 to 1) to account for the respective positions of the geographic areas through time while also accounting for potential geographical barriers to antlion dispersal (e.g., deserts, forests, high reliefs). A DR of 1.0 was set in the absence of barriers between two adjacent areas; a DR of 0.7 was set in the presence of a small or temporary barrier between two adjacent areas (e.g., cycles of fragmentations/reconnections of rainforests, high reliefs); a DR of 0.5 was set in the presence of a large barrier between two adjacent areas (e.g., deep sea, desert, rainforest); a DR of 0.3 was set to account for potential dispersal between two non-adjacent areas separated by one area; a DR of 0.1 was specified for long-distance dispersals, when two areas are separated by more than one area. Dispersal between Madagascar and other areas other than Eastern Africa was also considered as long-distance dispersal.
To better take into account changes through time three time slices were used. The first time slice ('TSI') runs from 120 to 56 Ma.
At that time India was initially connected to Madagascar (until 90 Ma) and progressively drifted towards Eurasia (Blakey, 2008); to account for that, a DR of 0.5 was specified between areas A and G, and between areas A and B. Several barriers between West and East Africa have been hypothesized to have existed during the Late Cretaceous and the Paleocene, first in the form of a large intercontinental desert during the Late Cretaceous (Carvalho, de Gasparini, Salgado, de Vasconcellos, & da Silva Marinho, 2010;DeConto, Hay, Thompson, & Bergengren, 1999), followed by the rise of a large pan-African rainforest during the Paleocene (Morley, 2000(Morley, , 2007. To account for these, a DR of 0.5 was set between areas C, D, E and F. A second time slice ('TSII') runs from 56 to 23 Ma. This timeframe corresponds to the Eocene and Oligocene. At that time India was still connected to Madagascar via a chain of islands and was no longer separated from Eurasia by an ocean, colliding with it ca. 45 Ma (Blakey, 2008); thus, a DR of 0.5 was set between areas A and G, and a DR of 0.7 was set between A and B. This period is also marked by the Eocene-Oligocene transition (EOT, 33 Ma;Zachos, Pagani, Sloan, Thomas, & Billups, 2001;Zachos, Dickens, & Zeebe, 2008) characterized by a global cooling leading to the first fragmentation of the pan-African rainforest (Couvreur et al., 2021;Morley, 2000). Therefore, a DR of 0.7 was set between areas C, D, E and F. This period is also marked by intense volcanic activity in Ethiopia during the Oligocene (Coulié et al., 2003;Couvreur et al., 2021), which likely acted as a barrier to dispersal; thus a DR of 0.7 was set between areas B and E. A third time slice ('TSIII') runs from 23 Ma to the present. At that time India was no longer connected to Madagascar, and to account for this, the DR between A and G was set to 0.1. Until the beginning of the Miocene ca. 23 Ma, cycles of rainforest fragmentation and reconnection brought on by changes from humid and hot environments to dry and cold conditions occurred (Couvreur et al., 2021;Morley, 2000). Also, this time slice corresponds to a period of major uplifts in the Afrotropics, and particularly in Eastern Africa (Guillocheau et al., 2018;Sepulchre et al., 2006). To take that into account, the DR between areas C, D, E and F was set to 0.5.
The apparition of three deserts also probably acted as barriers: Namib desert 16 Ma (DR of 0.5 between areas D and F), Saharan desert 7 Ma and Arabian deserts 3.5 Ma (DR of 0.7 between areas B and C, and DR of 0.5 between B and E). Finally, for comparison purposes, we also defined a model without any constraints (DR of 1.0 between all areas) nor time slices ('null model').
Analyses were carried out for each of these models both with and without the founder-event speciation parameter (+j) sensu Matzke (2014). Models were selected based on the Akaike information criterion corrected for sample size (AICc_wt), following Yu et al. (2020). As a guide tree, we used the dated phylogeny estimated with BEAST; this tree was modified in Mesquite ('prune clade' tool) by removing all species not belonging to the Palparini and Palparidiini tribes.

Phylogenetic analyses
The results of the ML analysis of the specimen-level dataset are presented in Figure S1. In the corresponding tree (likelihood score of À38310.305), all species for which multiple specimens were included are recovered monophyletic. Overall, the inferred relationships for Palparini and Palparidiini are similar to those inferred through the analysis of the species-level dataset (see below). The analyses carried out to test the impact of missing data yield branch support values that are generally lower and also result in the misplacement of several lineages (see Appendix S1 for details).
The best-scoring tree from ML analysis of the species-level dataset has a likelihood score of À67640.115 (see

Historical Biogeography analyses
Likelihood ratio tests indicate that +j models provide the best fit across all analyses. The reference model fitted the data better than the null model (lnL = À116.7 vs lnL = À122; Table S5). According to the AICc_w model selection, the DEC+j and DIVALIKE+j models are

Palparellus flavofasciatus (7)
(1) F I G U R E 3 Maximum Likelihood tree resulting from the analysis of the species-level dataset. Node support values are indicated using black circles when supported by three metrics (SH-aLRT ≥80%, uBV ≥95% and TBE ≥70%), by grey circles when supported by two metrics, by white circles when supported by one metric; the absence of a circle indicates that the node is not statistically well-supported (see Figure S2 for exact support values). Information on taxonomic ranks sensu Stange (2004)

Phylogenetic relationships and systematics
First, we want to stress that this study was not designed to investigate the relative placements of owlfly and antlion lineages, hence the fact that we recovered a monophyletic antlion family does not constitute a finding of particular interest; it was also anticipated given that Stilbopteryginae based on all of these results. In any case, our results are in line with those of Badano, Aspöck, Aspöck, and Haring (2017).
To answer this question definitively, a phylogenomic investigation based on a large sample of Palparinae and members of several Ascalaphidae and Myrmeleontidae lineages will likely be required. The paraphyly of Myrmeleontinae sensu Stange (2004) is due to the placement of members of the tribe Acanthaclisini. This result supports the views of Oswald and Penny (1991) and New (1985New ( , 2003, who do not recognize the acanthaclisines as a myrmeleontine tribe, and instead consider them as a subfamily of their own. Furthermore, this lineage is well-defined by the morphological characters of the larvae and adults and clearly constitutes a homogeneous clade (Hölzel, 1972;Insom & Carfì, 1992;Markl, 1954;Stange & Miller, 1985, 1990. Although it is tempting to propose a formal recognition of acanthaclisines as a subfamily, it is important to underline that the study of Machado et al. The high level of polyphyly evidenced for the species-rich genus Palpares (69 valid species; see Table S1) is not unexpected, as it reflects the complex nomenclatural history of Palparinae (e.g., see Mansell, 1990), where the majority of species were first described in Palpares or assigned to this genus after having been initially described in Myrmeleon (Mansell, 1990;Stange, 2004; see also Table S1). As a result of the definition of new genera, many species of Palpares were also secondarily transferred. However, since several of these new genera were based on the study of a limited number of species (e.g., the work of Insom & Carfì, 1989, who described five new genera while only examining 12 palparine species), this inflated the polyphyly of the genus Palpares, while also creating paraphyletic genera due to the non-inclusion of several species. When examining the phylogenetic relationships of the Palpares species sampled in our study, some of them could have been predicted based on morphological similarities; for example, a group proposed by Akoudjin and Michel (2011), composed of P. kalahariensis Stitz, P. longimaculatus Akoudjin & Michel, P. incommodus (Walker) and P. radiatus Rambur, is retrieved here except for P. kalahariensis. Several Palpares species that are morphologically closely related were also recovered sister in our analyses; it is the case for instance for P. caffer (Burmeister) and P. speciosus (L.) (see Banks, 1913;Mansell & Erasmus, 2002), or for P. digitatus Gerstaecker and P. umbrosus Kolbe (see Banks, 1913;Prost, 1995;Michel, 1999). Here, it remains to be seen whether some of the inferred Palpares lineages are supported by apomorphic characters (such as the shape of distal palpomere of labium) that could justify further definition of new genera.
With 12 species, Palparellus is the second most speciose genera of Palparini. Mansell (1996) hypothesized that "Palparellus, as constituted here, may be a paraphyletic assemblage, as no autapomorphy has yet been found to confirm the monophyly of the genus". Our phy- were sampled-retrieved a Palparellus sister to a Pamexis, hence providing some support for a close relationship between these two genera. That said, this result was not expected because Pamexis were originally thought to be the sister group of Pamares, as they are morphologically quite similar; the fact that they both have reduced eyes was considered to be a synapomorphy (Mansell, 1990), but it could also be an evolutionary convergence in relation with the diurnal activity of adults. The remaining Palparellus lineages can be identified based on male genitalia (shape of the gonarcal bulla) and similarity in wing patterns (see Mansell, 1996 for details) and correspond to the following groups: (i) a first clade grouping P. damarensis (McLachlan), P. dubiosus (Péringuey), P. pulchellus (Esben-Petersen), P. ulrike Mansell and a potential new species (labelled P. pr. damarensis in our study), (ii) a second clade grouping P. astutus (Walker), P. rothschildi (Van der Weele) and P. spectrum (Rambur) and (iii) a third clade grouping P. festivus (Gerstaecker), P. flavofasciatus (McLachlan), P. nyassanus (Navás) and P. ovampoanus (Péringuey). These three clades could potentially be assigned to specific genera; this would involve restricting the genus Palparellus to the clade including the type species P. spectrum, and defining two new genera in the future for members of the remaining two clades.
Lastly, our analyses confirm the recent transfer of Nosa tristis (sensu Whittington, 2002) back to the genus Palpares by Prost (2019) who performed a comprehensive morphological revision of the genus Nosa. They also support the assignation of P. papilionoides to the genus Parapalpares by Insom and Carfì (1989).
Biogeographic history Mansell (1992aMansell ( , 1996 (Morley, 2000(Morley, , 2007, but interestingly there is little direct fossil evidence to support the existence of a rainforest biome in East Africa (Jacobs, 2004;Linder, 2017). It has been also suggested in the review of Couvreur et al. (2021) that the Zambezian bioregion was too far south, given that Africa was located about 10 south of its current location, to enable the growth of rainforest vegetation under local climatic conditions. Lastly, the inferred colonization of Madagascar by overseas dispersal, which occurred in the Paleocene, is also of particular interest as it echoes patterns observed for other insect groups such as in Sericini beetles (Eberle, Fabrizi, Lago, & Ahrens, 2017 (Carmichael et al., 2017;Handley et al., 2012;Handley, Pearson, McMillan, & Pancost, 2008;Korasidis et al., 2022). During the Middle Eocene, the plant fossil record also indicates the presence of open arid woodland formation in Eastern Africa (Jacobs & Herendeen, 2004;Linder, 2017). We can therefore assume that southern and Eastern Africa have long been suitable environments for palparine antlions.
During the Eocene-Oligocene transition (EOT), a major climate shift from a greenhouse to an icehouse climate occurred (Zachos et al., 2001(Zachos et al., , 2008 The inferred colonization route through Madagascar is not unexpected because of the well-known existence of a chain of steppingstone islands that acted as a land bridge between Asia/India and Madagascar-Africa during the Cenozoic (Warren, Strasberg, Bruggemann, Prys-Jones, & Thébaud, 2010). Biogeographic reconstructions in many groups recover this dispersal route (Sanmartín & Ronquist, 2004), including in insects (e.g., Condamine, Sperling, & Kergoat, 2013;Toussaint, Fikáček, & Short, 2016 to live in or on the fringes of desert environments (Mansell, 1990;Prost, 1995), most cannot cope with overly arid conditions. Another potential obstacle is linked to the uplifts that took place during the We would like to thank the Senior Editor Bonnie Blaimer for several constructive comments and suggestions, and Shaun L.   Figure S4. Results of RASP analyses. Table S1. Palparidiini and Palparini species list along with additional information on taxonomy and distributional information. Distributional information for almost all species is based on the comprehensive catalogue of Stange (2004) but other studies are also listed when relevant. Table S2. Taxon sampling, including the species for which we relied on GenBank data. GenBank accession numbers for seven gene fragments are provided on the right (newly generated data is highlighted using bold fonts).