Adult specimens of M. tritaenia and M. metafiligera were collected from Baiyun Mountain (Luoyang, Henan Province, China) in July 2015 and from Bolonggong (Tibet, China) in August 2015, respectively. We preserved specimens in 95% ethanol in the field and stored them at – 20 °C until tissues were used for DNA extraction. Voucher specimens of M. tritaenia (No. Vhem–0008) and M. metafiligera (No. Vhem–0061) were deposited in Entomological Museum of Henan institute of Science and Technology (HIST), Henan Province, China. Specimens were identified by Wei-Hai Li (HIST). We extracted total genomic DNA from thoracic muscle tissue using QIAamp DNA Blood Mini Kit (Qiagen, Duesseldorf, Germany) according to the manufacturer’s protocols.

The mitogenomes were amplified and sequenced as described in previous studies (Song et al. 2016; Li et al. 2017; Wang et al. 2017a, 2017b; Wang et al. 2018a, 2018b). From the genomic DNA, an Illumina TruSeq library with an insert size of 480 bp was generated. We sequenced our library on a single lane of Illumina Hiseq 2500 with 500 cycles of paired-end sequencing (250 bp reads). High-quality reads were generated using Trimmomatic v0.30 (Lohse et al. 2012). Then, IDBA-UD (Peng et al. 2012) was used for de novo assembling with these reads. A similarity threshold of 98% and minimum and maximum k values of 80 and 240 bp, respectively, were used to build assemblies. COI and srRNA fragments were amplified as bait sequences using PCR (Li et al. 2012; Wang et al. 2014) and previously designed primers (Simon et al. 2006) were used to determine the mt genome components. Raw sequences from the mitochondrial genome of each species were assembled into contigs with BioEdit 7.0.5.3 (Hall 1999). We identified protein-coding genes and two ribosomal RNA genes using BLAST searches in NCBI and by alignment with homologous genes from published stonefly species. tRNA genes were identified using the ARWEN online service and checked manually (Laslett and Canbäck 2008). Each protein-coding gene was aligned individually based on codon-based multiple alignments using the MAFFT algorithm within the TranslatorX online platform (Abascal et al. 2010). Before back-translating to nucleotides, poorly aligned sites were removed from the protein alignment using Gblocks in the TranslatorX with default settings. Nucleotide substitution rates, base composition, and codon usage were analyzed with MEGA 6.0 (Tamura et al. 2013). The GC and AT skews were calculated using the formulae: AT skew= (A−T)/ (A+T) and GC skew= (G−C)/ (G+C) (Perna and Kocher 1995).

We calculated the value of Ka (the rates of non-synonymous substitutions), Ks (the rates of synonymous substitutions) using DnaSP 5.1.0 (Librado and Rozas 2009). We also identified tandem repeats in the A+T-rich region using the Tandem Repeats Finder server (Benson 1999), and predicted the stem-loop structure using the mfold web server (Zuker 2003).

The phylogenetic trees among three families of the superfamily Nemouroidea were reconstructed using DNA data from nine published and two newly sequenced mitogenomes. Two stonefly species, Pteronarcys princeps and Styloperla spinicercia, were used as outgroups (Table 1). We assembled the following datasets for phylogenetic analyses: 1) the “PCGs matrix” (total of 11,046 bp), inclusive of 13 PCGs; 2) the “PCG12 matrix” (total of 7,364 bp) which contains the first and second codon positions of protein-coding genes; 3) the “PCGR matrix” (total of 13,056 bp) which contains all three codon positions of protein-coding genes and two ribosomal RNA genes; 4) the “PCG12R matrix” (9374 bp) which contains the first and second codon of PCGs and two ribosomal RNA genes. The BI and ML analyses were carried out on the CIPRES Science Gateway (Miller et al. 2010). ML analyses were performed with RaxML-HPC2 on XSEDE 8.1.10 (Stamatakis 2006) using the GTRGAMMAI model, and node confidence was assessed with 1,000 bootstrap replicates. Bayesian analyses were carried out using MrBayes 3.2.6 on XSEDE (Ronquist et al. 2012). For BI analyses, GTR+I+G was the best-fit model for the nucleotide sequence alignments, using jModelTest 0.1.1 based on Akaike’s information criterion (AIC) (Posada 2008). We conducted with two simultaneous runs for 10 million generations. Each set was sampled every 1,000 generations with a burn-in rate of 25%. Stationarity was examined by Tracer v.1.6 (http://tree.bio.ed.ac.uk/software/tracer/), and was considered to be reached when the ESS (estimated sample size) value was above 200.

In the present study, two complete mitogenomes of M. metafiligera and M. tritaenia were sequenced and deposited in GenBank of NCBI under accession numbers MH085450–MH085451 (see Table 1 and Fig. 1). The M. metafiligera mitogenome was a 15,739 bp circular DNA molecule, which was smaller than M. tritaenia (15,778 bp) due to differences in the size of the A+T-rich region (Tables 1 and 2). The size of completely sequenced mitogenomes were medium-sized when compared with the mitogenomes of other stoneflies, which ranged from 16,602 bp (N. nankinensis) (Chen and Du 2017a) to 15,048 bp (Isoperla bilineata) (Chen et al. 2018). All genes identified in the two mitogenomes are typical animal mitochondrial genes with normal gene sizes. Mitochondrial gene order was the same as all previously published stonefly mitogenomes, as well as the ancestral gene order of insects, as exemplified by Drosophila yakuba (Clary and Wolstenholme 1985). When compared with the mitogenome of N. nankinensis (Chen and Du 2017a), the length variation was limited in the PCGs, tRNA, and rRNA genes, but very different in the A+T-rich region.

Figure 1. 

Map of the mitogenomes of M. metafiligera and M. tritaenia. Genes shown on the inside of the map are transcribed in a clockwise direction, whereas those on the outside of the map are transcribed counterclockwise. Different gene types are shown as filled boxes in different colors. Numbers show the sizes of intergenic spacers (positive values) and overlapping region between genes (negative values).

In the mitogenome of M. metafiligera, 50 overlapping nucleotides were located in 15 pairs of neighboring genes, while in the mitogenome of M. tritaenia, there were 51 overlapping nucleotides in 16 gene boundaries. These overlapping nucleotides varied in length from 1 to 8 bp (Table 2). When compared with the mitogenome of N. nankinensis, three conserved regions were found in overlapping regions of genes of each sequenced nemourid species and were also observed in some other stoneflies: AAGCCTTA (tRNA^Trp-tRNA^Cys), ATGATTA (ATP8–ATP6) and ATGTTAA (ND4–ND4L) (James and Andrew 2006; Qian et al. 2014; Wu et al. 2014; Elbrecht et al. 2015; Huang et al. 2015; Wang et al. 2017a, 2017b, 2018a). The ATP8/ATP6 and the ND4/ND4L gene pairs overlap are also found in many insect mitogenomes and thought to be translated as a bicistron (Stewart and Beckenbach 2005). The longest intergenic spacer we found located between tRNA^Ser(UCN) and ND1 in two nemourid species, ranging from 35 bp in M. metafiligera to 38 bp in M. tritaenia (Table 2). In Drosophila melanogaster, two conserved non-coding intergenic regions (tRNA^Glu-tRNA^Phe, and tRNA^Ser(UCN)–ND1) have been shown to be binding sites for a bidirectional transcription termination factor (DmTTF) (Beckenbach 2012). For the first binding site of tRNA^Glu-tRNA^Phe, this region was absent in two nemourid mitogenomes (Table 2). Studies show that this region is absent or incomplete in other plecopteran species as well as in other insect orders (Beckenbach 2012). However, the second binding site of tRNA^Ser(UCN)–ND1 is found in two nemourid species. This region is more widely conserved and similar non-coding sequences are present at this site in other insect orders (Beckenbach 2012).

Similarly to what observed in other insects, the nucleotide composition of two nemourid mitogenomes was clearly biased towards a higher content of A/T nucleotides. The base composition bias of the M. metafiligera mitogenome is 69.1% A + T content, made up of 66.9% in the PCGs, 71.3% in the tRNAs, 72.6% in the rRNAs and 84.1% in the A + T rich region. While the base composition bias of the M. tritaenia mitogenome is 68.6% A + T content, made up of 66.4% in the PCGs, 71.3% in the tRNAs, 72.0% in the rRNAs and 82.1% in the A + T rich region (Table 3). This is identical to the base composition biases reported in other stoneflies, which ranges from 62.5% (Dinocras cephalotes Curtis, 1827) to 71.5% (Pteronarcys princeps Banks, 1907) (James and Andrew 2006; Elbrecht et al. 2015). After comparing the nucleotide compositions of four major partitions (PCGs, tRNAs, rRNAs, and the A+T-rich region), we found that all four partitions were consistently biased towards A and T. Among the four partitions, the lowest A + T content was found in PCGs, whereas the highest in the A+T-rich region of the two nemourid specimens. The strand bias also can be measured as AT– and GC-skews (Perna and Kocher 1995). Both two nemourid species had positive AT-skew and negative GC-skew values for the entire mitogenome (Table 3), showing a biased use for the A and C nucleotides, which was common in most other insects (Wei et al. 2010).

The total length of all PCGs of both mitogenomes of M. metafiligera and M. tritaenia were 11,199 bp (Table 3). Most PCGs of the two mitogenomes initiated with the typical start codon ATN (Met/Ile), however TTG was proposed as the start codon for ND1 in two species (Table 2). Gene ND1 has been found to employ TTG as a start codon in some insects, thus minimizing intergenic spacing and avoiding overlap with adjacent genes (Bae et al. 2004; Sheffield et al. 2008; Wei et al. 2010). The ND5 gene of M. metafiligera has the unusual start codon GTG, as previously reported in other stoneflies (James and Andrew 2006; Elbrecht et al. 2015; Huang et al. 2015; Sproul et al. 2015; Chen et al. 2016; Elbrecht and Leese 2016; Chen and Du 2017b; Wang et al. 2017a, 2017b).

Similar to most other stoneflies, the most commonly used stop codon in two nemourid specimens was TAA, which was found in eight of the 13 PCGs (ATP6, ATP8, COI, COIII, ND1, ND4, ND4L and ND6) for both two nemourid mitogenomes. The stop codon TAG was used in ND2, ND3 and CytB. In M. metafiligera and M. tritaenia, both COII and ND5 terminate with incomplete stop codon T (Table 2). The phenomenon of incomplete stop codons is common in insect mitogenomes and it is likely that these truncated stop codons are completed by posttranscriptional polyadenylation (Ojala et al. 1981). Overall, the start and stop codons were similar in two nemourid specimens, and showed little differences between them in some PCGs (Table 2).

The genome-wide bias toward A+T content was also reflected in the codon usage by the PCGs. The relative synonymous codon usage (RSCU) showed high base compositional biases for AT in the mitogenomes of M. tritaenia and M. metafiligera (Suppl. material 1, Figure S1). The codons ending with A or U were preferred in both the four-fold and two-fold degenerate codons. However, in M. metafiligera mitogenomes, the RSCU value of GGG (Gly) was higher than GGT (Gly) (Suppl. material 1, Figure S1). This may be caused by the anticodon tRNA^Gly-TCC, which contain swinging nucleotides and can decode GGG (Mao and Dowton 2014).

In order to analyze the evolutionary patterns of PCGs in the three nemourid mitogenomes, we compared the ratio of Ka/Ks for each PCGs (Suppl. material 2, Figure S2). The observed average Ka/Ks ratios were consistently lower than 0.3, increasing from 0.025 for ATP6 to 0.202 for ND6 (Suppl. material 2, Figure S2). The uniformly low values of Ka/Ks ratios for COI, COII, COIII and CytB indicate strong evolution constraints in cytochrome c oxidase (Schmidt et al. 2001; Zsurka et al. 2010). The Ka/Ks values for all PCGs were below 0.3, indicating that these genes are evolving under purifying selection. Therefore, all mitochondrial PCGs could be used to investigate phylogenetic relationships within nemourid.

The sizes of 22 tRNA genes of M. metafiligera and M. tritaenia range from 63 bp to 71 bp, comprising 9.36% (1,473 bp) and 9.32% (1,471 bp) of the total mitogenomes, respectively (Table 2 and 3). By compared with N. nankinensis mitogenome, all 22 tRNA genes in the three nemourid species were identified and the secondary structures were shown in Fig. 2. Most tRNA genes could be folded into a classic clover-leaf secondary structure except for tRNA^Ser(AGN) due to a lack of the dihydrouridine (DHU) arm (Fig. 2). The loss of the DHU arm in tRNA^Ser(AGN) was a typical feature in metazoan mitogenomes (Negrisolo et al. 2011). However, we found that tRNA^Ser(AGN) in the three nemourid mitogenomes possessed an unusual anticodon stem (9 pairs of nucleotides) with an extended nucleotide (Fig. 2). This is an unusual phenomenon, but has also been observed in the DHU arm in other stoneflies (Chen and Du 2017b; Wang et al. 2017a; Wang et al. 2018a, 2018b).

Figure 2. 

Secondary structure of tRNA families in nemourid mitogenomes. The nucleotide substitution pattern for each tRNA family is modeled using as reference the structure determined for M. tritaenia. Red arrows correspond to insertions. Inferred Watson-Crick bonds are illustrated by lines, whereas GU bonds are illustrated by dots.

According to the secondary structures and sequence alignment, the most conserved tRNAs in nemourid mitogenomes were tRNA^Val, tRNA^Thr, tRNA^Ser(AGN) and tRNA^Leu(CUN) with no more than two nucleotides substitution in each gene (Fig. 2). Nucleotide insertion-deletion polymorphisms were mainly restricted to TΨC and DHU loops. As expected, tRNAs of M. tritaenia and M. metafiligera showed high level of sequence and structural identity, with five identical tRNAs (tRNA^Val, tRNA^Ser(UCN), tRNA^Ser(AGN), tRNA^Leu(UUR) and tRNA^Phe) (Fig. 2). The eighteen remaining tRNAs only had four insertion-deletion positions in total (one insertion position in the DHU loop of tRNA^His, each with one insertion position in the TΨC loop of tRNA^Glu and tRNA^Tyr, and one deletion position in the TΨC loop of tRNA^Ile), and the nucleotide substitutions of individual gene were mostly restricted to 1–3 nucleotide sites, with the exception of tRNA^His and tRNA^Met with 5 and 6 sites respectively (Fig. 2).

The large subunit ribosome gene (lrRNA) was located between tRNA^Leu(CUN) and tRNA^Val, while the small subunit ribosome gene (srRNA) was flanked by tRNA^Val and the A+T-rich region (Fig. 1). Two ribosome genes (srRNA and lrRNA) were identified at 794 and 1,329 bp in M. tritaenia and at 792 and 1,327 bp in M. metafiligera (Table 2). Our result shows that the size differences in both ribosomal subunits are not distinct among different species.

In this study, the secondary structures of the srRNA and lrRNA of M. tritaenia were constructed following the models proposed for other insects, with marked the conserved sites. The secondary structure of lrRNA contained five structural domains (I–II, IV–VI, domain III is absent in arthropods) and 44 helices (Fig. 3). The multiple alignments of nemourid lrRNA genes had 1,334 positions and contained 1,139 conserved positions (85.4%), 190 necleotide substitutions (14.2%) and 5 indels (0.4%), respectively. Nucleotide variability was unevenly distributed among domains and helices, mainly in domains I and II. Several helices (H821, H837, H946, H991, and H1196) in lrRNA, with high variability at the primary sequence level, showed conserved secondary structures. The secondary structure of srRNA contained three domains (I–IV) and 27 helices (Fig. 4). The multiple alignments of nemourid srRNA genes was inclusive of 795 positions and contained 737 conserved positions (92.7%), 56 necleotide substitutions (7.0%) and 2 indels (0.3%), respectively. Nucleotide variability was unevenly distributed among domains and helices, mainly in helices H47 of domain I.

Figure 3. 

Predicted secondary structure of the lrRNA in M. tritaenia. Key: red circle, nucleotide conserved in three nemourid mitogenomes; green circle, nucleotide not conserved. Roman numerals represent the conserved domain structures. Dashes (−) indicate Watson-Crick base pairing, and dots (·) indicate G–U base pairing I–VI indicate six domains in the secondary structure of lrRNA.

Figure 4. 

Predicted secondary structure of the srRNA in M. tritaenia. Key: red circle, nucleotide conserved in three nemourid mitogenomes; green circle, nucleotide not conserved. Roman numerals represent the conserved domain structures. Dashes (−) indicate Watson–Crick base pairing, and dots (·) indicate G–U base pairing I–VI indicate six domains in the secondary structure of srRNA.

The A+T-rich regions all locate between srRNA and tRNA^Ile in M. tritaenia and M. metafiligera mitogenome with 921 bp and 893 bp in length, respectively (Table 2). This region is considered to be the control region as it contains both an origin of transcription and replication (Taanman 1999; Cameron 2014). The A+T-rich region is the most variable region in mitogenome due to insertion and deletion events, variation in numbers of tandem repeats, and differences in variable domain length (Zhang et al. 1995; Taanman 1999). In addition, some structural elements are founded in the A+T-rich region, such as poly-T stretch, stem-loop (SL) structures, tandem repeats (TRs) and G+A-rich region, etc. (Zhang and Hewitt 1997). A comparison of molecular features in A+T-rich region between the two newly sequenced mitogenomes was shown in Fig. 5. Some essential elements were observed: (1) a leading sequence adjacent to srRNA; (2) large tandem repeats present as two or more copies; (3) the remainder of the A+T-rich region (Fig. 5A).

Figure 5. 

A+T-rich region of two nemourid mitogenomes A structure elements found in the A+T-rich region of two nemourid mitogenomes. TR is the abbreviation of tandem repeat units B putative stem-loop structures found in the A+T-rich region of M. metafiligera (MM indicate M. metafiligera) C putative stem-loop structures found in the A+T-rich region of M. tritaenia (MT indicate M. tritaenia).

Repeated sequences are common in the A+T-rich region for most insects, and length variations are decided to a large degree by the various numbers of repeats (Zhang and Hewitt 1997; Zhang et al. 1995). In the case of M. tritaenia, there is a large repeat region (positions: 15,243–15,483) which is 241 bp and contains two tandem repeat units plus a partial copy of the repeat. While in the A+T-rich region of M. tritaenia mitogenome, two tandem repeat regions (TR1 and TR2) were founded. TR1 region (positions: 15,265–15,428) is 164 bp and contains three tandem repeat units plus a partial copy of the repeat; In TR2 region (positions: 15,623–15,661), the two tandem repeats units were 21 and 18 bp (not strictly repeats). We also found some SL structures in two newly sequenced species. Two SL structures were predicted in the A+T-rich region of M. metafiligera mitogenome: SL–1 (positions: 15,122–15,146) and SL–2 (positions: 15,193–15,208) (Fig. 5B). In the A+T-rich region of M. tritaenia mitogenome, we found two long repeated sequences: 5’–TATATATATATATATAAATTAATAAATAATTTATTGTA–3’. Both of them can be folded into a stem-loop structure (SL–1 positions: 15,306–15,343; SL–2 positions: 15,352–15,389) (Fig. 5C). This special structure in the A+T-rich region of M. tritaenia mitogenome is quite different with that of other published stoneflies except for Acroneuria hainana Wu, 1938, and it may be able to adjust the replication speed of two replicate directions (Huang et al. 2015). The proposed “G(A)_nT” motif was detected in SL structures of M. metafiligera mitogenome, while it was modified as “GTA” in SL structures of M. tritaenia mitogenome. In the A+T-rich region of N. nankinensis mitogenome, five SL structures with three different motifs (“GTA”, “G(A)_nT” and “TGA”) were detected (Chen and Du 2017a). The stem-loop structure in the A+T-rich region is identified in many insects and it is thought to be the site of the initiation of minority strand (N-strand) synthesis in Drosophila (Clary and Wolstenholme 1987). By contrast, some tRNA-like structures also founded in N. nankinensis mitogenome, but these structures were not reported in our newly sequenced species.