Molecular phylogeny inferred from the mitochondrial genomes of Plecoptera with Oyamia nigribasis (Plecoptera: Perlidae)

In this study, the mitochondrial genome of the stonefly, Oyamia nigribasis Banks, 1920 (Plecoptera: Perlidae), was sequenced and compared with the mtDNA genomes of 38 other stoneflies and two Ephemerae. The O. nigribasis mitogenome is a circular 15,923 bp molecule that encodes a large, noncoding control region (CR) and 37 typical mtDNA genes; these include 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), and two ribosomal RNA genes (rRNAs), respectively. Most of the PCGs initiated with ATN and terminated with TAN. The dihydrouridine (DHU) arm of tRNASer (AGN) was missing, whereas the other 21 tRNAs all exhibited the typical cloverleaf secondary structure. Stem-loop (SL) structures and tandem repeats were identified in the CR. Phylogenetic analyses using Bayesian inference and maximum likelihood were undertaken to determine relationships between stoneflies. Results indicated that the Antarctoperlaria, which contains Gripopterygidae, was absolutely separated from Arctoperlaria; this finding agrees with morphology. Finally, the overall relationships could be summarized as follows ((((Notonemouridae + Nemouridae) + Leuctridae) + (Scopuridae + (Capniidae + Taeniopterygidae))) + (((Perlodidae + Chloroperlidae) + Perlidae) + (Pteronarcyidae + (Peltoperlidae + Styloperlidae))) + ((Diamphipnoidae + Eustheniidae) + Gripopterygidae)).

The Plecoptera order of stoneflies is a basal infraorder of the Neoptera; it is dispersed worldwide (except Antarctica) and contains ancient hemimetabolous insects 1 . Stonefly larvae frequently inhabit clean rivers and streams and are quite sensitive to dirty, polluted environments; thus stoneflies are an important bioindicator of water quality [2][3][4] .
The typical metazoan mitochondrial genome includes a noncoding sequence called the control region (CR) and 37 genes including 13 protein coding genes (PCGs), 22 tRNAs, two rRNAs 5,6 . The analysis of mitochondrial genomes has profoundly influenced the genetics and taxonomy of insects [7][8][9][10][11][12] . The development of whole genome sequencing technologies for insects has been slow in contrast to mtDNA; however, mitochondrial barcodes and sequences are commonly used for insect species identification [13][14][15][16] . The mtDNA sequences of representative species in all Insecta orders have been deposited in GenBank and exceed 1,500 entries 17 .
Despite previous studies, conclusions on the phylogeny and biogeography of stoneflies are inconsistent, especially concerning the family composition in the highest systematic categories [18][19][20][21][22][23] . There is some controversy on the phylogeny of Plecoptera; thus, we obtained the mtDNA sequence of Oyamia nigribasis and constructed phylogenetic trees based on PCG sequences to deduce the phylogenetic relationships of 39 stonefly species.

Results and discussion
Genome annotation and base composition. The O. nigribasis mitogenome is a circular 15,923 bp molecule, and contains the typical set of 37 mtDNA genes (13 PCGs, 22 tRNAs, and two rRNAs) along with a noncoding control region (CR) of 1022 bp. Among the 37 genes, nine PCGs and 14 tRNAs were majority strand (J-strand); four PCGs, eight tRNAs, and two rRNAs were minority strand (N-strand) (Fig. 1, Table 1). The gene arrangement in O. nigribasis mtDNA was highly conserved with other sequenced stoneflies and identical to the mitogenome of Drosophila yakuba, which is regarded as the putative ancestral arthropod 24  www.nature.com/scientificreports/ encompassed 48 intergenic nucleotides (IGNs). The shortest overlap was 1 bp (multiple sites), whereas the longest was 9 bp and mapped between trnTyr (Y) and cox1 ( Table 1). The shortest interval between genes was 1 bp (two sites) while the longest was a 16-bp intergenic region between trnSer2 (UCN) and nad1 (Table 1). Similar A + T contents were observed for the entire O. nigribasis mtDNA molecule, PCGs, tRNAs, rRNAs, and CR, which were 70.2%, 69.1%, 70.2%, 71.5%, and 72.7%, respectively ( Table 2). The lowest and highest A + T content was 62.4% for cox3 and 88.1% for trnGlu (E), respectively ( Table 1). The AT-and GC-skew expressed positively and negatively, respectively, which is consistent with other stonefly mitogenomes (Table 2).
Protein-coding genes. O. nigribasis PCG were similar in size (65-71 bp), whereas the A + T content varied from 62.7-88.1% (Tables 1, 2). The majority of PCGs possessed the standard start codon ATN (ATT, ATC or ATG); however, nad2 started with GTG, which has also previously been observed for other stonefly such as Taeniopteryx ugola 10 . Furthermore, most PCGs terminated with complete codons (e.g. TAA, TAG); however, the stop colon in cox1, cox2 and nad5 terminated in a single T, which have also been previously reported for many stoneflies like Leuctra sp., Nemoura nankinensis, Taeniopteryx ugola and Doddsia occidentalis 9,10,12 . Such translation termination could be completed by post-transcriptional polyadenylation 9,10 . The relative synonymous codon usage (RSCU) values of TTA (Leu), TCT (Ser), and CCT (Pro) were relatively high, whereas TCG (Ser) and ACG (Thr) were used less frequently than other codons (Fig. 2, Table 3). Most species of stoneflies show a high RSCU value of leucine while the usage of other amino acids are diverse from each other 8-12 . Transfer RNA genes. The typical set of 22 tRNA genes was observed in the O. nigribasis mitogenome, and the combined length and mean A + T content was 1426 bp and 70.2%, respectively ( Table 2). Fourteen tRNAs were encoded in a clockwise orientation, whereas the remaining eight were transcribed counterclockwise (Fig. 1, Table 1). Apart from trnSer (AGN), where the dihydrouridine (DHU) arm was absent, and that is a very common feature of mitochondrial tRNA-Ser conserved in mammals and some insects [9][10][11][12]25 , the other 21 tRNAs exhibited the representative cloverleaf secondary structure (Fig. 3) that is typical of other metazoan mitogenomes. The tRNAs contained some mismatched base pairs, and many of these contained G-U pairs (Fig. 3).  Table 2). The two rRNA genes (rrnL, rrnS) generally map between trnLeu (CUN) and the CR, and this location was conserved in the mtDNA of O. nigribasis (Fig. 1). The full-length, intact rrnL was 1362 bp with an A + T content of 72.8%, whereas the 819 bp rrnS was truncated and had a 70.2% A + T content ( Table 1).
The non-coding control region. Mitogenome control regions are highly variable and exhibit variable lengths and nucleotide composition. The O. nigribasis CR was slightly larger than the CR in Plecopteran insects,    (Fig. 4). Four tandem repeats mapped between 15,525-15,661 bp. The remaining sequences in the CR were A + T rich (Fig. 4).

Phylogenetic analyses. The concatenated sequences of 13 PCGs from 38 additional Plecopteran species
were downloaded from GenBank. The mtDNAs from two Ephemeroptera species, Parafronurus youi and Isonychia ignota, served as outgroups (Table 5) because they are relatively close to stonefly in classification. ClustalX was used to align the amino acid sequences of the 13 PCGs, and MrBayes v. 3.1.2 and IQ-Tree v. 1.6.12 were utilized to generate the topology by Bayesian inference (BI) and maximum likelihood (ML) analysis, respectively.
The two trees showed a high degree of similarity (Fig. 5, 6), excluding individual species such as the ones in Scopuridae and Pteronarcyidae. However, it is important to note that the nodal support values of the BI tree (Fig. 5) were more credible based on a previous analysis by ML 11 . Sequence data of selected southern hemisphere families were analyzed; the suborder Antarctoperlaria, including Gripopterygidae, Diamphipnoidae and Eustheniidae, was separated from other stonefly families that affiliated with Arctoperlaria, which is distributed in the northern hemisphere. This finding differs from the conclusion that Gripopterygidae could not be separated from other Arctoperlarian families in Shen & Du, 2019 12 , while consistent with Ding, 2019 23 . The two clades of Arctoperlaria, Euholognatha and Systellognatha, were strongly supported at the family level as monophyletic clades. In the infraorder Euholognatha, it was explicit that Leuctridae was clustered with the group of Nemouridae + Notonemouridae and Taeniopterygidae was recovered as the sister group of Capniidae. However, the phylogenetic relationship of Scopuridae, which need more datasets and independent evidence, was difficult to determine. Scopuridae was close to Taeniopterygidae and Capniidae based on BI analysis but clustered with other Euholognatha in the ML tree. From the perspective of Systellognatha, the monophyletic relationships  www.nature.com/scientificreports/ in the superfamily Perloideae could be highly advocated as (Perlidae + (Perlodidae + Chloroperlidae)), even though marginal divergence has been reported 11,16 . However, the phylogenetic relationship within the superfamily Pteronarcyoidea is more controversial. As shown in the BI tree (Fig. 5), Styloperlidae was more closely related to Peltoperlidae and clustered with Pteronarcyidae, which was consistent with morphology but inconsistent with phylogeny of ((Pteronarcyidae + Styloperlidae) + Peltoperlidae) 11 . Pteronarcyidae was included in the clade containing Perloideae and clustered with Styloperlidae and Peltoperlidae in the ML tree. Similar discrepancies have been reported in related studies 23 and are potentially caused by the use of different algorithms and models. Increasing numbers of stonefly mtDNAs are undergoing sequence analysis. Thus, it is likely that controversial phylogenetic relationships will eventually be resolved and the phylogeny of Plecoptera can be more accurately presented based on increased numbers of mitogenomes. It is worth looking forward to that more genes just like nuclear genes can also help to improve the phylogeny of stoneflies.

Methods
Sample preparation and mitogenome amplification. This study was conducted without harming protected or endangered species, and all research activities were authorized. Specimens of O. nigribasis were collected from Benxi (Liaoning Province, China; July, 2018) and preserved in 100% ethanol. DNA extraction was performed using instructions supplied with the Column mtDNAout kit (Tianda Beijing, China). Universal or specifically-designed primers were used to amplify mitochondrial genes in long overlapping fragments (Table 4). LA-PCR and consecutive specific PCR amplifications were conducted using conditions described previously 10