Adult specimens of Haploperla japonica were collected from Yoshino River, Minocho, Miyoshi City, Tokushima Prefecture, Japan, and adult samples of Sweltsa sp. were collected from a forested stream along the road No.328, Kasatori Mountain, Kumakogen, Ehime Prefecture, Japan. Species-level identification of Sweltsa species is unsure, due to taxonomic problems in Japanese representatives of the genus (Shimizu et al. 2005). All specimens were fixed in 100% ethanol and stored at –20°C. The examined specimens of H. japonica (No. VHem-0114) and Sweltsa sp. (No. VHem-0116) are deposited in the Entomological Museum of Henan Institute of Science and Technology (HIST), Henan Province, China. The full genomic DNA was isolated from thoracic muscle using the QIAamp DNA Blood Mini Kit (Qiagen, Germany) and maintained at –20°C prior to PCR analysis.

The full mitogenomes of both Chloroperlidae species were amplified and sequenced as described in previous studies (Li et al. 2017; Wang et al. 2017a, 2017b; Liu et al. 2018) and uploaded in GenBank with the accession numbers: OL351265 (Haploperla japonica) and OL351266 (Sweltsa sp.) (Table 1). The program BioEdit version 7.0.5.3 was used to assemble raw sequences into contigs (Hall 1999). tRNA genes were identified with the software ARWEN using default setting (Laslett and Canbäck 2008). PCG and two rRNA genes were identified using BLAST searches in NCBI (http://www.ncbi.nlm.nih.gov), and confirmed by alignments of homologous genes from previously published stonefly mitogenomes. MEGA 5.0 was used to calculate nucleotide composition and codon usage for PCGs (Tamura et al. 2011). Strand asymmetry was measured by the following formulas: AT skew= [A−T]/ [A+T] and GC skew= [G−C]/ [G+C] (Perna and Kocher 1995). Tandem repeats in the A+T-rich region were identified using Tandem Repeats Finder (http://tandem.bu.edu/trf/trf.advanced.submit.html).

We used forty previously published and the two newly sequenced mitogenomes for phylogenetic analyses (Table 1). Two species from the family Nemouridae, Amphinemura yao Mo, Yang, Wang and Li, 2017 and A. longispina (Okamoto, 1922), were used as outgroup taxa (Table 1). Each PCG was aligned separately using the MAFFT algorithm implemented the TranslatorX online platform (Abascal et al. 2010). Poorly aligned regions were masked from the protein alignment before back-translation to nucleotides using GBlocks (implemented in the TranslatorX). Each rRNA gene was aligned separately by MAFFT v7.0 using G-INS-I strategy (Katoh and Standley 2013), and ambiguously aligned positions masked by GBlocks v0.91b (Castresana 2000).

We assembled the “PCG13 matrix” (11,049 bp in total) for the phylogenetic analyses, including 13 PCGs. GTR+I+G was the best-fit model for the nucleotide sequence alignments (jModelTest 0.1.1, Posada 2008). Bayesian analyses were conducted on the PCG13 dataset partitioned by gene was conducted Bayesian inference (BI), using MrBayes 3.2.6 and maximum likelihood (ML) using RAxML-HPC2 8.1.11 (Stamatakis 2006; Ronquist et al. 2012). For the BI analyses, two simultaneous runs of 10 million generations were conducted, sampling every 1,000 generations and a burn-in rate of 25%. Stationarity was assessed using Tracer v.1.5 (Estimated sample size > 200). Bootstrap ML analyses of 1,000 replicates was performed with the fast ML method conducted in RAxML using the GTRGAMMA model for nucleotide sequences.

In this study, the full mitogenomes of two stonefly species in the subfamily Chloroperlinae: Haploperla japonica and Sweltsa sp., were studied for the first time (Table 2, Fig. 1). The new complete mitogenomes were 15,893 (Sweltsa sp.) and 16,012 bp (H. japonica) in length (Fig. 1, Table 2). Genome sizes of these mitogenomes were in the middle of the range compared to mitogenomes of other stoneflies, that range from 15,048 bp (Isoperla bilineata) (Chen et al. 2018) to 16,602 bp (Nemoura nankinensis) (Chen and Du 2017a). Control region variation is the primary factor influencing size variation in Plecoptera mitogenomes (Table 2). The typical gene content found in most insects were found in these species: 37 mitochondrial genes (13PCGs, 22tRNAs and 2 rRNAs) and a large control region (Fig. 1). Twenty-three genes (9 PCGs and 14 tRNAs) were transcribed on the majority strand (J-strand), and the remaining fourteen genes (4 PCGs, 8 tRNAs, and 2 rRNAs) were encoded on the minority strand (N-strand) (Fig. 1).

Overlap between gene regions ranged from one to eight bp in length. Both species had seven bp overlaps in (ATP6–ATP8) and (ND4–ND4L), “ATGATAA” and “TTAACAT”, respectively. This phenomenon was also found in the mitogenome of many insects (Yang et al. 2017; Cao et al. 2019a). An intergenic spacer sequence was found to be located between tRNA^Leu(UUR) and COII in the two species, ranging from 12 bp in Sweltsa sp. to 31 bp in Haploperla japonica (Table 3).

The A+T content in the Haploperla japonica and the Sweltsa sp. mitogenome is detailed in Table 4. The A+T content of Sweltsa sp. is similar to the base compositional biases reported for other chloroperlid species, that range from 66.7% (Suwallia talalajensis Zhiltzova, 1976) to 69.6% (Sweltsa longistyla (Wu, 1938)) (Chen and Du 2015b; Wang et al. 2018a), whereas A+T content in H. japonica is larger than other chloroperlids. The four major genome partitions (PCGs, tRNAs, rRNAs, and CR) and the whole mitogenome of both Chloroperlinae species were distinctly biased towards A and T. Among the four major partitions across Chloroperlinae species, PCGs consistently had the lowest A+T content, while the control region had the highest. A+T content of PCGs encoded on the minority strand (PCGs–N) was also higher than that of PCGs encoded on the majority strand (PCGs–J) (Table 4).

Both Chloroperlinae species showed positive AT-skew and negative GC-skew for the whole mitogenome. Total mitogenome of Haploperla japonica, the AT-skew and GC-skew were 0.065 and –0.223, respectively, whereas in Sweltsa sp., the AT-skew and GC-skew were 0.087 and –0.254, respectively. Slight T-skews and moderate C-skews were found in the majority strand PCGs while the reverse phenomenon was found for minority strand PCGs: both T-skews and G-skews were marked. For tRNAs and rRNAs, G-skews varied from slight to marked, while T-skews were slight to moderate. Modest A-skew and marked C-skew were found in the CR (Table 4).

For the majority strand PCGs, most known metazoan mitogenomes have positive AT-skew and negative GC-skew (Wei et al. 2010a). However, both Chloroperlinae species had negative AT-skews and negative GC-skews for the majority strand PCGs (Table 4). This difference has been noted for other published stonefly mitogenomes (Wu et al. 2014; Chen and Du 2017a, 2017c; Wang et al. 2017a, 2017b), as well as in other insects such as Philopteridae (bird lice), Aleyrodidae (whiteflies) and Braconidae (wasps); this phenomenon may be a result from replication direction charges (Wei et al. 2010a). For both Chloroperlinae species, the four major partitions (PCGs, tRNAs, rRNAs and CR) all had the same trends in the A+T content, AT and GC-skews which displayed similar patterns and were consistent with the common base composition biases of insect mitogenomes.

The total length of PCGs in the studied mitogenomes was 11,247 bp (Haploperla japonica) and 11,244 bp (Sweltsa sp.) (Table 2). In both mitogenomes from Chloroperlinae, typical start codons ATN (Met/Ile) were the most commonly used start codons. Additionally, TTG was inferred as a start codon for ND1, and GTG as the start codon for ND5 in both species. This phenomenon has been found in other stoneflies, and TTG is used in other insect species (Bae et al. 2004; James and Andrew 2006; Sheffield et al. 2008; Wei et al. 2010b; Elbrecht and Leese 2015; Elbrecht et al. 2015; Huang et al. 2015; Sproul et al. 2015; Chen et al. 2016; Chen and Du 2017c; Wang et al. 2017a, 2017b; Cao et al. 2019b). In both species, ATP6, ATP8, COII, COIII, CYTB, ND2, ND4 and ND4L started with ATG (Met), while COI, ND3 and ND6 started with ATT (Ile) (Table 5).

The most commonly used stop codon in both Chloroperlinae species was TAA, found in ATP6, ATP8, COI, COIII, ND2, ND4, ND4L and ND6. In both species, the stop codon TAG was used by CYTB, ND1 and ND3 (except CYTB of Sweltsa sp., which used TAA). COII and ND5 used a partial stop codon T in both species (Table 5). This phenomenon has also been recorded from other published stoneflies (James and Andrew 2006; Qian et al. 2014; Elbrecht and Leese 2015; Elbrecht et al. 2015; Huang et al. 2015; Wu et al. 2014; Chen et al. 2016; Chen and Du 2017c; Wang et al. 2017a, 2017b; Cao et al. 2019b). Truncated stop codons occur in many arthropod mitogenomes, and it is corrected by post-transcriptional polyadenylation (Ojala et al. 1981; Wolstenholme 1992). The two species had the same start and stop codons in ten genes (COI, COII, COIII, ATP8, ATP6, ND1, ND2, ND4, ND4L and ND5), while the remaining three genes (ND3, ND6 and CYTB) only shared either the same start or stop codon (Table 5).

The genome-wide AT bias is reflected in codon usage. RSCU (relative synonymous codon usage) was calculated to identify the predominant synonymous codon (Grantham et al. 1980; Meganathan et al. 2012) (Fig. 2). The high AT base compositional biases was shown in the RSCU statistics for both Chloroperlinae mitogenomes. Codons ending with an A or U were favoured both four-fold and two-fold degenerate codons. RSCU demonstrates that codon usage patterns of these two stonefly species was consistent with other published Plecoptera species. All most commonly used codons: AAU (Asn), AUU (Ile), UUA (Leu1), AAA (Lys), AUA (Met), UUU (Phe) and UAU (Tyr), were always biased towards A and T at third codon positions (Fig. 2, Table 5).

In the present study, the tRNA genes in both species ranged from 65 to 71 bp, and comprised 9.26% (1,483 bp in H. japonica) and 9.36% (1,488 bp in Sweltsa sp.) of the complete mitogenomes (Tables 2 and 3). tRNA^Ser(AGN), which could not be detected by ARWEN software, was identified by alignment with homologous genes from other published stoneflies (James and Andrew 2006; Chen et al. 2016; Wang et al. 2017a, 2018a, 2018b; Cao et al. 2019b; Guo et al. 2021). Most tRNAs had the typical cloverleaf structure, except tRNA^Ser(AGN) that lacks the dihydrouridine (DHU) arm that formed a simple loop. In most insect mitogenomes, it is common for tRNA^Ser(AGN) to lack the DHU stem (Cameron and Whiting 2007; Wan et al. 2013). The length of the different nucleotides of tRNAs between H. japonica and Sweltsa sp. ranged from 1 bp (tRNA^Ser(AGN), tRNA^Ser(UCN), tRNA^Val) to 13 bp (tRNA^His). By comparison, the tRNA^Cys and tRNA^His had the most variation with more than ten indels or substitutions; the nucleotide insertion-deletion appears in tRNA^Cys, tRNA^His, tRNA^Phe, tRNA^Thr, tRNA^Tyr (Fig. 3).

Based on the secondary structures of tRNA, eight mismatched base pairs are found in H. japonica, which are U-U (2 bp), A-A (3 bp) and A-C (3 bp) (Fig. 3). And they have 24 G-U pairs, 7 located in the AA, 7 in the DHU, 8 in the AC and 2 in the TΨC stem. In Sweltsa sp., there are nine mismatches in tRNAs: U-U (1 bp), A-A (5 bp), C-C (1 bp) and A-C (2 bp) (Fig. 3); the 29 G-U pairs are located in the AA (5 bp), the DHU (8 bp), the AC (12 bp) and the TΨC stem (4 bp).

Similar to the mitogenomes of most other insects, the large and small rRNA subunits (lrRNA and srRNA) in both species were located between tRNA^Leu(CUN) – tRNA^Val (lrRNA) and tRNA^Val – the control region (srRNA) (Fig. 1, Table 3). The lengths of lrRNA and srRNA of Haploperla japonica were 1324 bp and 776 bp, and the lengths of lrRNA and srRNA of Sweltsa sp. were 1325 bp and 799 bp (Table 2). Secondary structures of both lrRNA and srRNA were highly similar between these species (Pair-wise sequence identity was 87.82% (lrRNA) and 87.98% (srRNA)).

Haploperla japonica was used as model for comparison of rRNA secondary structures. lrRNA consists of 5 structural domains (I–II, IV–VI), while domain III was absent, as usual in arthropods (Fig. 4) (Cannone et al. 2002). Domains I, II and IV were variable between these species, while eight helices (H2064, H2347, H2395, H2455, H2507, H2520, H2547 and H2588) within domain V had high similarity (Wang et al. 2017a).

srRNA consists of three domains (Fig. 5). Domains I and II were more variable than domain III in both species, while Helix 1399 region was the most conserved (Wang et al. 2017a). Furthermore, Sweltsa sp. possessed 21 more bp nucleotides in Helix 577 than Haploperla japonica.

In a previous study, the A+T-rich region (CR) was found to contain essential elements for the initiation of transcription and for replication (Zhang and Hewitt 1997). Haploperla japonica mitogenome had an A+T-rich region of 1,118 bp, with an A+T content of 80.6%; while Sweltsa sp. it was 1,011 bp long, with an A+T content of 80.0% (Tables 2 and 4). The A+T-rich region of both species is located between the srRNA and tRNA^Ile gene cluster (Fig. 1). The following structural elements were found in the control region of H. japonica and Sweltsa sp. (Fig. 6): (1) a leading sequence (588 bp (H. japonica), 496 bp (Sweltsa sp.)) adjacent to srRNA; (2) two large tandem repeats (TRs), while (3) the remaining control region was 321 bp (H. japonica) and 372 bp (Sweltsa sp.).

Phylogenetic relationships were reconstructed for the clade Systellognatha using 42 species: 36 complete (including the two species herein studied) and 6 nearly complete mitogenome sequences (Table 1). Amphinemura yao Mo, Yang, Wang and Li, 2017 and A. longispina (Okamoto, 1922) were used as outgroup taxa. Bayesian inference (BI) and maximum likelihood (ML) analyses, of the aligned 13 PCGs dataset, generated trees with identical topologies (Fig. 7).

Systellognatha includes seven families: superfamily Perloidea with Chloroperlidae, Kathroperlidae, Perlidae, Perlodidae; and superfamily Pteronarcyoidea with Peltoperlidae, Pteronarcyidae, Styloperlidae; all except Kathroperlidae (no mitogenomes yet sequenced) were represented in our study. We obtained all families as monophyletic (colors in Fig. 7). Our analyses strongly supported Peltoperlidae as the sister group of a clade containing the five other families. Similar relationships were recovered in previous phylogenetic analyses (Chen et al. 2018; South et al. 2020, 2021; Cao et al. 2021). A sister-group relationship between Perloidea (monophyletic herein: Fig. 7) and Pteronarcyoidea was not supported, as the latter was obtained as paraphyletic, which is inconsistent with morphological hypotheses (Zwick 2000). Phylogenetic analyses based on morphological characters proposed a clade Peltoperlidae + Styloperlidae, with Pteronarcyidae as their sister group (Zwick 2000). According to the present studies, Styloperlidae is sister to Pteronarcyidae rather than Peltoperlidae albeit with poor nodal support (as in Ding et al. 2019Wang et al. 2019; Cao et al. 2019a; 2021). Perlidae was sister to Chloroperlidae + Perlodidae with high support (Fig. 7) (confirming Ricker 1952; Illies 1965; Chen et al. 2018; Wang et al. 2019, 2018c; South et al. 2020, 2021). Our analyses found relationships within Systellognatha to be (((Chloroperlidae + Perlodidae) + Perlidae) + (Styloperlidae + Pteronarcyidae)) + Peltoperlidae (Fig. 7). Recent phylogenetic analyses based on mitochondrial data have inferred the same topology (Chen et al. 2018; Cao et al. 2021), however, our results differ from the older morphology-based classification. The different results may be due to the limited mitogenomic data yet published about the Systellognatha, being especially scarce in Pteronarcyidae and Styloperlidae, and lacking in Kathroperlidae.

At the generic level, one pattern is the disjunct phylogenetic positions for Asian and North American members of the same genus (Fig. 7). For example, regarding Isoperla Banks, 1906, Pseudomegarcys japonica Kohno, 1946 was sister to the Asian I. eximia Zapekina-Dulkeit, 1975, rather than the North American I. bilineata (Say, 1823), the type species of genus Isoperla. Additionally, regarding Acroneuria Pictet, 1841, Sinacroneuria dabieshana Li and Murányi, 2014 was sister to the Asian species A. hainana Wu, 1938, rather than the North American A. carolinensis (Banks, 1905). Regarding their morphology, both Pseudomegarcys and Sinacroneuria are rather distinctive genera, in addition, Pseudomegarcys belongs to the subfamily Perlodinae Klapálek, 1909, while Isoperla belongs to Isoperlinae Frison, 1942. Inclusion of these genera in Isoperla and Acroneuria, respectively, is not likely, and additional molecular studies are needed to clarify their phylogenetic position.