Two Complete Mitochondrial Genomes From Leuctridae (Plecoptera: Nemouroidea): Implications for the Phylogenetic Relationships Among Stoneflies.

The family-level relationships within Plecoptera have been a focused area of research for a long time. Its higher classification remains unstable, and the phylogenetic relationships within Plecoptera should be re-examined. Here, we sequenced and analyzed two complete mitochondrial genomes (mitogenomes) of Paraleuctra cercia and Perlomyia isobeae of the family Leuctridae. We reconstructed the phylogeny of Plecoptera based on 13 protein-coding genes (PCGs) from published stoneflies. Our results showed that the Bayesian inference and maximum-likelihood tree had similar topological structures except for the positions of two families, Peltoperlidae and Scopuridae. The Plecoptera is divided into two clades, the suborder Antarctoperlaria and the suborder Arctoperlaria. The two suborders subsequently formed two groups, Eusthenioidea and Gripopterygoidea, and Euholognatha and Systellognatha, which is consistent with the results of morphological studies. In addition, the Leuctridae is the earliest branch within the superfamily Nemouroidea. But the monophyly of Perloidea and Pteronarcyoidea are still not well supported.

Stoneflies (order: Plecoptera) are an ancient group that includes approximately 4,000 extant species worldwide (DeWalt et al. 2020). Their larval habitats include cool lakes and pools, clean rocky streams and rivers, or sandy places. Larval instars are sensitive to water pollution; therefore, stoneflies serve as biological indicators of clean water (Stewart and Stark 2002, Fochetti and Figueroa 2008, Stewart and Stark 2008). The monophyly of the order Plecoptera is supported by some derived features, such as gonads forming loops, anterior ends of left and right ovaries and testes (Zwick 1973). The order is considered to be divided into two suborders, Antarctoperlaria and Arctoperlaria (Zwick 2000). The family-level relationships within Plecoptera have been studied for more than a century. The unstable relationship in higher classification has been raised by long-term controversial taxonomic hypotheses based on primitive (archaic) and derived (specialized) characters (Klapálek 1905; Tillyard 1921; Ricker 1950, 1952; Illies 1965), and fossil wing venation and proportions of body parts (Sinitshenkova 1987). Zwick (2000) summarized the phylogeny of Plecoptera and reconstructed the generally accepted phylogenetic relationship within Plecoptera based on the morphological data. However, more conflicting hypotheses have been proposed when using the single or multiple gene and mitochondrial genome (mitogenome) data in these phylogenetic analyses. For instance, recent studies have produced conflicting results on the relative position of the family Leuctridae (Thomas et al. 2000, Terry and Whiting 2003, Ding et al. 2019), whereas others have suggested that superfamily Pteronarycoidea is not monophyletic (Terry and Whiting 2003; Davis 2013; Chen and Du 2017a,b, 2018; Chen et al. 2018; Veale et al. 2019; Wang et al. 2019).

Currently, mitogenomes are widely used for phylogenetic analysis of insect lineages. Compared with single- or multi-gene analyses, mitogenomes can provide ample nucleotide sites and genome-level characters, such as the relative position of different genes, RNA secondary structures, and modes of control of replication and transcription (Boore 2006, Masta and Boore 2008). Collection of more nucleotide sites may help to increase resolution in phylogenetic analyses (Okajima and Kumazawa 2010). Nowadays, many new mitochondrial genomes of stonefly species (covering all the families) have been added to Plecoptera.

In this article, we sequenced and described the microgenome of two Leuctridae (Paraleuctra cercia and Perlomyia isobeae). We established a family-level phylogeny within Plecoptera, including 50 mitogenomes of stonefly species representing 16 families based on two data sets (PCGs: 13 protein-coding genes; PCG12R: the first and second codons of PCGs and two rRNAs).

Materials and Methods

Specimens and DNA Extraction

Adult male specimens of P. cercia were collected from Yuzawa town (Niigata Prefecture, Japan; April 2016) and of P. isobeae from Judzuru-wa-ikoi Park (Minamiawaji city, Japan; April 2016). Specimens (P. cercia with no. VHem-00091 and P. isobeae with no. VHem-00093) were stored in 100% ethanol and maintained at −20°C. All the voucher specimens are preserved in Henan Institute of Science and Technology (HIST), Henan Province, China. The thoracic muscles of specimens were used to extract the genomic DNA. Finally, we extracted the genomic DNA using the Qiagen DNAeasy Blood & Tissue Kit (Qiagen, Germany) according to the manufacturer’s instructions.

Sequencing, Mitogenome Assembly, and Annotation

Illumina Hiseq 2500 with 500 cycles of paired-end sequencing (250 bp reads) was performed at BerryGenomics Co., Ltd., Beijing, China. Two newly sequenced mitogenomes have been deposited in GenBank (Table 1). Raw data files were proofread and then assembled into contigs using BioEdit v7.0.5.2 (Hall 1999). Transfer RNA genes (tRNAs) were identified by MITOS WebServer (Bernt et al. 2013), with the parameter of genetic code set to invertebrate. All the PCGs and ribosomal RNA genes (rRNAs) were verified manually by the comparing sequences with other published stonefly mitogenomes. The control region was identified by the boundary of the rRNA genes. MEGA version 6.06 software was used to calculate the nucleotide composition (Tamura et al. 2013). The AT and GC skews were calculated using the following formulae: AT skew = (A − T)/(A + T); GC skew = (G − C)/(G + C) (Perna and Kocher 1995).

Table 1.

List of taxa used in this research

Order	Superfamily	Family	Species	Accession number
Plecoptera	Perloidea	Perlidae	Acroneuria hainana	NC_026104
			Caroperla siveci	MG677942
			Sinacroneuria dabieshana	MK492253
			Togoperla sp.	KM409708
			Kamimuria wangi	NC_024033
			Kamimuria chungnanshana	NC_028076
			Dinocras cephalotes	NC_022843
		Chloroperlidae	Suwallia bimaculata	KX091865
			Suwallia teleckojensis	MF198253
		Perlodidae	Perlodes sp.	MF197377
			Isoperla bilineata	MF716959
			Isoperla eximia	MG910457
			Pseudomegarcys japonica	MG910458
	Pteronarcyoidea	Pteronarcyidae	Pteronarcys princeps	NC_006133
			Pteronarcella badia	NC_029248
		Peltoperlidae	Cryptoperla stilifera	KC952026a
			Peltoperlopsis cebuano	MK387068
			Soliperla sp.	MF716958
		Styloperlidae	Styloperla sp.	KR088971a
			Styloperla spinicercia	KX845569
			Cerconychia sapa	MF100783a
	Nemouroidea	Capniidae	Apteroperla tikumana	NC027698
			Capnia zijinshana	KX094942
			Mesocapnia arizonensis	KP642637a
			Mesocapnia daxingana	KY568983a
		Leuctridae	Rhopalopsole bulbifera	MK111419a
			Paraleuctra cercia	MK492251
			Perlomyia isobeae	MK492252
		Nemouridae	Nemoura nankinensis	KY940360
			Nemoura papilla	MK290826
			Amphinemura longispina	MH085446
			Amphinemura yao	MH085447
			Indonemoura jacobsoni	MH085448
			Indonemoura nohirae	MH085449
			Mesonemoura metafiligera	MH085450
			Mesonemoura tritaenia	MH085451
			Protonemura kohnoae	MH085452
			Protonemura orbiculata	MH085453
			Sphaeronemoura grandicauda	MH085454a
			Sphaeronemoura hamistyla	MH085455a
		Notonemouridae	Neonemura barrosi	MK111418a
		Taeniopterygidae	Taeniopteryx ugola	MG589786
			Doddsia occidentalis	MG589787
		Scopuridae	Scopura longa	MH510071
	Gripopterygoidea	Gripopterygidae	Zelandoperla fenestrata	NC_034997
			Antarctoperla michaelseni	MK111413
		Austroperlidae	Klapopteryx armillata	MK111414a
	Eusthenioidea	Eustheniidae	Neuroperla schedingi	MK111415
		Diamphipnoidae	Diamphipnoa annulata	MK111416a
			Diamphipnopsis sp.	MK111417
Ephemeroptera (outgroup)	Ephemeroidea	Ephemeridae	Ephemera orientalis	NC_012645
	Heptagenioidea	Heptageniidae	Parafronurus youi	NC_011359
		Isonychiidae	Isonychia ignota	HM143892

Incomplete genome sequence. Two new sequenced species are marked in bold.

Phylogenetic Analysis

In this study, a total of 53 insect mitogenomes, including 50 plecopteran mitogenomes as ingroup and 3 ephemeropteran mitogenomes (Ephemera orientalis, Isonychia ignota, and Parafronurus youi) as outgroup, were selected for phylogenetic analyses (Table 1). We implemented the MAFFT within the TranslatorX online platform for each PCG (excluding stop codons) alignment, and the final nucleotide alignment was generated by using the Gblocks within the TranslatorX (Abascal et al. 2010). rRNA genes was aligned using the MAFFT v7.0 online server with the G-INS-I strategy (Katoh and Standley 2013). Phylogenetic analyses were based on the two data sets: 1) the ‘PCGs matrix’, including all connected PCGs, and 2) the ‘PCG12R matrix’, including the first and second codons of PCGs and two rRNAs. We conducted the maximum-likelihood (ML) analysis from IQ-TREE web server (Trifinopoulos et al. 2016), and the Bayesian inference (BI) analysis in MrBayes from the CIPRES Science Gateway v3.3 (Miller et al. 2010). The best-fit nucleotide substitution model (GTR + I + G) of each gene was determined by using ModelFinder from IQ-TREE with Bayesian information criterion (Trifinopoulos et al. 2016). We constructed the ML tree using an ultrafast bootstrap approximation approach with 10,000 replicates. The conditions of BI analysis were set as follows: two simultaneous runs of 10 million generations with sampling every 1,000 generations, and a burn-in rate of 25% trees.

Results and Discussion

General Features of the Two Mitogenomes

The complete mitogenome of P. cercia and P. isobeae were 15,625 and 15,795 bp long, respectively (Supp Fig. S1 [online only], Supp Table 1 [online only]). Both mitogenomes had 37 genes (including 13 PCGs, 22 tRNAs, and 2 rRNAs) and a control region, and retained the so-called ancestral mitogenome arrangement of insects (Boore 1999).

The A+T content was also much higher than the G+C content in all the gene categories of the mitogenomes of P. cercia and P. isobeae, respectively. Consistent with all 24 species of Nemouroidea sequenced to date, the A+T content, AT and GC skews, exhibited similar features (Supp Table 2 [online only]).

In the P. cercia mitogenome, all PCGs started with ATN (ATT and ATG) codons, except for ND6 and ND1 that initiated with ATG and TTG, respectively (Supp Table 1 [online only]). Similar to the P. cercia mitogenome, the most commonly used start codon in the P. isobeae mitogenome was ATN, whereas only ND1 and ND5 used TTG and GTG, respectively (Supp Table 1 [online only]). All PCGs used TAA and TAG as stop codons except COII and ND5. In P. cercia, both COII and ND5 stopped with an incomplete codon T. However, COII used a single T, and ND5 used TAA as stop codon in P. isobeae (Supp Table 1 [online only]).

Phylogenetic Relationships Within Plecoptera

The BI and ML tree of PCGs had similar topological structures except for the relative location of two families, Peltoperlidae and Scopuridae (Figs. 1 and 2), whereas BI and Ml analyses of PCG12R generated the same tree topologies (Supp Fig. 2 [online only]). Both BI and ML analyses using two data sets recovered two monophyletic suborders, Antarctoperlaria and Arctoperlaria, with suborder Arctoperlaria subsequently divided into Euholognatha and Systellognatha. Meanwhile, the suborder Antarctoperlaria was divided into two superfamilies, Eusthenioidea and Gripopterygoidea. The phylogenetic relationship within Antarctoperlaria was as follows: (Austroperlidae + [Gripopterygidae] + [Eustheniidae + Diamphipnoidae]). These results support the most widely accepted organization proposed by Zwick (2000).

Fig. 1.

Phylogenetic relationships among stoneflies inferred by Bayesian inference. Values at nodes are Bayesian posterior probabilities using the 13PCG data set. The tree was rooted with three outgroups (Ephemera orientalis, Parafronurus youi, and Isonychia ignota).

Fig. 2.

Phylogenetic relationships among stoneflies inferred by maximum-likelihood analysis. Values at nodes are bootstrap values using the 13PCG data set. The tree was rooted with three outgroups (Ephemera orientalis, Parafronurus youi, and Isonychia ignota).

Currently, the classification system of Euholognatha includes one superfamily (Leuctridae, Taeniopterygidae, Capniidae, Nemouridae, and Notonemouridae) and one family (Scopuridae) based on morphological characters. Both Zwick (2000) and Nelson (1984) listed Nemouridae and Notonemouridae as sister families, and Nemouridae + Notonemouridae is then placed as a sister group to Capniidae + Leuctridae. Meanwhile, the monotypic Scopuridae and the large superfamily Nemouroidea are supported to be the sister group (Zwick 2000). In this study, BI analyses of PCG data set show that the Leuctridae is the earliest branch within the infraorder Euholognatha (post-probability values [PP] = 1.00), whereas Scopuridae forms a sister group with two families (Capniidae and Taeniopterygidae) of the superfamily Nemouroidea (PP = 0.7; Fig. 1). For BI and ML analyses of PCG12R and ML analyses of PCGS, Leuctridae is the sister group to Scopuridae (bootstrap values [BS] = 48/67, PP = 0.93), and this lineage is assigned to be the sister of the remaining four families of Nemouroidea (Fig. 2). Our results are incongruent with traditional views. For example, Zwick (2000) placed two families of Scopuridae and Taeniopterygidae near each other as the two most basal lineages of Euholognatha. However, the position of Leuctridae in our analysis is similar with the molecular analyses of Terry and Whiting (2003) and Ding et al. (2019). However, we only included a single scopurid sequence in our study. Additional Scopuridae sequences should be required to better resolve the phylogenetic position of this family.

Our analysis supports the sister group relationship of Taeniopterygidae and Capniidae (BS = 94/80, PP = 1.00/0.83; Figs. 1 and 2 and Supp Fig. S2 [online only]). This result is inconsistent with previous morphological hypotheses (both Illies 1965 and Zwick 2000 placed Leuctridae as a sister taxon to Taeniopterygidae). Meanwhile, our analysis was not consistent with the molecular analysis of Thomas (2000), which supported Nemouridae as a sister group to the remainder of Plecoptera. Nevertheless, it is congruent with the results of our previous study (Cao et al. 2019, Wang et al. 2019) and molecular analyses of Chen and Du (2018) and Ding et al. (2019). In addition, the clade of Capniidae + Taeniopterygidae has a closer relationship with the group of Nemouridae + Notonemouridae, which is similar with another study generated by mitogenomic data (Ding et al. 2019). However, less support was present in BI and ML analyses (BS = 42/47, PP = 0.89/0.56), demonstrating that the phylogenetic relationship of these two groups may be unstable.

The infraorder Systellognatha comprises two superfamilies, Perloidea (Perlidae, Perlodidae, and Chloroperlidae) and Pteronarcyoidea (Pteronarcyidae, Peltoperlidae, and Styloperlidae). Based on morphological features, Zwick (2000) proposed the family-level phylogenetic relationship within the infraorder Systellognatha as (Pteronarcyidae + [Peltoperlidae + Styloperlidae]) + (Chloroperlidae + Perlidae + Perlodidae), and the phylogeny of the superfamily Perloidea was unresolved. Many molecular studies have attempted to reconstruct the phylogenetic relationships within Systellognatha, and more conflicting opinions were raised. For example, Thomas et al. (2000) recovered the relationship of Systellognatha as Pteronarcyidae + (Peltoperlidae + [Chloroperlidae + (Perlodidae + Perlidae)]) based on a single gene (18S) and wing structure data sets. Terry and Whiting (2003) generated similar results based on multiple makers, and they proposed a monophyletic Perloidea and a paraphyletic Pteronarcyoidea, recovered as Styloperlidae + (Peltoperlidae + [Pteronarcyidae + (Chloroperlidae + [Perlidae + Perlodidae])]). Davis (2013), using data on transcriptomes, found that Perloidea recovered as (Perlidae + [Chloroperlidae + Perlodidae]). In addition, mitogenome data have also been used in some studies and failed to reconstruct the monophyly of Perloidea and Pteronarcyoidea (Elbrecht et al. 2015; Chen et al. 2016, 2018; Chen and Du 2017b, 2018; Wang et al. 2019). In our study, the ML analysis of PCGs highly supported a monophyletic Pteronarcyoidea and a monophyletic Perloidea (BS = 91); however, the BI tree of PCGs recovered a paraphyletic Pteronarcyoidea. Both BI and ML trees of PCGs showed the family-level relationship of Perloidea is Perlidae + (Perlodidae + Chloroperlidae). However, both BI and ML analyses of PCG12R still do not support the monophyly of Perloidea and Pteronarcyoidea. Our results are consistent with some previous studies (Chen and Du 2017b, 2018; Chen et al. 2018; Wang et al. 2018a,b, 2019), suggesting that more thorough sampling of taxa and other analytical methods should be added so as to obtain more accurate phylogenic relationships.

Click here for additional data file.