The complete mitochondrial genome of Flavoperla biocellata Chu, 1929 (Plecoptera: Perlidae) and the phylogenetic analyses of Plecoptera

Of the roughly 400 species of Perlidae in the world, most species are widely distributed in the northern hemisphere, but a few can be found in South Africa and South America. There are only five species in the genus Flavoperla of the family Perlidae in China. To gain a better understanding of the architecture and evolution of mitochondrial genome in Flavoperla, the entire mitochondrial genome (mitogenome) of a Chinese Flavoperla biocellata Chu, 1929 from family Perlidae (Insecta: Plecoptera) was sequenced. The 15,805-bp long mitochondrial genome of F. biocellata contained 37 genes, including 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs) and a putative control region (CR). The gene arrangement of F. biocellata was identical with that of other stoneflies and with the fly Drosophila yakuba. Most PCGs of F. biocellata used the standard ATN start codons and complete TAN termination codons. Twenty-one of the 22 tRNA genes exhibited cloverleaf secondary structures, but the dihydrouridine (DHU) arm of trnSer (AGN) was completely reduced. Phylogenetic analyses with both Bayesian inference (BI) and maximum likelihood methods (ML) generated similar topology, both supporting the monophyly of all stonefly families and the infraorder Systellognatha. The phylogenetic analysis based on mitochondrial genomic data from 30 stonefly species recovered a well-supported tree resolving higher-level relationships within Plecoptera. The northern hemisphere suborder Arctoperlaria divided into two groups, Euholognatha and Systellognatha. The southern hemisphere suborder Antarctoperlaria formed two clades: Eustheniidae+Diamphipnoidae and Austroperlidae+ Gripopterygidae; consistent with relationships proposed based on morphology. The final relationships within Plecoptera were recovered as (((Perlidae+(Perlodidae+Chloroperlidae))+(Pteronarcyidae+(Peltoperlidae+Styloperlidae))) +(Taeniopterygidae+(Capniidae+(Nemouridae+Notonemouridae))))+ (Gripopterygoidae+Eusthenioidae).


INTRODUCTION
The Plecoptera (stoneflies) is an ancient order of aquatic insects that is an important bioindicator of water quality (William & Felmate, 1992). It is also important for the phylogenetic reconstruction of insects. Plecoptera contains two suborders, Antarctoperlaria and Arctoperlaria, including 16 extant families (Zwick, 1973). However, the phylogenetic relationship among the families of the Plecoptera is still unresolved (Gibson et al., 2004;Cameron, 2014). The area of origin and sister-group of Plecoptera are inconclusive (Zwick, 2009;Yoshizawa, 2011;Simon et al., 2012;Letsch & Simon, 2013;Song et al., 2016). In 2009, Zwick tried to compare the hypotheses concerning the placement of stoneflies based on the morphological evidence, but failed (Zwick, 2009). Most molecular studies support Dermaptera as the sister group of Plecoptera, however, the transcriptome data suggests that Plecoptera may be one of the basal branches of Polyneoptera, and sister to a clade composed of all other orders except Dermaptera and Zoraptera (Ishiwata et al., 2011;Yoshizawa, 2011;Letsch & Simon, 2013;Wu et al., 2014;Song et al., 2016). Low taxon sampling within Plecoptera which reduced the reliability of their proposed phylogenies (Ding et al., 2019). To develop a better understanding of the phylogeny of Plecoptera, we sequenced the mitochondrial genome of F. biocellata, which was the first species sequenced in the genus Flavoperla from the family Perlidae of the Arctoperlaria. There are only five species of Flavoperla in China, and thus the sequence of the mitochondrial genome of F. biocellata will help us understand the architecture and evolution of the genus and will supplement data on the Perlidae. Recently, the number of complete mitochondrial genome sequences available has increased dramatically due to the development of next-generation sequencing technology. But the information on the Perlidae is still limited, with only five complete mitochondrial genomes sequenced, two from Kamimuria, and one each from Togoperla, Dinocras and Acroneuria. Mitochondrial genomes have been used in phylogenetic and evolutionary studies of insects, including the Plecoptera (Chen & Du, 2017a;Chen & Du, 2017b;Chen & Du, 2018;Chen et al., 2018;Wu et al., 2014;Wang, Cao & Li, 2017;Wang et al., 2019). To date, over thirty mitochondrial genomes of Plecoptera have been sequenced and reported, and their nucleotide composition and gene arrangement are highly conserved and identical with that of Drosophila yakuba (Wei & Chen, 2011;Chen & Du, 2015;Chen, Wu & Du, 2016;Chen & Du, 2017a;Chen & Du, 2017b;Chen & Du, 2018;Chen et al., 2018;Huang et al., 2015;Zhou et al., 2016;Wu et al., 2014;Elbrecht et al., 2015;Wang, Ding & Yang, 2016a;Sproul et al., 2015;Wang, Cao & Li, 2017;Wang, Wang & Yang, 2016b;Wang et al., 2019). The mitochondrial genome of Plecoptera has been verified as a double strand circular molecule that is nearly 16 kb in length and contains 13 PCGs, 22 tRNAs, two rRNA and a large non-coding control region, similar to other mitochondrial genomes (Clary & Wolstenholme, 1985;Steinberg & Cedergren, 1994;Gibson et al., 2004;Jia & Higgs, 2007;Cameron, 2014;Chen & Du, 2018;Chen & Du, 2018;Wang et al., 2019). We analyzed the nucleotide compositions of mitochondrial organizations, codon usages, tRNA structures of F. biocellata, and performed phylogenetic analyses of Plecoptera on the strength of the nucleotide sequences of available stonefly mitochondrial genomes.

Sample collecting and DNA extraction
Fresh specimens of F. biocellata were collected from Liyang, Jiangsu Province China (31.43 • N, 119.48 • E) in 2018. No endangered or protected species were involved in the study. All the samples were collected by sweeping the trees near streams, and identified using a microscope to identify the samples and were then stored in 100% ethanol and held at −20 • C. Genomic DNA was extracted with the Column mtDNAout kit (Tianda, Beijing, China) from the identified adults and stored at −20 • C for subsequent PCR amplifications.

PCR reaction and DNA sequencing
The mitochondrial genome of F. biocellata was amplified using a similar sequencing strategy to Chen & Du (2018). Firstly, the large and small ribosomal genes (rrnL and rrnS), cox1 and cox2 were amplified using four pairs of newly designed primers. Based on the four obtained sequences, two pairs of LA-PCR primers were designed to amplify two overlapping fragments covering the whole mitochondrial genome (rrnL-cox2 and cox1-rrnS). The remaining gaps were filled with specially designed primers. LA-PCR amplifications were performed with LA Taq DNA polymerase (Takara, Japan) using the following conditions: 2 min of initial denaturation at 93 • C, followed by 40 cycles at 92 • C for 10 s; 30 s of annealing at 54 • C; and 8 min of elongation at 68 • C (20 cycles), which increased 20 s/cycle in the final 20 cycles; and 10 min of final elongation at 68 • C. Table 1 summarizes all primers used in this study. Electrophoresis with 1.0% agarose gels was used to separate the products according to gene size of PCR reactions. Axygen DNA Gel Extraction Kit was used to purify the PCR products. All obtained purified PCR fragments were sequenced using sanger sequencing method by Map Biotech Company (Shanghai, China).

Mitochondrial genome annotation
The mitochondrial genome was assembled with CodonCode Aligner (http://www. codoncode.com/aligner/). BioEdit was used to align all used mitochondrial genomes. Alignment with other published Plecoptera mitochondrial genomes were performed to identify the PCGs and rRNAs. ORF finder was used to examine the boundaries of the PCGs. The mitochondrial genome map was depicted by CGView (http://stothard.afns.ualberta. ca/cgview_server/). ARWEN (http://mbio-serv2.mbioekol.lu.se/ARWEN/) was used to detect the tRNAs and show their secondary structures. MEGA v. 6.0 was used to calculate the nucleotide composition. Composition skew values were obtained with AT-skew =

Phylogenetic reconstructions
In this study, 30 stonefly mitochondrial genomes were used for the phylogenetic analyses; the outgroup was the mayfly Parafronurus youi (Insecta Ephemeroptera) whose GenBank accession number is EU349015 (Table 2). MAFFT was used to align the amino acid sequences of the 13PCGs and concatenate them excluding the stop codons (Vaidya, Katoh & Standley, 2013). PartitionFinder v.2.1.1 with BIC (Bayesian Information Criterion) was used to determine the best nucleotide substitution models (GTR+I+G) and the best partitioning schemes (partition by gene types) using a greedy search algorithm with unlinked branch lengths (Lanfear et al., 2016). BI analysis (Bayesian inferences) was conducted with MrBayes v.3.2.6, conditions are as follows: 20 million generations of runs sampled every 1,000 generations; four chains (three hot chains and one cold chain) with a burn-in of 25% trees (Ronquist & Huelsenbeck, 2003). Tracer v.1.5 was used to examine the stationarity of all runs (ESS > 200). Maximum likelihood (ML) analysis was performed using RAxML v.8 and had 1,000 bootstrap replicates (Stamatakis, 2014). FigTree v.1.4.2 was used to visualized all phylogenetic trees.

Mitochondrial genome structure and nucleotide composition
The F. biocellata mitochondrial genome was a double strand circular 15,805 bp molecule. The typical set of 37 genes was identified along with a 852 bp-long control region, of the total, 23 genes were found on the J-strand (majority strand) and 14 genes were on the N-strand (minority strand) (Fig. 1, Table 3). Altogether 54 overlapping nucleotides were found in 14 gene pairs, and their length ranged from 1 to 17 bp (Table 3). There were 84 intergenic nucleotides (IGNs) dispersed in 13 gene gaps; and the longest IGN was located between 16s (rrnL) and trnVal and was 28 bp long, similar to other stoneflies. In F. biocellata, the whole mitochondrial genomes, the PCGs, tRNAs and rRNAs have  (Table 3). A+T content of the 37 genes ranged from 53.0% in trnTyr to 83.6% in trnGlu (Table 3). The mitochondrial genome of F. biocellata has a positive AT-skew value and negative GC-skew value, biased for the A and C nucleotides.

PCG, tRNA and rRNA genes
When compared with other stoneflies, each PCG of F. biocellata had a similar A+T content and length (Table 3). Eleven PCGs began with the standard start codon ATN (or ATT or ATG), however, for nad5 GTG was the start codon, while for nad1 it was TTG. Eleven PCGs of the mitochondrial genome had complete termination codons (TAA or TAG), but cox1, cox2 terminated with an incomplete stop codon T. The RSCU (relative synonymous codon usage) value was calculated for F. biocellata, suggesting that the most frequently used codon was TTA (Leu), while CGC (Arg) was the least often used codon (Fig. 2, Table 4). The F. biocellata mitochondrial genome had 22 tRNA genes, these tRNA genes were a total of 1,490 bp long with an average A+T content of 70.4% (similar to that in other stoneflies). Each tRNA gene ranged from 64 to 72 bp. Most tRNAs of F. biocellata had a typical cloverleaf secondary structure (Fig. 3), however, the dihydrouridine (DHU) arm of trnSer (AGN) was reduced to a small loop. Anticodons of all tRNAs in F. biocellata   were identical with other species in Plecoptera. Thirty-three mismatched base pairs were detected in the tRNAs of F. biocellata, all of which were G-U pairs (Fig. 3). Two rRNAs in F. biocellata had a total length of 2,195 bp and an A+T content of 70.6%. The large ribosomal RNA gene (rrnL) was 1,351 bp with an A+T content of 70%, the small ribosomal gene (rrnS) was 844 bp long, with an A+T content of 71.1% (Table 3).

The control region
In the mitochondrial genome of insects, the most variable region is the control region, which is difficult to sequence because of its complicated secondary structures. The control region of F. biocellata was 852-bp long with an A+T content of 73.2%. All stoneflies have the conserved control region and the location was between rrnS and trnIle (Fig. 1).

Phylogenetic analyses
Based on the concatenated sequences of 13 PCGs from 30 stonefly mitochondrial genomes, we reconstructed the phylogenetic relationship of Plecoptera; using the mayfly species Parafronurus youi as the outgroup. BI and ML were used to generate two phylogenetic trees which have generally identical topological structures (Figs. 4 and 5). In both analyses, the monophyly of each family is highly supported. The monophyly of suborders Antarctoperlaria and Arctoperlaria were consistently recovered in both trees. The phylogeny within Antarctoperlaria shows two superfamilies: Eusthenioidea (including Eustheniidae and Diamphipnoidae) and Gripopterygoidea (including Austroperlidae and Gripopterygidae), which is consistent with the study of Ding et al. (2019). In the suborder Arctoperlidae, the sister-group relationship between Nemouridae and Notonemouridae were consistently supported. In Systelognatha, the superfamily Perloidea is monophyletic; and the Perlodidae and Chloroperlidae are sister groups, which is congruent with recent studies (Chen & Du, 2018;Wang et al., 2019). The superfamily Pteronarcyoidea is monophyletic, with Styloperlidae being the sister group to the Peltoperlidae.

DISCUSSION
The mitochondrial arrangement of the F. biocellata mitochondrial genome was identical with other stoneflies and the model insect Drosophila yakuba (Clary & Wolstenholme, 1985), which was regarded as the ancestral mitochondrial genome due to the conserved gene order in non-insect hexapods (Nardi et al., 2003) and crustaceans (Crease, 1999). The nucleotide skew was consistent in Plecoptera and most other insects (Wei et al., 2010).
The trnSer was had a small loop, which is a common phenomenon in stoneflies and other metazoans (Garey & Wolstenholme, 1989). Zwick proposed two suborders: Antarctoperlaria and Arctoperlaria, and Arctoperlaria further divided into Systellognatha and Euholognatha (Zwick, 1973). In our study, the mitogenomic phylogeny strongly supported this classification (Figs. 4 and 5). In our study, the phylogeny within Antarctoperlaria supports the morphological phylogeny proposed by Zwick (Zwick, 2000). Ding et al. (2019) used four gene inclusion/exclusion datasets: protein-coding genes alone, with third codon positions (PCG123), PCGs without third In both BI and ML, Nemouridae + Notonemouridae confirmed the northern origin of the currently southern hemisphere restricted Notonemouridae. Nemouridae currently is only found in the northern hemisphere, while Notonemouridae is only distributed in the southern hemisphere. Our phylogenetic analysis recovered the sister group relationship between Nemouridae and Notonemouridae, which supports the hypothesis that Notonemouridae was historically distributed in the northern hemisphere, but then migrated to the southern hemisphere and is now extinct from the northern hemisphere.
The difference in phylogenetic studies may have also resulted from the low numbers of samples; more sampling of the three families in future studies is expected to reconstruct a more robust phylogenetic relationship among the three families. Recent phylogenetic studies aimed to resolve the phylogeny of Plecoptera have demonstrated that the mitogenomic data were reliable and informative in inferring the inner phylogenetic relationships of Plecoptera. However, more comprehensive sampling especially for stoneflies from the southern hemisphere is needed to better resolve the mitochondrial phylogeny of Plecoptera.

CONCLUSIONS
The entire mitochondrial genome of F. biocellata was sequenced and analyzed. The gene arrangement of F. biocellata was very conserved and identical with other stoneflies and the mitochondrial genome of Drosophila yakuba. The phylogenetic reconstructions with Bayesian inference (BI) and maximum likelihood methods (ML) had similar phylogenetic topology. More sequences should be obtained in future works to better resolve the phylogeny of Plecoptera.