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

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 = [A–T]/[A+T] and GC-skew = [G–C]/[G+C]. The mitochondrial genome sequence of F. biocellata was added to GenBank under accession number MK905206.

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, Lohman & Meier, 2010; 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 very similar A+T contents:67.2%, 66.4%, 70.4%, and 70.6%, respectively (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.