The evolution of Zoraptera

Zoraptera is one of the most enigmatic and least understood orders in insects. Based on a wide taxon sampling from all continents where the group is known, we applied a phylogenetic approach using multiple DNA sequences to elucidate species‐level relationships. The resulting phylogeny shows that Zoraptera is divided into three major clades, and that two comprise species distributed on different continents. The monophyly of these clades is at least partly supported by shared derived morphological features. The divergence age estimation and ancestral distribution area reconstruction suggest an ancient origin and early radiation initiated in the Permian. Plate tectonics theory suggests that the present distribution of Zoraptera was mainly established by vicariance, rather than dispersal. The three major clades probably originated on the Pangaea supercontinent, or alternatively on the linked Gondwana and Laurasia supercontinents. Their ancient origin explains previously found conspicuous interspecific variation of the genital apparatus, sperm structure and mating behaviour, in striking contrast to a highly conserved general body morphology. We compiled data of available reproductive features and reconstructed the character evolution. Our analyses revealed repeated acquisitions and/or losses of a hyperelongated intromittent organ, mating hooks and tergal protuberances.


Introduction
Zoraptera is the third smallest order in Insecta after Mantophasmatodea and Grylloblattodea. The group is mainly distributed in subtropical and tropical regions (Grimaldi & Engel, 2005;Beutel et al., 2014;Mashimo et al., 2014c;Choe, 2018). Its phylogenetic position has been controversial (reviewed in Mashimo et al., 2014c;Kjer et al., 2016;Beutel et al., 2017), with consensus apparently reached recently with Zoraptera being placed in a monophyletic Polyneoptera based on different sources of evidence (Yoshizawa, 2011;Misof et al., 2014;Wipfler & Pass, 2014;Mashimo et al., 2014aMashimo et al., , 2015Matsumura et al., 2015;Wipfler et al., 2019). Recently, Wipfler et al. (2019) reconstructed the morphology of the common ancestor of Polyneoptera and subsequent evolutionary developments, based on a robust phylogenetic hypothesis congruent with Misof et al. (2014). They recovered Zoraptera as the sister group of Dermaptera, and both orders were placed as sister to the rest of the polyneopteran orders. They also postulated that the last common ancestor of Polyneoptera was 'a ground-dwelling insect with a largely unmodified body relative to the last common ancestor of winged insects' . Considering the winged morphs of Zoraptera (e.g. Friedrich & Beutel, 2008;Mashimo et al., 2014c;Matsumura et al., 2015), it is reasonable to postulate that extant species are relatively similar in their morphology to the aforementioned 'ancestral' polyneopteran. The combination of mostly plesiomorphic morphological features with few autapomorphies partly explains the difficulty of placing this order in a phylogenetic context (Mashimo et al., 2014c).
In contrast to the species diversity found in major polyneopteran orders (e.g. Phasmatodea, Mantodea, Blattodea, Orthoptera), to date only 44 extant species and 16 extinct species are described in Zoraptera (Mashimo et al., 2018(Mashimo et al., , 2019Chen & Su, 2019). Kukalová-Peck & Peck (1993) proposed six genera within Zoraptera based on the wing venation. However, this character complex is known to vary even within a species (Choe, 1989). Consequently, these genera were synonymized with Zorotypus by Engel (2000) who considered 'their homologies tenuous and their system unstable'. Likewise, a genus described by Chao & Chen (2000) was synonymized with Zorotypus (Engel, 2000). Since then, Engel's monogeneric classification has been widely accepted (Rafael & Engel, 2006;Terry & Whiting, 2012;Mashimo et al., 2013;Yin et al., 2015;Wang et al., 2016;Choe, 2018). The described species of Zorotypus are relatively small, typically < 3 mm, and live cryptically, mainly in rotten trees. They lack any conspicuous features with the exceptions of a very distinct dimorphism (apterous and alate morphs) and extremely varying genitalia (Mashimo et al., 2014c;Choe, 2018). Considering the assumed reconstructed common ancestor of Polyneoptera , Zoraptera have apparently acquired a lifestyle characterized by cryptic habitats, gregarious behaviour and miniaturization. However, their evolutionary origin and morphological transformations over time are still largely obscure.
The striking diversity of genital structures is in strong contrast to the external homogeneity among species. The highly diversified male genitalia have been investigated intensively, with detailed anatomical reconstructions (Hünefeld, 2007;Matsumura et al., 2014), but also in the framework of taxonomic studies (e.g. Gurney, 1938). In some species, a spiral-shaped elongated male genital structure was reported (e.g. Gurney, 1938;New, 1978New, , 2000Mashimo et al., 2013). To our best knowledge, this is a unique character state in Polyneoptera. Different types of elongation of the intromittent organs are also known within Zoraptera, looped, for instance, in Z. zimmermani (Gurney, 1939) and straight in Z. barberi (Gurney, 1938). Some species possess asymmetrical genitalia (Gurney, 1938;Paulian, 1949Paulian, , 1951Hwang, 1974Hwang, , 1976Rafael & Engel, 2006;Hünefeld, 2007;Rafael et al., 2008;Mashimo et al., 2013Mashimo et al., , 2018Wang et al., 2016;Kočárek et al., 2017;Yin & Li, 2017), in some species remarkably differing in their specific features. Another conspicuously diversified character system is the structure of elements of the terminal abdominal segment, the presence or absence and size of spines and mating hooks (e.g. Gurney, 1938). This variation may be related to the mating posture to some extent. Males are probably connected to females by these terminal structures, laying in an upside-down position in the majority of the species in which the mating posture is known (Shetlar, 1978;Choe, 1994Choe, , 1995Mashimo et al., 2011;Dallai et al., 2013). The disparity between a far-reaching uniformity of the general morphology and an extreme diversity of genital features is obviously a fascinating phenomenon and a challenging topic in evolutionary biology. However, a reliable evaluation has not been possible so far due to the lack of formal phylogenetic analyses on the species level (see Engel, 2003).
The primary aim of our study is to reconstruct the phylogeny within Zoraptera using molecular data. The taxon sampling covers species from all continents. The evolutionary history is evaluated by means of divergent time estimation and based on the plate tectonics theory (see Seton et al., 2012). The character evolution with a special focus on reproductive structures was reconstructed based on the phylogenetic trees.

Materials and methods
Most of the specimens were collected for this study and fixed with 80-99.5% ethanol (Fig. 1). Type specimens of Zorotypus novobritannicus were borrowed from the Arthropod Collection, Brigham Young University, Provo, UT, U.S.A. In total, 31 individuals belonging to 21-22 species were included (Fig. 2). We failed to trace the identity of sample YK16-10 collected from Ecuador and named it Z. sp. 6. We collected only Z. huxleyi and Z. hamiltoni from the same locality at the same time, and Z. sp. 6 probably belongs to the latter species. Outgroups were selected from all polyneopteran orders and some species from Psocodea, Hemiptera, Ephemeroptera, Odonata, Zygentoma and Archaeognatha. The tree was rooted with Archaeognatha. The sequences of the outgroup taxa were obtained from GenBank (metadata of the samples are listed in Tables S1, S2).

Model selection and phylogenetic estimation
The best substitution models and partition schemes were estimated using partitionfinder 2.3.3 (Lanfear et al., 2017), with the greedy algorithm. The codon positions for Histone 3 (three partitions) and rRNA (three partitions) were predefined for the partitionfinder analyses. The best-fit partition scheme and models were described in the nexus formatted data matrix (File S1).
Previous studies showed that the phylogenetic relationships of polyneopteran orders cannot be estimated accurately using a limited number of gene sequences (e.g. Kjer, 2004;Yoshizawa & Johnson, 2005;Ishiwata et al., 2011;Misof et al., 2014;Wipfler et al., 2019). Therefore, we constrained the phylogenetic relationships among orders according to Misof et al., (2014) and Wipfler et al., (2019) for the following phylogenetic analyses. In addition, unconstraint analyses were also performed to test the monophyly of Zoraptera (see File S1). We estimated a maximum likelihood tree using iq-tree 1.6.3 (Nguyen et al., 2015), with 10 000 replicates of an ultrafast likelihood bootstrap with -bnni option (Hoang et al., 2018) to obtain bootstrap branch support values. To see the stability of results, ten independent iq-tree analyses were performed. All the analyses resulted in a concordant result, and we selected the tree obtained from the last run for Fig. 3. A Bayesian analysis was performed using mrbayes (Ronquist & Huelsenbeck, 2003). We performed two runs each with four chains for 1 000 000 generations, and trees were sampled every 1000 generations. The first 10% of sampled trees was excluded as burn-in, and a 50% majority-consensus tree was computed to estimate posterior probabilities. tracer in the beast software package (Bouckaert et al., 2014) was used to check that the Markov chain Monte Carlo (MCMC) runs reached a state of convergence.

Divergence time estimation
For divergence date estimation, a Bayesian method was adopted using the software mcmctree in the paml 4.8 software package (Yang, 2007) and beast 2.6 (Bouckaert et al., 2014). The following fossil ages were used as soft minimum bounds according to Misof et al. (2014) and Tong et al. (2015): 160 Ma for the deepest divergence of Plecoptera; 130 Ma for the deepest divergence of Isoptera; and 99 Ma for the deepest divergence of Embioptera (Table 1). For all fossil calibrations, the age of Rhynie chert (412 Ma) was used as the hard maximum bound according to Evangelista et al. (2019). In addition, the hard maximum bound 450 Ma was also adopted for Zygentoma-Pterygota divergence age according to Misof et al. (2014) and Tong et al. (2015).
For the mcmctree analysis, we first estimated the substitution rate prior using the divergence date 419 Ma for the Polyneoptera-Paraneoptera branching according to Tong et al. (2015). Based on the result, a gamma prior for the substitution rate was estimated using baseml in the paml software package. The GTR + G model was adopted with = 0.5, which was a close approximation of the best substitution model estimated by jmodeltest (Posada, 2008) for the entire dataset. We performed a run for 1 000 000 generations, and the values were sampled every 50 generations. The first 10% of the obtained values were excluded for burn-in. We ran two independent analyses to check that the MCMC runs reached a state of convergence (dos Reis et al., 2017).
For the beast analysis, we used the clade ages package (Matschiner et al., 2017). The following options were selected: BEAST Model Test for the site model, relaxed clock log normal for the clock model, and birth-death model for the priors. We performed a run for 20 million generations, and the first 10% of the obtained values were excluded for burn-in. tracer in the

Biogeographical analysis
Ancestral area reconstruction was performed using the dated tree obtained from the mcmctree analysis. Outgroup samples were excluded from the analysis. We used a dispersal-extinction-cladogenesis (DEC) model (from Lagrange; Ree & Smith, 2008) as implemented in the software rasp 3.2 (Yu et al., 2015). Dispersal-vicariance analysis (DIVA; Ronquist, 1997) was a potential alternative to the DEC model. However, a previous study showed that DIVA wrongly identifies ancestral areas with complex patterns of dispersals and within-area speciation events (Kodandaramaiah, 2010). Five geographical realms were defined: Afrotropical (AF), Indomalaysian (IM), Nearctic (NA), Neotropical (NT) and Australasian. The maximum number of areas allowed for ancestral distributions at each node was set to two, and dispersal between all pairs of distributional areas was equally weighted. For extant species, there are no species distributed in two or more biogeographical regions. The biogeographical region coding of each sample was based on the known distributional range of the species.

Character evolution
We focused on the following eight features: (i) absence or presence of an elongated intromittent organ; (ii) symmetry of genitalia; (iii) absence or presence of basal plate; (iv) absence or presence and size of mating hook; (v) absence or presence of records of males; (vi) absence or presence and size of protuberances on abdominal tergites 9-11; (vii) modifications of subgenital plate; and (viii) absence or presence of hairy patch on vertex ('fontanelle gland' in the literature). The relevant data are provided in many taxonomic studies, and most of the information was obtained from the literature. If necessary, we examined specimens under a stereomicroscope Olympus SZX12 (Olympus Corporation, Tokyo, Japan) to obtain additional information. For an overview of the diversity focus stacking images of the caudal view were taken using a stereomicroscope Leica M205 A equipped with a Leica DFC420 camera and the software las 3.8 (Leica Microscopy GmbH, Wetzlar, Germany). Relevant information for the outgroups was obtained from the following studies: Tuxen (1970), Helm et al. (2011) and Klass et al. (2003). Character evolution was reconstructed with the software mesquite 3.6 (Maddison & Maddison, 2018).

Data accessibility
All supporting data are available electronically in File S1.

Phylogeny, dating and biogeography
The aligned sequences consisted of 1753 bp (of which 75 bp were excluded from the analyses), and the obtained maximum likelihood (ML) and Bayesian trees (Figs 3, S1-S3) were congruent except for one weakly supported branch. Neither analysis, with and without constraining phylogenetic relationship among outgroup orders (Misof et al., 2014;Wipfler et al., 2019), yielded different phylogenetic relationships within Zoraptera (Figs 3, S1-S3). Therefore, the influence of the analytical methods on the evolutionary history discussion is negligible, and we focused on the ML tree thereafter.
The divergence age estimation based on a beast model (Fig. S4) showed very similar results. The divergence ages of each node were estimated as slightly older than those obtained from a mcmctree analysis, while the estimated divergence ages of clade 1 were younger. However, in all cases they largely overlapped.

Evolution of reproductive character states
Based on our original observations and a literature survey, mainly from taxonomic studies, information on the eight characters described earlier related to reproduction was available for the majority of the described species and is summarized in Table S3. Consensus on structural homologization in zorapteran genitalia was not evident, and researchers used different terminologies based on varying interpretations in the literature. A reliable assessment of homologies is still pending, a point also emphasized in Boudinot's (2018) first comprehensive synthesis of insect male genitalia. As this turned out to be too ambiguous for characters (3) and (6), they were not scored for terminals outside of Zoraptera.
The ML reconstruction was performed to estimate evolutionary histories of the eight features (Fig. 5). Hyperelongation of the intromittent organ occurred at least twice, possibly even three times (Fig. 5A). As clade 2a includes many species with unknown character states, the independent occurrence of the spiral in clades 2a and 2b is not conclusive. An asymmetric configuration of the genitalia is a commonly observed feature in Polyneoptera (Huber et al., 2007), but a symmetrical condition is considered as ancestral for the group (Helm et al., 2011;Boudinot, 2018). The ancestral state of Zoraptera was probably symmetric according to our estimate (Fig. 5B).
The asymmetric state was probably acquired in the ancestor of clade 1.
The last common ancestor of Zoraptera had probably acquired a mating hook, or at least a small mating hook was present in the common ancestor of clades 2 + 3. Its enlargement and modifications occurred independently in several lineages (Fig. 5C). Modifications of the marginal area of the subgenital plate are  known in some species that form the monophyletic subunit Z. shannoni (NT) + Z. asymmetristernum (AF) + Z. sp. 1 (AF) of clade 1 (Fig. 5D). It is likely that a bifurcated margin has evolved in the common ancestor, with subsequent transformation of the subgenital plate in Z. asymmetristernum.
The presence or absence of a basal plate in the zorapteran groundplan remains equivocal (Fig. 5E), as the identity of the basal plate observed in species of clades 2 and 3 is not yet confirmed. However, it is clearly shown that the character states separate clades 1 from 2 + 3. This implies that this character was either lost or completely modified in clade 1, or alternatively newly developed in clade 2 + 3. The high diversity among species is usually visible in caudal view of the abdomen (Fig. 6). Easily visible differences among species are the presence or absence of protuberances of tergites 9-11, and also different sizes of these structures. The homology of these protuberances was completely unclear, and the evolutionary history was reconstructed based on data only acquired from Zoraptera. Varying conditions of these surface structures were found in different clades, and it seems that repeated acquisitions and losses occurred in Zoraptera (Fig. 5F).
The lack of any records of males does not necessarily mean that the concerned species are parthenogenetic. However, it is confirmed that females of Z. brasiliensis (NT, c2c) and Z. gurneyi (NT) can reproduce parthenogenetically (Silvestri, 1947;Choe, 1997). Our results must be considered as preliminary. However, it appears that species with unknown males are not closely related, suggesting possible independent losses (Fig. 5G), with parthenogenesis evolving several times independently in Zoraptera. Courtship feeding through a hairy patch on the vertex is known in Z. barberi (NT, c3) (Choe, 1995). The setal patches are known from several additional species (Table S3). Although the presence of a gland is not confirmed, we mapped the externally visible character state on the phylogeny. The ancestral state was probably absent, and it appears likely that acquisitions occurred repeatedly (Fig. 5H).

Discussion
The present study confirms the monophyly of Zoraptera by means of formal phylogenetic analyses based on a broad sampling of zorapteran species, as previously shown by Yoshizawa & Johnson (2005). Traditionally, the following features including reductions were considered as potential autapomorphies of the order Zoraptera: (i) distinct dimorphism (apterous and alate morphs); (ii) strongly simplified wing venation and a capability of shedding the wings; (iii) two-segmented tarsi without adhesive structures; and (iv) correlations of presence/absence of compound eyes, ocelli and distinct pigmentation (Beutel & Gorb, 2001;Beutel et al., 2014;Mashimo et al., 2014c). Holocentric chromosomes reported from Z. hubbardi (NA) (Kuznetsova et al., 2002) could be an additional autapomorphy of Zoraptera. Comparing with the reconstructed ground-dwelling ancestor of Polyneoptera , which intuitively resembles a grasshopper, distinct miniaturization and partial structural simplification must have occurred in the common ancestor of Zoraptera, possibly due to the habitat specialization, a preference for subcortical spaces (under bark) of fallen trees where spatial size is extremely limited. Our molecular phylogenetic approach revealed the further evolutionary history of Zoraptera.

Phylogeny, dating and biogeography
The results of our analyses of molecular data suggest that extant Zoraptera form three major clades. The early splits presumably occurred in the early Permian (Fig. 3) or possibly the Carboniferous (Fig. S5), when the continents were united as Pangaea, or at least a connection existed between the supercontinents Gondwana and Laurasia (Smith et al., 2004). The heterogeneous distribution ranges found in clades 1 and 2 may be mainly due to their old origin and subsequent vicariance between the contemporary continents.
Most species of clade 1, i.e. Z. hubbardi (NA), Z. impolitus (IM), Z. shannoni (NT), Z. sp.1 (AF) and Z. asymmetristernum (AF), presumably originated before the continents rifted around 80-100 Ma (Seton et al., 2012). The split of Z. hubbardi (NA) and Z. impolitus (IM) occurred in the early Jurassic, when Laurasia still existed (Seton et al., 2012). The rest of clade 1 diverged in the late Cretaceous, and the ancestral distribution was estimated to be Afrotropic + Neotropic, which corresponds to Gondwana (Seton et al., 2012). The break-up of Pangaea probably took place in the Jurassic and Cretaceous (100-160 Ma; Seton et al., 2012), and this may explain the split of the two lineages of clade 1. The split of Z. shannoni (NT) and Z. asymmetristernum (AF) (85 Ma) is possibly also a result of vicariance. Although the contemporary Neotropical and Afrotropical regions had probably rifted c. 100 Ma (Seton et al., 2012), the estimated divergence age contains an estimation error. However, it is also possible that the two species arose in the Afrotropical region and Z. shannoni dispersed by drifting through the South Atlantic Ocean, and indeed the divergence age estimated by a beast analysis was relatively young (24 Ma).
The split between clades 2 and 3 was estimated at c. 236 Ma, and clade 2 diverged c. 183 Ma. The first split probably occurred in the regions corresponding to the contemporary Neotropical region. This happened before the separation of Gondwana and Laurasia (Smith et al., 2004;Seton et al., 2012). Therefore, it is conceivable that the ancestral species was distributed in the corresponding southern part of Pangaea. Clade 2 includes three major lineages (c2a, c2b, c2c), and clade 2b comprises solely Indomalaysian and Australasian species. This lineage presumably arose 161 Ma and diverged 138 Ma. During this period, it is assumed that the Indian subcontinent + Australasia started to rift from Gondwana, and Australasia started to separate from the Indian subcontinent 120 Ma (Seton et al., 2012). The time of the zorapteran radiation and continental break-up is an estimation and prone to errors. Our estimations did not always suggest that speciation and lineage splits occurred before the estimated continental break-up. However, considering the comprehensive information mentioned earlier, the species distribution can be explained by vicariance rather than dispersal. Consequently, it can be assumed that the speciation mainly occurred on individual continents. The only exception among the studied species is Z. mexicanus (NA, c2a), whose origin was dated as 30 Ma. The formation of the Panama-Costa Rica Arc is estimated around 60-90 Ma (different hypotheses are discussed in Seton et al., 2012). This suggests that Z. mexicanus (NA) is derived from the South American lineage, which invaded Central America after the formation of the Arc.
The early Permian origin of the major zorapteran lineages explains the enormous disparity of the mating behaviour (Choe, 1994(Choe, , 1995Dallai et al., 2013) as well as the impressive variation in the genitalia and sperm morphology among species (Dallai et al., , 2012. The recently observed external sperm transfer of Z. impolitus (IM, c1)  was the first report for a pterygote insect. The exceptional divergence of characters linked to reproduction stands in stark contrast to a far-reaching uniformity in the general body morphology, which is preserved since the late Palaeozoic for reasons not yet understood. Similar phenomena did not evolve in other groups with an origin in the same period, e.g. in the presumptive sister taxon Dermaptera. Sperm morphology can be useful for estimating phylogenetic relationships in some cases (Gottardo et al., 2016). Dallai et al. (2014a,b) proposed a hypothesis for sperm evolution in Zoraptera, suggesting that: (i) those species used in Dallai et al. (2011Dallai et al. ( , 2012Dallai et al. ( , 2014a arose before the fragmentation of Gondwana in the mid-Cretaceous period (the time was estimated due to available fossil records, e.g. Poinar Jr, 1988;Engel & Grimaldi, 2002;Kaddumi, 2005); (ii) Z. caudelli (IM, c2b), Z. magnicaudelli (IM, c2b), Z. huxleyi (NT, c2c) and Z. weidneri (NT, c2c) form a monophyletic unit; and (iii) that Z. shannoni (NT, c1), Z. hubbardi (NT, c1) and Z. impolitus (IM, c1) definitely belong to different lineages. Mashimo et al. (2015) also found a possible synapomorphic feature of eggs of Z. impolitus (IM, c1) and Z. hubbardi (NA, c1). These interpretations are congruent with our molecular phylogeny and corroborate our evolutionary hypotheses.

Classification
All extant species of Zoraptera are now classified under the single genus, Zorotypus. The ancient origin, genetic divergence and the unusual diversity of genitalia and sperm arguably suggest a division into several supraspecific subunits. As mentioned in the Introduction, Kukalová-Peck & Peck (1993) established six genera based on the wing venation and biogeographic distribution (Old vs New World). However, only one Old World zorapteran was included in their study, and phylogenetic relationships among the studied species was not reconstructed with a formal approach (Kukalová-Peck & Peck, 1993). In addition to this, Chao & Chen (2000) established another Old World genus from Taiwan due to an unusual appearance. The taxonomic treatments of Kukalová-Peck & Peck (1993) and Chao & Chen (2000) did not meet the criteria for the erection of supraspecific taxa outlined by Komarek & Beutel (2006), especially the claim that all supraspecific units (not only the newly erected one) should be monophyletic. Our results also clearly reject the idea that the biogeographic distribution of Zoraptera is useful for the classification of the order. Nevertheless, the concept of Kukalová-Peck & Peck (1993) appears at least partly justified. Four of the six included species were also analysed in our study, i.e. Z. barberi (NT, c3), Z. brasiliensis (NT, c2c), Z. caudelli (IM, c2b) and Z. hubbardi (NA, c1), each of them designated as type species of a separate genus in Kukalová-Peck & Peck (1993). These species were each recovered in different lineages in our analyses, and clade 1, including Z. hubbardi (NA, c1), showed specific trends of character state evolution as discussed in detail in the following. Clade 1 features asymmetric genitalia, without a hyperelongated intromittent organ and without a basal plate. The asymmetric condition of this subunit is apparently an autapomorphy, and the remaining species are also monophyletic. Therefore, we consider it a potential option to resurrect one of the genera proposed by Kukalová-Peck & Peck (1993).
For a further taxonomic step, the position of the type species of the genus Zorotypus, i.e. Z. guineensis (AF, not included here), would have to be clarified. The original description of Silvestri (1913) is relatively detailed, but we could not extract sufficient information from it. Dallai et al. (2014b) re-evaluated Silvestri's study and the original histological samples. They confirmed that males lack an elongated intromittent organ. Although it is not explicitly mentioned, the figures show neither asymmetric genital sclerites found in species of clade 1 nor any basal plate-like structure typical of species of clade 2 + 3. Furthermore, Dallai et al. (2014b) discussed that the male reproductive system (documented with histological sections) displays features that are probably similar to conditions found in Z. magnicaudelli (IM, c2b), Z. caudelli (IM, c2b) and Z. huxleyi (NT,c2c), rather than in Z. hubbardi (NA, c1) and Z. impolitus (IM,c1). However, there are also features resembling those of species of clade 1. For instance, features of the hind femur of Z. guineensis are very similar to those found in Z. shannoni (NT, c1) (see Silvestri, 1913;Gurney, 1938), and similarly the hairy area on the vertex. Based on the evidence at hand, Z. guineensis (AF) could belong to any clade recognized in our study. Therefore, we refrain from further taxonomic steps in our study.

Morphological evolution
Our analyses suggest independent origins of hyperelongated intromittent organs, and that these derived states originated from symmetric genitalia. From a morphological point of view, this is also supported by obvious differences between the straight elongated intromittent organ of Z. barberi (NT, c3) (Gurney, 1938) and the spiral-shaped elongated one found in many species of Zoraptera (Table S3). However, considering the very specific and complex anatomy of the male genital apparatus of Z. caudelli (IM, c2b), studied in detail by Matsumura et al. (2014), it appears unlikely that the type with a complex, spiral-shaped element has evolved several times independently. The entire apparatus is exceptionally complicated, with structures of unclear homology [e.g. Z. hubbardi (NA, c1)] (Hünefeld, 2007). Moreover, it is highly unlikely that complex structures with very specific and complicated configurations have repeatedly evolved in the same way in different species. The ambiguity of the scenario is increased by species of clades 2a and 2b with males not known yet. Recently, Rafael et al. (2017) reported a gynandromorph of Z. brasiliensis (NT, c2c), whose males were previously unknown, containing both male and female characteristics. They discussed possible thelytokous parthenogenesis, with unfertilized eggs yielding females but not males. As another example, Z. gurneyi produces males, but parthenogenetically reproducing populations also occur (Choe, 1997). As Choe (2018) stated in a recent review, divergent mating systems are exhibited even between sympatric species. Therefore, it is still debatable whether the spiral was present or absent in the groundplan of clade 2a + 2b. Additional laboratoryand field-based observations of reproductive modes from different populations are necessary, as pointed out by Rafael et al. (2017).
Asymmetric genitalia appear to have evolved in clade 1. The homology of the sclerites composing this type of genitalia is not yet clarified. However, using the available literature (Table S3) we identified 14 out of 28 species with an asymmetric genital apparatus. In the present study we treated asymmetric genitalia as one category, even though structural differences were reported between Z. hubbardi (NA, c1) and Z. shannoni (NT, c1) (illustrated in Gurney, 1938). Detailed morphological data on male genitalia, including musculature and related membranes, are available for only two species, Z. hubbardi and Z. caudelli (Hünefeld, 2007;Matsumura et al., 2014). Although zorapterans are rather small, technical problems caused by size reduction play a minor role in state-of-the-art insect anatomy, if at all . New detailed anatomical studies will probably help to clarify homology issues, and in a second step to unveil the evolution of the genital structures. This also applies to symmetric genitalia, to clarify not only the origin of the elongated intromittent organs, but also the homology of the basal plate. This issue is related to the challenging interpretation of the tergite numbering in Zoraptera, with distinctly different interpretations suggested by various authors, as shown in columns B and F of Table S3. Mashimo et al. (2014a,b) elegantly established the tergite numbering for Z. caudelli (IM,c2b). Therefore, this issue should be relatively easy to solve by carefully comparing abdominal segments for representative species.
Structural diversity as typically seen in a caudal view of the abdomen has seemingly evolved through repeated development, retrogress or loss of the mating hook and protuberances on tergites 9-11. Presence of the mating hook is estimated as a possible plesiomorphic state in our analyses. Although its function is not yet known, it may indeed work as a hook during copulation. Except for the external sperm transfer of Z. impolitus (IM, c1), the known mating posture is that males are coupled to a female through the genitalia and lay upside down [Z. hubbardi: NA, c1 (Gurney, 1938); Z. barberi: NT, c3 (Choe, 1995); Z. gurnery: NT, not included (Choe, 1994); Z. huxleyi and Z. weidneri: NT, c2c (J.A. Rafael et al., personal observation); Z. caudelli: IM, c2b ; Z. magnicaudelli: IM, c2b ]. Any clasping structure or hook would probably be helpful for this type of mating posture. Therefore, it is surprising that Z. impolitus has one of the most developed mating hooks , although they attach sperm to females externally . The mating hooks might have an additional function, e.g. opening the female genitalia forcefully to deposit a spermatophore in the female genital tract. This needs verification by detailed observations of intertwining male and female genitalia. In addition, characteristic mating behaviour is also reported in Apachyus chartaceus (Dermaptera: Apachyidae) (Shimizu & Machida, 2011), and it may prove worthwhile to compare their genital coupling with that of zorapterans.
A function of the hairy patch on the male vertex in Z. barberi (NT, c3) is secreting nutritious fluid for the females as a nuptial gift (Choe, 1995). Superficially similar structures were observed in some studied species, and also in five out of 25 species with information available in the literature (Table S3). Studies using histological sections of the head of the relevant species are necessary to clarify the presence or absence of gland tissue to confirm its function in the other species.
The knowledge on Zoraptera has increased rapidly in the last decade. However, to further understand the evolution of the group, additional investigations are necessary. Future studies should have a main focus on the detailed morphology of genital organs, interactions of the male and female genitalia, mating behaviour of each species, and sperm morphology of representative species from additional early split clades, also including the type species Z. guineensis.

Author contributions
YM, RGB, KY and JAR, and JTC conceptualized independently; JAR, JTC, YM and KY collected samples; YM, IY and KY performed molecular experiments; KY analysed molecular data; YM and KY analysed evolutionary histories; and all authors interpreted the results. YM and KY wrote the draft and RGB, JAR, SPL, JTC and IY revised it. All authors contributed to and approved the final draft of the manuscript.