Mitochondrial Genome of the Stonefly Kamimuria wangi (Plecoptera: Perlidae) and Phylogenetic Position of Plecoptera Based on Mitogenomes

This study determined the mitochondrial genome sequence of the stonefly, Kamimuria wangi. In order to investigate the relatedness of stonefly to other members of Neoptera, a phylogenetic analysis was undertaken based on 13 protein-coding genes of mitochondrial genomes in 13 representative insects. The mitochondrial genome of the stonefly is a circular molecule consisting of 16,179 nucleotides and contains the 37 genes typically found in other insects. A 10-bp poly-T stretch was observed in the A+T-rich region of the K. wangi mitochondrial genome. Downstream of the poly-T stretch, two regions were located with potential ability to form stem-loop structures; these were designated stem-loop 1 (positions 15848–15651) and stem-loop 2 (15965–15998). The arrangement of genes and nucleotide composition of the K. wangi mitogenome are similar to those in Pteronarcys princeps, suggesting a conserved genome evolution within the Plecoptera. Phylogenetic analysis using maximum likelihood and Bayesian inference of 13 protein-coding genes supported a novel relationship between the Plecoptera and Ephemeroptera. The results contradict the existence of a monophyletic Plectoptera and Plecoptera as sister taxa to Embiidina, and thus requires further analyses with additional mitogenome sampling at the base of the Neoptera.


Introduction
The order Plecoptera is comprised of 16 families and includes 3,497 species [1]; these occur on all continents except Antarctica [2]. Plecoptera is a small order of hemimetabolous insects that are primarily associated with clean, cool, running water and cool, damp terrestrial environments [3]. The nymphs generally have their highest density in riffle areas of streams where rocks, gravel, snags, and accumulated leaves are abundant. Their potential use as biological indicators of water quality is well-known [4].
The insects mitochondrial genome is a circular molecule ranging from 15 to 20 kp; it generally encodes two rRNA genes, 22 tRNA genes, 13 protein-coding genes (PCGs), and an A+T-rich region varying in length among taxa [5,6]. Recently, whole mitochondrial DNA (mtDNA) sequences have been used in phylogenetic analyses of various insect orders. After the complete mitochondrial genome of Pteronarcys princeps was published [7], increasing numbers of researchers used this data for phylogenetic analyses, which resulted in various theories on the phylogenetic position of Plecoptera. Zhang et al. and Lin et al. [8,9] used 13 PCGs of mtDNA data to support relatedness between Ephemeroptera and Plecoptera within the Neoptera.
The stoneflies are an ancient group of insects. Their fossil record extends into the Lower Permian [10]. But the phylogenetic position of Plecoptera has long been under debate. Some of the competing evolutionary hypotheses place the stoneflies in a single basal group, diverging early after the split between the Paleoptera and Neopteran ancestors, and before other Neopteran diversification [11,12]. At the same time, the relationship about Styloperlidae, Peltoperlidae, Chloroperlidae, Perlidae and Perlodidae was controversial [13,14]. The aim of this study was to discuss the phylogenetic position of Plecoptera based on mtDNA.
In this work, we present the complete mitochondrial genome of the stonefly, Kamimuria wangi Du, 2012 (Perlidae) and reconstructed phylogenies of 17 major basal Pterygote lineages based on 13 PCGs obtained from mitochondrial genomes.

Genome organization, structure and composition
The full-length mtDNA of K. wangi is a circular molecule of 16,179 nucleotides, just slightly longer than that of P. princeps (16004 bp). The K. wangi mtDNA contains 37 genes (13 PCGs, 22 tRNA genes, two rRNA genes) and an AT-rich control region. The number and arrangement of these genes is consistent with insect mitochondrial DNA [5,15] and identical to the P. princeps mitogenome ( Figure 1). In 12 locations of the mtDNA, genes overlapped by 1-47 bp (excluding the A+T-rich region; see Table 1). The mtDNA of K. wangi also contained intergenic regions (1-28 bp in length) in 13 locations.
Thirteen PCGs in K. wangi utilize the conventional translational start and stop codons for invertebrate mtDNA. For example, nine PCGs (ND2, COII, ATP6, COIII, ND5, ND4, ND4L, ND6 and CytB) contained ATG codons. Three PCGs (COI, ATP8 and ND1) initiated with ATA codons, and one PCG (ND3) contained ATT as the start site. However, in P. princeps, ND2 and ND5 initiated with GTG, ND6 and ATP8 with ATT, and ND1 gene initiated with TTG. Thirteen PCGs used the typical termination codons TAA and TAG in K. wangi. In contrast, only two PCGs encode conventional termination codons in P. princeps, e.g. TAA for ND4 and TAG for ND1.
The relative synonymous codon usage (RSCU) values showed a biased use of A and T nucleotides in K. wangi ( Table 2). The high incidence of anticodons NNA and NNU indicated partiality for AT in the anticodon of PCGs. Ile (8.12%), Leu (11.61%), Ser (11.99%) and Thr (7.03%) were the most frequent amino acids in K. wangi mtDNA PCGs, and Ser showed the highest incidence. In P. princeps PCGs, Ile (7.51%), Leu (12.23%), Asn (7.48%) and Ser (12.15%) were the most common amino acids, and Leu occurred most frequently.
The overall AT content in K. wangi mtDNA was 69.6%, which was slightly lower than that of P. princeps (70.6%). The average AT content across all PCGs in K. wangi was 68.0%, which also lower than that of P. princeps (70.5%). The nucleotide composition of the insect mitogenome was biased toward A and T nucleotides. The AT content of the third codon (73.3%) was much higher than the first (63.5%) and second codon positions (67.1%); this suggested that both higher mutation rates and the increased AT occurrence are related and dependent on relaxed selection at the third codon position. However, the A+T% at the second codon (65.7%) was lower than the first (71.0%) and third codon positions (74.9%) in P. princeps. Furthermore, we present the detailed comparison in the nucleotide composition, AT-skew and GC-skew between the two stonefly species, K. wangi and P. princeps ( Table 3).
The mitogenome of K. wangi had 22 traditional tRNA genes interspersed with rRNAs or PCGs. The tRNA genes were encoded by a total of 1484 nucleotides; tRNAs ranged from 64 to 71 bp and the A+T content was 72.4%. Fourteen tRNAs were encoded by the H-strand and the remaining eight were encoded by the Lstrand. All tRNA genes had the typical cloverleaf secondary structure except for the tRNA Ser(AGN) gene; in this molecule, the stable stem-loop structure of the dihydrouridine (DHU) arm was missing ( Figure 2), a feature that has been observed in other metazoan mitogenomes [16,17].  There were two rRNAs in K. wangi with a combined length of 2221 bp and an A+T content of 71.0%. The lrRNA gene (rrnL) had a length of 1345 bp, an A+T content of 73.1%, and was positioned between tRNA Leu(CUN) and tRNA Val . The srRNA gene (rrnS) was 867 bp, had an A+T content of 67.6%, and was located between the tRNA Val and A+T-rich region. Gene sizes and map positions were consistent with those observed in the mitogenomes of other insects. Furthermore, both the AT-and GC-skews were slightly positive in the tRNA and rRNA genes, which correlates with the results obtained for P. princeps. The secondary structures of lrRNA and srRNA were consistent with the models proposed for these genes in other insects [18][19][20]. In K. wangi, the lrRNA gene contained six domains (labeled I, II, III, IV, V and VI) with 44 helices (Figure 3). The srRNA gene contained three domains (labeled I, II, III) and 32 helices ( Figure 4).
Previously, the A+T-rich region was reported to contain elements essential to the initiation of replication and transcription [21]. The A+T-rich region of the K. wangi mitogenome was 1251 bp and mapped between the srRNA and tRNA Ile -tRNA Gln -tRNA Met gene cluster ( Figure 1). The A+T content of the AT-rich region was 78.2%, slightly lower than the corresponding region in P. princeps (81.3%). Additionally, there was a 10-bp poly-T stretch observed in the A+T-rich region of K. wangi. Two regions located downstream of the poly-T stretches were identified and designated stem-loop 1 (positions 15848-15651) and stem-loop 2 (positions 15965-15998) based on their potential ability to fold into stemloop structures ( Figure 5).

Phylogeny of basal Neoptera
In this study, the amino acid sequences of the 13 PCGs were concatenated to construct phylogenetic relationships, which may result in a more complete analysis than analyzing each sequence separately. We incorporated species from Orthoptera, Ephemeroptera, Blattodea, Mantodea, Mantophasmatodea, Phasmatodea, Lepidoptera, Embioptera, Hemiptera, and Hymenoptera in the phylogenetic analysis (Table 4). Zoraptera and Dermaptera were not included in the analyses because mtDNA data are incomplete for these groups. Odonata was included as an outgroup.
Analysis using BI and ML resulted in two trees with the same topology except for some variation in node confidence values ( Figure 6). Numbers at each node indicate bootstrap support, percentages of Bayesian posterior probabilities (first value) and ML bootstrap support values (second value).
The results support the Embioptera was a sister group of a clade containing Hemiptera and Hymenoptera. S. immanis (Ephemeroptera) grouped with the mitogenomes of P. princeps and K. wangi (Plecoptera) in Polyneoptera; this grouping differed from conclusions based on morphological classifications and some molecular studies, but was consistent with the findings of Zhang et al. [8] and Lin et al. [9].
Plecoptera are traditionally associated with the lower Neoptera and often placed within Polyneoptera. There has been very little consensus regarding the placement of Plecoptera within this group, which is comprised of ten other orders: Blattodea, Dermaptera, Embiidina, Grylloblattodea, Isoptera, Mantodea,   Mantophasmatodea, Orthoptera, Phasmatodea and Zoraptera [22][23][24]. Boudreaux placed Plecoptera within Polyneoptera as the sister taxon to Embiidina [11]. Henning recognized all polyneopterous orders except Plecoptera as a monophylectic group and assigned it to the Paurometabola [12]. However, he was unable to assign Plecoptera and left it as a lineage disconnected from the overall topology. Kukalova-Peck depicted Polyneoptera as paraphyletic and placed Plecoptera as the sister taxon to the ''Orthopteroid orders'' [25]. However, Kristensen leaves Polyneoptera (including Plecoptera) as a largely unresolved polytomy at the base of Neoptera [22]. Molecular studies generally place Plecoptera as a sister to Dermaptera based on the gene sequences encoding the small subunit of nuclear ribosomal DNA [26]. Wheeler et al. used both molecular and morphological data to support a monophyletic Polyneoptera and placed Plecoptera as sister taxon to Embiidina [24]. However, Terry & Whiting [27] used sequences derived from 18S rDNA, 28S rDNA and histone 3 to assign Plecoptera as a sister lineage to Dermaptera and Zoraptera. The evolution of hemocyanin subunits follows the widely-accepted phylogeny of the Hexapoda and provides strong   evidence for the monophyly of the Polyneoptera (Plecoptera, Dermaptera, Orthoptera, Phasmatodea, Mantodea, Isoptera, Blattaria), which subsequently places Plecoptera as a sister taxon to Embiidina and Zoraptera [28]. The resulting trees strongly support monophyly of Neoptera and Polyneoptera and place Plecoptera as the sister taxon to Dermaptera based on three nuclear protein-coding genes [29]. However, our phylogenetic analyses show an unexpected monophyletic Plecoptera and Ephemeroptera clade sister in Polyneoptera. This result is so interesting, because the two taxa share aquatic life history of naiads and a few studies have used stoneflies as models for understanding the origin of flight in early winged insects [30].
The phylogenetic clustering of the Plecoptera with Ephemeroptera was unexpected, may be a result of taxon sampling at the base of the pterygotes. Broad and sufficient taxon sampling is an important factor in phylogenetic analyses of insect mitochondrial genomes. Accordingly, it is necessary to increase the number of representative species within the target taxa, excluding the highly divergent genomes with variable genome rearrangements, elevated substitution rates, or base compositional bias, which likely cause long branch attraction. Our phylogenetic trees cannot illustrate phylogenetic system of Neopteran insects completely and clearly, because the missing mtDNA sequences from Dermaptera and Zoraptera.

Ethics statement of Animals
There were no special permits for the insect collection for this study in Henan Province, China. The insect samples used in this experiment were collected from a clean creek flowing from the mountain. The field studies did not involve endangered or protected species. The species in the genus of Kamimuria are common small aquatic insects and are not included in the ''List of Protected Animals in China''.

Sample preparation and DNA extraction
Specimens of K. wangi were collected at Luoyang in Henan Province (33.70u N, 111.84uE) in July 2009, China. All specimens were identified and then preserved in 100% ethanol and stored at 270uC until DNA extraction was performed. Genomic DNA was extracted from individual audult samples using the Axygen DNA kit (Axygen Biotechnology Hangzhou, China) and used as a template for LA-PCR.

PCR amplification and sequencing
Two pairs of LA-PCR primers were used to amplify two overlapping portions of mtDNA in K. wangi. The LA-PCR primer sequences were as follows: A1: 59-CGAGCWTACTTTACTT-CAGCMACWATAATTA-39and A2: 59-GCAAATAR RAA-TRATCATTCATTCTGGTTGGAT-39; B1: 59-TACACCAAA-CAGGATCAAAT AAYCCMWTAGG-39and B2: 59-GCTC-CAACATRTTTCTRCATTGACCAAATA-39. The LA-PCR amplifications were performed with LA Taq DNA polymerase (Takara, Japan) in an ABI thermal cycler using the following program: initial denaturation at 93uC for 2 min, followed by 40 cycles at 92uC for 10 s, annealing at 54uC for 30 s, elongation at  Additional PCR primers (shown in Table S1) were designed by comparing mtDNA sequences from 30 Neopteran insect species (including Plecoptera) with universal primers of insect mtDNA [5]. These amplifications were performed using templates from LA-PCR as follows: primary denaturation at 94uC for 2 min, 35 cycles at 94uC for 30 s, annealing at 40-55uC for 30 s, extension at 72uC for 90 s, and final extension at 72uC for 7 min. All PCR fragments from K. wangi were sequenced after separation and purification.
Purified PCR products were ligated into pGEMT Easy Vector (Promega). Recombinant clones were sequenced in both direction using the BigDye Terminator Sequencing Kit (Applied BioSystems) and the ABI 3730XL Genetic Analyzer (PE Applied Biosystems, San Francisco, CA, USA) with two vector-specific primers and internal primers for primer walking. The mitogenome of K. wangi was amplified from overlapping PCR fragments. In order to guarantee the quality of sequences, each PCR fragment contained five copies and the sequences were checked by the Chromas. All well sequenced fragements were assembled using ContigExpress and DNAMAN and then were blasted in NCBI to identify the purpose fragements.

Genome annotation and secondary structure prediction
We used CodonCode Aligner for sequence assembly and annotation. Transfer RNA genes were initially identified using the tRNAscan-SE software available online at http://lowelab.ucsc. edu/tRNAscan-SE/ [31]. PCGs and rRNA genes were identified by comparison with the P. princeps mitochondrial genome. The software of XRNA 1.2.0.b (http://rna.ucsc.edu/rnacenter/xrna/ xrna.html) was used to draw the secondary structure of tRNA genes. The secondary structures of rrnL and rrnS genes were inferred based on models developed for other insect species [32][33][34][35]. To infer the secondary structures of tRNA and rRNA molecules, we used a commonly accepted comparative approach to correct for unusual pairings with RNA-editing mechanisms that are well known in arthropod mitogenomes [36,37]. AT-skew and GC-skew were calculated to determine the AT and GC biases of PCGs in K. wangi and P. princeps, using the following formulas: ATskew = [A2T]/[A+T]; and GC-skew = [G2C]/[G+C] [38]. Translated amino acid sequences of the PCGs were aligned using ClustalW [39] in MEGA 5 [40] after eliminating initiation codons. The AT content and codon usage were calculated using MEGA 5.

Phylogenetic analysis
To examine the phylogenetic position of stonefly, 12 additional Neopteran mitogenomes were downloaded from GenBank and compared using the mitogenome of Odonata (Euphaea formosa) as the outgroup ( Table 4). The amino acid sequences of the 13 PCGs were aligned with ClustalX using default settings and concatenation [41]. The best-fitting model by Modeltest using likelihood ratio tests, was then used to perform Bayesian inferences (BI) and maxi-mum likelihood (ML) analysis using the program MrBayes 3.1.2 (http://morphbank.ebc.uu.SE/mrbayes/) and MEGA version 5.0 software with maximum likelihood (ML) methods [42]. DNA alignments were inferred from the amino acid alignment of 13 PCGs using default settings in ClustalX and MEGA version 5.0, which could alternate between DNA and amino acid sequences within alignments. Alignments of individual genes were then concatenated without the stop codon. The best fit model for nucleotide alignments was determined by Modeltest 3.7. According to the Akaike information criterion, the GTR + I + G paradigm was the most ideal model for analysis using nucleotide alignments. The BI analyses were conducted under the following conditions: 1,000,000 generations, four chains (one cold chain and three hot chains) and a burn-in step for the first 10,000 generations. The confidence values of the BI tree were expressed as the Bayesian posterior probabilities in percentages. The ML methods bootstrap analysis was done with 1000 replications, and values were calculated using the 50% majority rule.