Increasing 28 mitogenomes of Ephemeroptera, Odonata and Plecoptera support the Chiastomyaria hypothesis with three different outgroup combinations

Ethical statement

No invertebrate specimens used in this study are protected under the provisions of the laws of People’s Republic of China on the protection of wildlife. Hence, there is no ethical problem with animal sampling. The study protocol was reviewed and approved by the Committee of Animal Research Ethics of Zhejiang Normal University.

Samples and sequencing

Mitochondrial genomes from 28 species of mayflies, dragonflies, damselflies and stoneflies were sequenced. All samples were collected between 2005 and 2014; information on species collected are available in Table 1. All mayflies and stoneflies were collected at the larval stage, whereas the dragonflies (or damselflies) were collected at the adult stage. All samples were preserved in 100% or 85% ethanol and stored at −40 °C in the Zhang‘s lab, College of Life Science and Chemistry, Zhejiang Normal University, China. All samples were identified by JY Zhang based on morphological characteristics. DNA was extracted from either a whole individual (for smaller species) or from the legs of an individual (for larger species) using DNeasy Tissue Kits (Qiagen, Hilden, Germany). For each sample we amplified short 400–1,500 bp mt gene fragments using degenerate primers designed to match specific arthropod mt genes according to the method of Simon et al. (2006). PCR reaction conditions and procedures were as described in Zhang, Song & Zhou (2008a), Zhang et al. (2008b), Yu et al. (2016). Amplification parameters were not stringent (46−58 °C annealing temperatures). Short fragments were sequenced using the ABI 3730 system with bidirectional sequencing. Species-specific primers were then designed from these short fragments and used to amplify longer 2–3 kb segments of the mt genome from each sample. Long PCR products were sequenced directly using a primer-walking strategy initiated with specific primers with bidirectional sequencing. Raw sequence files were proofread manually and assembly of the nearly complete genome sequence was performed in SeqMan of Lasergene version 5.0 (Burland, 1999). Thirteen protein-coding genes (PCGs) and 2 rRNA (16S rRNA and 12S rRNA) genes were identified by comparisons with homologous sequences of known insect mt genes. Most tRNA genes were identified by their cloverleaf secondary structure using tRNAscan-SE 1.21 (Lowe & Eddy, 1997). If selected tRNA genes could not be determined by tRNAscan-SE, they were identified by comparison with the homologous insect tRNA genes.

Taxa and alignment

In previous phylogenetic analyses some insect species have been reported to suffer from long-branch attraction and to have unstable phylogenetic positions (Nardi et al., 2003). Although many mt genomes from all orders of Insecta can now be acquired from GenBank data, we did not include data from those orders whose species showed long-branch attraction in previous phylogenetic analyses or high gene rearrangements or high base compositional biases which can cause erroneous results (Bergsten, 2005; Castro & Dowton, 2005; Ma et al., 2014; Wan et al., 2012). To reduce the computational analysis time we choose 1–5 species of the Diptera, Dermaptera, Lepidoptera, Orthoptera, Mantodea, Phasmida, Blattodea, Isoptera (=Blattodea: Termitoidae), Megaloptera, Neuroptera, Mecoptera, Grylloblattodea, Mantophasmatodea and Coleoptera as the ingroup. The outgroups selected were the four available bristletail sequences (Archaeognatha) ((Nesomachilis australica (Cameron et al., 2004), Pedetontus silvestrii (Zhang, Song & Zhou, 2008a), Petrobius brevistylis (Podsiadlowski, 2006), Trigoniophthalmus alternatus (Carapelli et al., 2007)) and three silverfishes (Thysanura) (Atelura formicaria (Comandi et al., 2009)), Tricholepidion gertschi (Nardi et al., 2003), Thermobia domestica (Cook, Yue & Akam, 2005)). We combined the 28 newly-sequenced mt genomes with 85 species from 19 orders of Insecta obtained from GenBank (Table 2). The ATP8 gene was not used in the subsequent analyses due to its shortness (about 50 amino acid residues) and its poor conservation (Zhang et al., 2008b). We aligned amino acid sequences of each mt protein-encoding gene separately. The dataset was composed of 113 taxa, including 4 bristletails, 3 springtails, 13 odonatans, 16 mayflies, 10 stoneflies and 39 neopteran species, all retrieved from GenBank (Table 2). We used two trivially different alignment methods (Cluster W and Muscle in Mega 7.0 (Kumar, Stecher & Tamura, 2016)). To reduce bias we deleted the nonconserved regions in each of gene alignment using GBlocks 0.91b (Castresana, 2000) with block parameters set as default. The two alignments and block identification procedures ensured that the conserved regions were reliably homologous. Each of the conserved alignments were concatenated for both nucleotide and amino acid translated datasets, a final alignment 7011 nucleotides and 2337 amino acids residues. A saturation analysis of the concatenated DNA was performed for first, second and third codon positions using DAMBE4.2.13 (Xia & Xie, 2001). Third codon positions were saturated so they were excluded from the final alignment and the final alignment 4674 including first and second codon positions of 113 sequences were obtained (named 113 dataset).

To compare the effect of outgroup choice, we added and analyzed two additional combinations of outgroup species: 114 dataset selecting five bristletails (Nesomachilis australica, Pedetontus silvestrii, Petrobius brevistylis, Trigoniophthalmus alternatus, and Petrobiellus puerensis) as outgroups according to the Metapterygota hypothesis supported by Cai et al. (2018); 119 dataset selecting five diplurans (Campodea lubbocki, C. fragilis, Japyx solifugus, Lepidocampa weberi, and Occasjapyx japonicus) as outgroups according to the Palaeoptera hypothesis supported by Song et al. (2019). The conserved region and saturation analyses are similar to the method used for the 113 dataset. The dataset with five bristletails as outgroups and the dataset with five diplurans as outgroups were named the 114 dataset and 119 dataset, respectively.

Phylogenetic analyses

PartitionFinder ver. 1.1.1 (Lanfear et al., 2012) was used to select partitions using the Bayesian Information Criterion (BIC). The 24 parts of three nucleotide datasets were divided into 9 partitions. The best model of three datasets was identified as a general time-reversible model (GTR+I+G) in all partitions excluding a TVM+I+ model in the second position of the COIII gene partition, whereas the second best model of this was the GTR+I+G using ModelGenerator v0.85 (Keane et al., 2006). To deduce differences in phylogenetic analyses we used a GTR+I+G model with first and second positions of 12 mt genes in three nucleotide datasets. We used MrBayes v. 3.2 (Ronquist et al., 2012) to estimate phylogenies by Bayesian inference (BI) for each dataset. Two independent runs of four incrementally heated Markov chain Monte Carlo (MCMC) chains (one cold chain and three hot chains) were simultaneously run for ten million generations depending on the datasets, with sampling conducted every 100 generations. The convergence of MCMC, which was monitored by the average standard deviation of split frequencies, reached below 0.01 within ten million generations depending on the dataset, and the initial 25% of the sampled trees were discarded as burn-in. Bayesian posterior probabilities were calculated from the remaining set of trees. Maximum likelihood (ML) analysis was conducted using PhyML (Guindon et al., 2005), specifying four substitution rate and the corresponding BI tree as the start tree. The confidence values of the ML tree were evaluated via a bootstrap test with 100 iterations.

Genome organization of mtDNA

Of the 28 complete or nearly mt genomes of Ephemeroptera, Odonata and Plecoptera, four species of Ephemeroptera were found to contain a gene rearrangement, and four species of Ephemeroptera contained a tRNA gene duplication (Fig. 2). Four species from Heptageniidae (Cinygmina obliquistrita, Epeorus sp., Heptagenia sp. and Iron sp.) had an extra tRNA-Met gene located between tRNA-Ile and tRNA-Gln. Four species from Ephemerellidae (Drunella sp., Serratella sp., Uracanthella sp. and Ephemerellidae sp.) had a tRNA-Ile gene translocated and inversed to a position between 12S RNA and the control region. One copy of tRNA-Ile was present in Drunella sp. and Serratella sp., whereas two copies of this gene were present in Uracanthella sp. and three copies in Ephemerellidae sp.

Phylogenetic analyses of three datasets

We employed BI and ML analyses to construct phylogenetic trees from three datasets with different outgroups (Fig. 3 and Figs. S1–S4). We found that different outgroups can affect the relationships between Ephemeroptera, Odonata and Neoptera. The Chiastomyaria hypothesis was strongly supported in BI analyses of three outgroup datasets and ML analysis of the 113 dataset (Fig. 3 and Figs. S1–S2), whereas different results in ML analyses of the 114 and 119 datasets were found (Figs. S3–S4). ML analyses in 114 and 119 datasets failed to support the Chiastomyaria hypothesis. No support for either the Palaeoptera or Metapterygota hypotheses was found in the BI and ML analyses in this study.

Comparing the phylogenetic trees derived from the three datasets with different outgroups, we found that the 113 dataset using 4 bristletails and 3 silverfishes as the outgroups is the stable phylogeny with high posterior probability and bootstrap value. In the 113 dataset, the monophyly of Odonata, Ephemeroptera and Neoptera was strongly supported in both BI and ML analyses (Fig. 3). Odonata as the sister of the remaining Pterygota had high posterior probability (1.00) and bootstrap value (100). Ephemeroptera as the sister of Neoptera also showed high support (BI/ML: 1.00/100) and Plecoptera was sister to the remaining other Neoptera with high support (0.92/100).

Within Odonata, the monophyly of the suborders Zygoptera and Anisoptera was well supported. At the family level, the monophyly of Gomphidae, Calopterygidae, Euphaeidae, and Libellulidae was well supported but the monophyly of Coenagrionidae was not supported since Platycnemididae clustered within Coenagrionidae.

Within Ephemeroptera, Siphluriscus chinensis (Siphluriscidae) was sister to the other mayflies, which were divided into two clades. All nodes had high posterior probabilities excluding (Heptageniidae + Isonychiidae) which was 0.59. The monophyly of all families excluding Siphlonuridae was well supported. The monophyly of Siphlonurus (Siphlonuridae) was not supported because Siphlonurus immanis (Siphlonuridae) (FJ606783) clustered within Ephemera (Ephemeridae) whereas Siphlonurus sp. (Siphlonuridae) was a sister clade to Ameletus sp. (Ameletidae).

Within Neoptera, the monophyly of Plecoptera, Mantodea, Blattodea (including Isoptera), Diptera, Coleoptera, Neuroptera, Megaloptera, Orthoptera, Trichoptera, Lepidoptera, Hemiptera, and Dermaptera were all supported. Polyneoptera was not monophyletic as Plecoptera was sister to the remaining Neoptera and in turn Dermaptera was sister to Neoptera except Plecoptera. Orthoptera failed to cluster with the other orders typically assigned to Polyneoptera and was sister to Amphiesmenoptera (Lepidoptera + Trichoptera). The relationship of (Plecoptera + (Dermaptera + other Neoptera)) was strongly supported. The other two datasets generated from the different outgroups in BI analyses (Figs. S1–S2) further supported the relationship of (Odonata + (Ephemeroptera + (Plecoptera + other orders of Neoptera)). All the topology of inter-order and intra-order in BI analyses of 114 and 119 datasets was coincident with the 113 dataset except posterior probabilities and branch length so we just showed the topology of intra-order level.

The mtDNA rearrangement and rearrangement mechanisms

In the previously published mt genomes of Ephemeroptera, an extra tRNA-Met gene was found in the heptageniids, Parafronurus youi (Zhang et al., 2008b), Epeorus sp. (Tang et al., 2014), and Epeorus herklotsi (Gao et al., 2018), but was not found in Paegniodes cupulatus (Zhou et al., 2016). In the current study, Cinygmina obliquistrita, Epeorus sp., Heptagenia sp. and Iron sp. (Heptageniidae) all had the extra tRNA-Met. Based on the all species excluding P. cupulatus possessing an extra tRNA-Met (Fig. 3), we deduce that the common ancestor of the Heptageniidae had an extra tRNA-Met gene, which can be explained as well as P. youi through the tandem duplication-random loss (TDRL) model (Zhang et al., 2008b). However, the extra tRNA-Met gene in P. cupulatus was apparently deleted in an independent random loss thereby restoring the IQM tRNA gene cluster (Fig. 4).

The transposition and inversion of tRNA-Ile genes in the mt genome of the Ephemerellidae can also be explained by the recombination and tandem duplication-random loss (TDRL) model (Fig. 5). Firstly, the tRNA-Ile inversion in Ephemerellidae may have been caused by a recombination in the original location, involving the breakage and rejoining of participating tRNA-Ile (Lunt & Hyman, 1997). Secondly, a tRNA-Ile inversion-CR arrangement probably occurred by duplication of the CR-tRNA-Ile inversion regions, resulting in a CR-tRNA-Ile inversion-CR-tRNA-Ile inversion arrangement, and then deleted by a subsequent random loss of the first copy CR and the second copy tRNA-Ile gene. At this point, it can be explained by a typical TDRL model (Moritz, Dowling & Brown, 1987). One copy of the tRNA-Ile gene inversion was formed in the ancestor of Drunella sp. and Serratella sp. as well as Ephemerella sp. (Ephemerellidae) (Tang et al., 2014), Cincticostella fusca (Ephemerellidae) (Li et al., 2020) and species of Serratella (Ephemerellidae) (Xu et al., 2020b), but three copies of the tRNA-Ile inversion in Ephemerellidae sp. as well as species of Torleya (Ephemerellidae) (Xu et al., 2020a; Xu et al., 2020b) and two copies of the tRNA-Ile in Uracanthella sp. were formed through different random duplications of the tRNA-Ile gene inversion.

Phylogenetic analyses and the Chiastomyaria hypothesis

Each of the three hypotheses, Palaeoptera, Chiastomyaria or Metapterygota, have been supported in different studies using molecular datasets and various analysis methods (see Table S1). In our results, the use of increased numbers of mayfly, dragonfly, damselfly, and stonefly mt genomes sheds light on the phylogenetic origins of winged insects. Although outgroup choice may affect the phylogenetic relationships of Ephemeroptera and Odonata, we recovered support for the Chiastomyaria hypothesis with the choice of different outgroups. The hypothesis of Odonata as the origin of winged insects was also supported by Kjer (2004), Wan et al. (2012), Misof et al. (2014), Simon et al. (2012), Li, Qin & Zhou (2014). But some researchers supported the Metapterygota hypothesis using the mitochondrial genomes method (Tang et al., 2014; Cai et al., 2018; Zhang, Song & Zhou, 2008a; Zhang et al., 2008b).

Suitable outgroups are very important to determining phylogenetic relationships in the rapid evolution of insects. Hence, we must carefully choose the appropriate outgroups in order to discuss the relationships of Ephemeroptera and Odonata. According to the three outgroups tested in the current study, we propose that four bristletails and three silverfishes together is an ideal outgroup to use to evaluate the phylogenetic relationship of Ephemeroptera and Odonata. When we chose four bristletails and three silverfishes as the outgroups (113 dataset), the Chiastomyaria hypothesis was supported in ML and BI analyses (Fig. 3). When we chose five bristletails and three silverfishes as the outgroups (114 dataset), the Chiastomyaria hypothesis was only supported in BI analysis (Fig. S1), whereas the phylogenetic relationship of Ephemeroptera and Odonata failed to recover and formed an different relationship in PhyML analysis (Fig. S3). However, Cai et al. (2018) supported the Metapterygota hypothesis using the same outgroups (five bristletails and three silverfishes). When we used five diplurans as outgroups we also found strong support for the Chiastomyaria hypothesis in BI analysis (Fig. S2), whereas the phylogenetic relationship of Ephemeroptera and Odonata failed to recover in PhyML analysis (Fig. S4). Song et al. (2019) supported the Palaeoptera hypothesis using two species of Collembola and five species of Diplura as outgroups. We compared the outgroups used in published papers (Table S1) and found that few studies use the same outgroups. Thomas et al. (2013) thought that using few suitable outgroups or only distant outgroups may cause systematic error even in large datasets. Despite testing three outgroup combinations and finding support for Chiastomyaria hypothesis, the possibility still exists of bias in outgroup choice.

The phylogenetic family-level relationships in Odonata and Ephemeroptera

The phylogenetic relationships revealed among Odonata families in the current study are similar to the results of Yong et al. (2016) and Dijkstra et al. (2014). We demonstrated that the mt genome was a suitable marker to discuss the family-level phylogenetic relationships of Odonata, as also reported by Yong et al. (2016). However, the monophyly of Coenagrionidae was not supported in our study.

In the phylogenetic relationships at the family level, Siphluriscus chinensis (Siphluriscidae) was supported as sister to the rest of the Ephemeroptera as in the previous results of Ogden et al. (2009), Li, Qin & Zhou (2014), and Zhou & Peters (2003). Isonychiidae was a sister group to Heptageniidae, making the Suborder Setisura (=Heptagenioidea) as supported by McCafferty (1991), whereas Ogden & Whiting (2005), Ogden et al. (2009), Gao et al. (2018), Ye et al. (2018), Cao et al. (2020), and Xu et al. (2020a) supported Isonychiidae as sister to most Ephemeroptera (excluding Baetidae and Siphluriscidae).

The monophyly of Ephemerellidae with different copies of the tRNA-Ile inversion is well supported in this study. Our study identifies Ephemerella sp. as the sister group to the clade of Drunella sp. and Serratella sp. with one copy of the tRNA-Ile inversion. Ephemerellidae sp. with three copies of the tRNA-Ile inversion is the sister group to Uracanthella sp. with two copies. Ogden et al. (2009) found that Ephemerella and Serratella were not supported as monophyletic. In the present study we found a tRNA-Ile translocation and inversion in Ephemerellidae and a different copy of the tRNA-Ile inversion in Ephemerella and Serratella. This gives us more evidence to identify the monophyly of Ephemerella and Serratella in future studies. However, the mt genomes of Ephemeroptera are suitable markers to discuss the family-level phylogenetic relationship of Ephemeroptera.

The phylogenetic relationship of Plecoptera and Dermaptera

In previous studies of Neoptera, one Plecoptera (Pteronarcys princeps) and one Dermaptera (Challia fletcheri) species were included (Li, Qin & Zhou, 2014). Results indicated that P. princeps was sister to Orthoptera whereas Challia fletcheri with a very long-branch was sister to either the rest of Polyneoptera or Ephemeroptera. When only one species of Plecoptera (P. princeps) was included, but Dermaptera (C. fletcheri) was excluded (Lin, Chen & Huang, 2010; Zhang et al., 2010), Plecoptera was sister to Ephemeroptera. However, when we increased the sampling of Plecoptera and Dermaptera, we found that Plecoptera was sister to the remaining Neoptera, and Dermaptera was sister to Neoptera excluding Plecoptera. The position of Plecoptera as sister to the rest of Neoptera and not within Polyneoptera was very interesting. Although Beutel & Gorb (2006) proposed Plecoptera to be in this position, Cai et al. (2018) and Song et al. (2019) supported the position of Plecoptera as sister to the rest of Neoptera in the phylogenetic trees of BI and ML analyses. Some researchers also supported Plecoptera as likely occupying a position near the root of the Neoptera or Polyneoptera, and the aquatic life history stage to be an ancestral feature of winged insects (Marden & Kramer, 1994; Marden, 2013; Zwick, 2009). Most other studies placed Plecoptera within Polyneoptera, often as sister to Dermaptera (e.g., Kjer, 2004; Misof et al., 2007; Simon et al., 2012; Von Reumont et al., 2009). However, Matsumura et al. (2015) supported Plecoptera as the sister group to Zoraptera. Song et al. (2016) reported that Plecoptera was sister to Dermaptera but their ingroup data included only Polyneoptera (no Holometabola and Paraneoptera). Wipfler et al. (2019) found Plecoptera was sister to Polyneoptera excluding Zoraptera and Demaptera and suggested the ancestor of winged insects did not evolve in an aquatic environment. In our study, we recovered the relationship Odonata + (Ephemeroptera + (Plecoptera + (Dermaptera + other Neoptera))) with high posterior probabilities and bootstrap values (Fig. 3). So, our results suggest that the common ancestors of the winged insects and of Neoptera were aquatic, supporting the idea that wings did evolve in an aquatic environment.

The monophyly of Polyneoptera was not supported in this study because of the placement of Plecoptera, Dermaptera and Orthoptera. Although recent studies have supported the monophyly of Polyneoptera using mt genomes or transcriptome data (Song et al., 2016; Wipfler & Pass, 2014; Wipfler et al., 2019), we found most of studies included just the Polyneopteran orders. It is suggested that such analyses should always add more species representing all insect orders, not only Polyneoptera but also including Holometabola in order to discuss the monophyly of Polyneoptera. Our analysis of the 113 dataset indicated that Orthoptera was sister to Amphiesmenoptera (Lepidoptera + Trichoptera) which may be caused by long-branch attraction of Hemiptera (Fig. 3). However, we failed to support the monophyly of Orthoptera in BI analysis of the 119 dataset.

In conclusion, we found that the mt genome is a suitable marker to investigate the phylogenetic relationship of the inter-order and inter-family relationships of insects but that the outgroup choice is very important for deriving phylogenetic relationships among winged insects. We highly recommend that we should choose suitable species from Archaeognatha and Zygentoma together as the outgroups in future research and discussions of the phylogenies of Ephemeroptera and Odonata.