Comparative Mitochondrial Analysis of Cnaphalocrocis exigua (Lepidoptera: Crambidae) and Its Close Relative C. medinalis.

Rice leaffolders are important pests on rice in Asia, Oceania, and Africa, causing serious loss to rice production. There are two main rice leaffolders in China, namely Cnaphalocrocis medinalis (Guenée) and C. exigua (Butler) with the former having the ability of long-distance migration. To reveal the differences in the mitochondrial genomes (mitogenome) between them, we compared the completed mitogenome of C. exigua with three C. medinalis individuals. Although phylogenetic analysis based on the mitogenomic data strongly supported the close relationship between these two species, many differences were still being revealed. The results showed that the mitogenome of C. exigua was shorter in length (15,262 bp) and slight lower in AT content than that of C. medinalis. Except for the different start codons of nad3 and nad6 gene, we also found the cox1 gene had a typical start codon 'ATG' which suggested that the starting position of this gene must be reconsidered in the entire superfamily Pyraloidea. All tRNAs have a typical clover-leaf structure, except for the dihydrouridine (DHU) stem losing of trnS1, which has the atypical anticondon 'TCT' instead of 'GCT' in C. medinalis and most Pyraloidea species. Two intergenic regions (between trnY and cox1, nad3 and trnA) featured by AT repeats were only found in C. medinalis and even rarely appeared in reported Pyraloidea species. Furthermore, regardless of interspecific comparison or intraspecific comparison of these two species, protein coding genes, especially the atp8 genes, had quite different evolutionary rates.

Two main rice leaffolders, Cnaphalocrocis exigua and C. medinalis (Lepidoptera: Pyraloidea: Crambidae) are widely distributed in the tropical and temperate areas of Asia, Oceania, and Africa, which cause heavy loss on the rice production. They create longitudinal white and transparent streaks on the leaves by attaching the leaf margins together with silk strands and feeding inside (Khan et al. 1988; Park et al. 2010, 2014; Padmavathi et al. 2013; Yang et al. 2015).

Because the two species have similar morphology and habits, they are generally confused in the field in China (Pan 1984, 1985). The rice leaffolder C. medinalis had drawn much attention from researches worldwide for its ability to migrate and re-migrate over long distances (Gao et al. 2008, Huang et al. 2010, Wang et al. 2010). To date, only few studies were proceeded on C. exigua, and there was no evidence that support C. exigua can migrate yet (Feng et al. 2017, Liao et al. 2017). Though C. exigua is not as widely distributed as C. medinalis in China, but it has dominated paddy fields in many places of Sichuan Basin in southwest China which suggests it may become a principal rice pest (Pan 1984, Gao et al. 2008, Yang et al. 2015).

The mitogenome has been widely applied as an useful molecular marker for diverse evolutionary studies among species including phylogenetic, population genetics, and comparative and evolutionary genomics (Harrison 1989, Boore 1999, Babbucci et al. 2014, Cameron 2014). Mitochondria are the energy-producing organelles in eukaryotic cells (Scott et al. 2011). The excellent flight ability of C. medinalis suggests it is a qualifier with developed energy metabolic level. Thus, we infer the mitogenome of C. medinalis differs from that of C. exigua to some extent.

The Pyraloidea, with more than 15,570 species described worldwide, is one of the largest superfamilies in Lepidoptera and shows the most diverse life history adaptations (Nieukerken et al. 2011, Regier et al. 2012). Up to now, only about 40 complete or nearly complete mitogenomes from Pyraloidea species have been sequenced. A better understanding of Pyraloidea mitogenomes also requires an expansion of the taxon and genome samplings.

In this study, we report the complete mitogenome of C. exigua and compare with that of its close relative C. medinalis. Detailed mitogenomic information of these two important rice pests may help us further understanding the mechanisms of migration of C. medinalis and creating new ways to control these pests. Moreover, phylogenetic analysis was also conducted to further study the phylogenetic relationships within Pyraloidea species.

Materials and Methods

DNA Sample Extraction

The overwintering C. exigua pupae were collected from rice stubble fields in Qianwei county (Leshan, Sichuan Province, China; E103° 93′, N29° 21′) in March 2016 and reared in an incubator (LAC-250HPY-2, Shanghai Longyue Instruments, Shanghai, China) under constant conditions (26 ± 1°C, 80 ± 5% RH, and a photoperiod of 14:10 (L:D) h) for two generations (Liao et al. 2017). The newly emerged adult were preserved in ethanol (100%) and stored at −80°C till the isolation of DNA. Total genomic DNA was isolated from an individual of C. exigua using the MiniBEST Universal Genomic DNA Extraction Kit (Ver. 5.0 TaKaRa, Japan), according to the manufacturer’s protocol. The quality of isolated DNA was examined by 0.7% (w/v) agarose gel electrophoresis and then the DNA was used to amplify the entire mitogenome of C. exigua.

PCR Amplification, Cloning, and Sequencing

To obtained the whole genome sequence, two-pair of primers (Ce12Sr: GAA AGC GAC GGG CAA TAT GT and CeCO1r722: TAA ACT TCT GGA TGW CCA AAA AAT CA, and 12Sai: AAA CTA GGA TTA GAT ACC CTA TTA T (Simon et al. 1994) and CeCOIf14: AYT CWA CAA ATC ATA AAG ATA TTG G) were designed to amplify one long fragment cox1-rrnS and one short fragment rrnS-cox1, respectively. PCR reaction mixture (50 µl in total) included 10 µl 5× PrimeSTAR GXL Buffer, 5 µl dNTP mixture (2.5 mM each), 1 µl PrimeSTAR GXL DNA Polymerase (1.25 U/µl, TaKaRa, Japan), 1 µl primer (10 µM) each, 2 µl template DNA, and 30 µl double distillated water. The amplification program for this two fragments was performed in Easycycler Gradient 96 (Analytik Jena AG, Germany) under the following conditions: an initial denaturation at 98°C for 10 s, followed by 32 cycles, at 98°C for 10 s, 56°C for 15 s, 68°C for 8 min (5 min for the short fragment rrnS-cox1), and a final extension one cycle at 68°C for 5 min. The PCR products were subjected to 0.7% (w/v) agarose gel electrophoresis and purified using agarose Gel DNA Extraction Kit (MiniBEST Ver. 4.0, TaKaRa, Japan). Then, these two purified fragments were directly sequenced step by step through BGI Genome-sequencing firm. To minimize in vitro amplification errors, each PCR had four replications. The complete mitochondrial genome sequence of C. exigua obtained in this study was deposited in GenBank under the accession numbers MN877384.

Genome Annotation and Secondary Structure Prediction

Protein-coding genes (PCGs) and rRNA genes were identified based on homologous regions of previously sequenced Pyraloidea mitogenomes using the Clustal X version 2.0 (Larkin et al. 2007), and the GeneDoc version 2.6 software. The base composition and relative synonymous codon usage (RSCU) values were measured by MEGA version 7.0 program (Kumar et al. 2016). Composition skewness was calculated by the following formulas: AT skew =[A−T]/[A+T], GC skew=[G−C]/[G+C] (Perna and Kocher 1995). The tRNA genes were identified by the tRNAscan-SE software (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan 2016), or predicted by sequence features of being capable of folding into the typical clover-leaf secondary structure with legitimate anticodon. The Map of the mitogenome was drawn by using the GCView Server (http://stothard.afns.ualberta.ca/cgview_server/).

Comparative Analysis of Two Leaffolders Mitogenomes

Firstly, we compared the mitogenome of C. exigua (Ce) with three mitogenomes of C. medinalis (Cm1, collected from Yangzhou, Jiangsu Province, China, GenBank JN246082; Cm2, collected from Wulong county, Chongqing city, China, GenBank JQ305693; Cm3, collected from Cheongwon, Chunnam Province, Korea, GenBank JQ647917) (Chai et al. 2012, Wan et al. 2013, Yin et al. 2014) in respect of genome structure, base composition, PCGs, tRNA genes, rRNA genes, intergenic spacers, gene overlaps, and control regions. Subsequently, in order to analyze the intraspecific and interspecific genetic distances of mitochondrial genes of C. exigua and C. medinalis, the nucleotide and amino acid sequence were analyzed in MEGA version 7.0 program using the Kimura two-parameter and P-distance model, respectively. Lastly, the software packages DnaSP 5.0 (Librado and Rozas 2009) was used to calculate the number of synonymous substitutions per synonymous site (Ks) and the number of nonsynonymous substitutions per nonsynonymous site (Ka) for each PCG of the two leaffolders’ mitogenomes.

Phylogenetic Analysis

In addition to the sequenced mitogenome of C. exigua, 39 mitogenome sequences of 37 Pyraloidea species were downloaded from the NCBI database (Supp Table 1 [online only]). Two Tortricoidea species, Cydia pimonella (GenBank JX407107) and Adoxophyes orana (GenBank JX872403) were served as outgroups. The amino acid and nucleotide sequences of each of the 13 PCGs were aligned by Muscle through MEGA 7.0 software (Kumar et al. 2016), and then the concatenated set of nucleotide sequences from all 13 PCGs was used for phylogenetic analyses using Bayesian inference (BI) and maximum likelihood (ML) methods. The nucleotide substitution model test performed in MEGA version 7.0 showed that the GTR+I+G was the best-fitting model for the concatenated data. The BI analyses were implemented in MrBayes 3.1.2 (Huelsenbeck and Ronquist 2001) with four MCMC chains running for 2,000,000 generations and sampling every 2,000 generations. After discarding the first 25% samples as burn-in, Bayesian posterior probability values were calculated in a consensus tree. The ML analyses were performed using MEGA version 7.0 with 500 bootstrap replicates.

Results and Discussion

General Features of the C. exigua Mitogenome

The mitogenome of C. exigua was a closed-circular molecule of 15,262 in length (Table 1), which was within the range observed in the sequenced other Pyraloidea mitogenomes available with the size ranging from 15,110 bp in Maruca testulalis (Zou et al. 2016) to 15,594 bp in Ephestia kuehniella (Zhu et al. 2018). Annotation of the C. exigua mitogenome indicated that the structure and orientation of the 13 PCGs, 22 tRNA genes, 2 rRNA genes, and a noncoding region (A+T-rich control region) were typical of and similar to the other closely related pyraloid species (Yang et al. 2018, Zhu et al. 2018). Nine PCGs (nad2, nad3, nad6, cox1-3, atp6, atp8, cytb), 14 tRNAs, and the control region were located on the major J-strand. Fourteen genes (4 PCGs: nad5, nad4, nad4l, and nad1; both rRNA genes; and 8 tRNA genes) were located on the minor N-strand (Fig. 1). Distinct to the predicted ancestral gene order for insects, trnM-trnI-trnQ instead of trnI-trnQ-trnM was recognized in C. exigua, which was consistent to most Ditrysia species of Lepidoptera (Timmermans et al. 2014).

Table 1.

Summary of the mitogenomes of Cnaphalocrocis exigua (Me) and C. medinalis (Cm)

Feature	Strand	Start–end				Size				Intergenic nucleotides				Anti/start codon		Stop codon
		Ce	Cm1	Cm2	Cm3	Ce	Cm1	Cm2	Cm3	Ce	Cm1	Cm2	Cm3	Ce	Cm	Me/Cm
tRNA-Met (M)	J	1–68	1–68	1–68	1–68	68	68	68	68	0	0	0	0	CAT	CAT
tRNA-Ile (I)	J	69–135	69–133	69–133	69–133	67	65	65	65	0	0	0	0	GAT	GAT
tRNA-Gln (Q)	N	133–201	131–199	131–199	131–199	69	69	69	69	-3	-3	-3	-3	TTG	TTG
ND2	J	268–1,281	274–1,290	276–1,289	282–1,295	1,014	1,017	1,014	1,014	66	74	76	82	ATT	ATT	TAA
tRNA-Trp (W)	J	1,288–1,355	1,296–1,363	1,295–1,362	1,301–1,368	68	68	68	68	6	5	5	5	TCA	TCA
tRNA-Cys (C)	N	1,348–1,414	1,356–1,423	1,355–1,422	1,361–1,428	67	68	68	68	-8	-8	-8	-8	GCA	GCA
tRNA-Tyr (Y)	N	1,424–1,489	1,429–1,495	1,428–1,494	1,434–1,500	66	67	67	67	9	5	5	5	GTA	GTA
COX1	J	1,494–3,027	1,537–3,067	1,552–3,082	1,556–3,086	1,534	1,531	1,531	1,531	4	41	57	55	ATG	CGA	T-
tRNA-Leu (UUR)	J	3,028–3,094	3,068–3,134	3,083–3,149	3,087–3,153	67	67	67	67	0	0	0	0	TAA	TAA
COX2	J	3,095–3,776	3,135–3,816	3,150–3,831	3,154–3,835	682	682	682	682	0	0	0	0	ATG	ATG	T-
tRNA-Lys (K)	J	3,777–3,847	3,817–3,887	3,832–3,902	3,836–3,906	71	71	71	71	0	0	0	0	CTT	CTT
tRNA-Asp (D)	J	3,856–3,921	3,899–3,964	3,914–3,984	3,918–3,988	66	66	71	71	8	11	11	11	GTC	GTC
ATP8	J	3,922–4,080	3,965–4,132	3,985–4,149	3,989–4,153	159	168	165	165	0	0	0	0	ATT	ATT	TAA
ATP6	J	4,074–4,748	4,126–4,800	4,143–4,817	4,147–4,821	675	675	675	675	-7	-7	-7	-7	ATG	ATG	TAA
COX3	J	4,748–5,536	4,800–5,588	4,824–5,612	4,828–5,616	789	789	789	789	-1	-1	6	6	ATG	ATG	TAA
tRNA-Gly (G)	J	5,539–5,603	5,591–5,655	5,615–5,679	5,619–5,683	65	65	65	65	2	2	2	2	TCC	TCC
ND3	J	5,604–5,957	5,656–6,009	5,680–6,033	5,684–6,037	354	354	354	354	0	0	0	0	ATC	ATT	TAA
tRNA-Ala (A)	J	5,967–6,032	6,056–6,121	6,066–6,131	6,082–6,147	66	66	66	66	9	46	32	44	TGC	TGC
tRNA-Arg (R)	J	6,032–6,098	6,123–6,189	6,133–6,199	6,149–6,215	67	67	67	67	-1	1	1	1	TCG	TCG
tRNA-Asn (N)	J	6,099–6,164	6,190–6,256	6,200–6,265	6,216–6,281	66	67	66	66	0	0	0	0	GTT	GTT
tRNA-Ser (AGN)	J	6,165–6,230	6,259–6,324	6,266–6,331	6,282–6,347	66	66	66	66	0	2	0	0	TCT	GCT
tRNA-Glu (E)	J	6,235–6,300	6,326–6,393	6,333–6,399	6,349–6,415	66	68	67	67	4	1	1	1	TTC	TTC
tRNA-Phe (F)	N	6,299–6,366	6,392–6,460	6,398–6,466	6,414–6,482	68	69	69	69	-2	-2	-2	-2	GAA	GAA
ND5	N	6,367–8,101	6,461–8,195	6,467–8,201	6,483–8,217	1,735	1,735	1,735	1,735	0	0	0	0	ATT	ATT	T-
tRNA-His (H)	N	8,102–8,169	8,196–8,261	8,202–8,267	8,218–8,283	68	66	66	66	0	0	0	0	GTG	GTG
ND4	N	8,170–9,509	8,262–9,601	8,268–9,607	8,284–9,623	1,340	1,340	1,340	1,340	0	0	0	0	ATG	ATG	TA-
ND4L	N	9,510–9,803	9,601–9,894	9,607–9,900	9,623–9,916	294	294	294	294	0	-1	-1	-1	ATG	ATG	TAA
tRNA-Thr (T)	J	9,806–9,870	9,897–9,961	9,903–9,968	9,919–9,984	65	65	66	66	2	2	2	2	TGT	TGT
tRNA-Pro (P)	N	9,871–9,937	9,962–10,027	9,969–10,034	9,985–10,050	67	66	66	66	0	0	0	0	TGG	TGG
ND6	J	9,940–10,476	10,030–10,566	10,037–10,573	10,053–10,589	537	537	537	537	2	2	2	2	ATC	ATT	TAA
CYTB	J	10,481–11,629	10,573–11,721	10,580–11,728	10,596–11,744	1,149	1,149	1,149	1,149	4	6	6	6	ATG	ATG	TAA
tRNA-Ser (UCN)	J	11,628–11,693	11,723–11,790	11,727–11,792	11,743–11,808	66	68	66	66	-2	1	-2	-2	TGA	TGA
ND1	N	11,709–12,647	11,807–12,745	11,809–12,747	11,825–12,763	939	939	939	939	15	16	16	16	ATG	ATG	TAA
tRNA-Leu (CUN)	N	12,649–12,717	12,747–12,813	12,749–12,817	12,765–12,833	69	67	69	69	1	1	1	1	TAG	TAG
l-rRNA	N	12,718–14,074	12,814–14,202	12,818–14,192	12,834–14,207	1,357	1,384	1,375	1,374	0	0	0	0
tRNA-Val (V)	N	14,075–14,140	14,203–14,268	14,193–14,258	14,208–14,273	66	66	66	66	0	0	0	0	TAC	TAC
s-rRNA	N	14,141–14,922	14,269–15,049	14,259–15,037	14,274–15,025	782	781	779	752	0	0	0	0
Control Region	J	14,923–15,262	15,050–15,388	15,038–15,377	15,026–15,368	340	339	340	343	0	0	0	0

Cm1 indicates GenBank JN246082; Cm2 indicates GenBank JQ305693; Cm3 indicates GenBank JQ647917.

Fig. 1.

Circular mitochondrial genome map of C. exigua and C. medinalis. Genes encoded on the heavy strand or light strand are represented inside or outside of the circular mitochondrial genome map. The tRNA genes are abbreviated by triple letter, with Leu1 = CUN, Leu2 = UUR, Ser1 = AGN, and Ser2 = UCN. The second circle shows the GC content and the third shows GC skew calculated as (G-C)/(G+C). The GC content and GC skew are plotted as the deviation from the average value of the entire sequence. The circular mitochondrial genome map of C. medinalis was drawn based on the sequence GenBank JN246082.

The commonest start codon of PCGs was the ATG (in 8 PCGs: cox1, cox2, atp6, cox3, nad4, nad4l, cob, nad4), followed by three for ATT (nad2, atp8, nad5) and two for ATC (nad3, nad6) (Table 1). Nine PCGs had TAA stop codon, the remaining genes had incomplete stop codons (TA in nad4; T in cox1, cox2, and nad5).

Similar to other pyraloid species, the C. exigua mitogenome nucleotide composition was bias toward adenine and thymine (accounting for 81.6%: A = 40.5%, T = 41.1%, C = 10.9% G = 7.5%) (Yang et al. 2018, Zhu et al. 2018). The nucleotide bias was also reflected in the codon usage. Base composition at each codon position of concatenated 13 PCGs showed that the A+T content of the third codon positions were significantly higher than the first and second positions. In particular, T in each codon position of PCGs was over represented (Supp Table 2 [online only]).

All the 22 identified tRNA genes in C. exigua mitogenome were folded into the typical clover-leaf structure of mitochondrial tRNAs, except for trnS1 (AGN) in which the dihydrouridine (DHU) arm failed to form a stable stem-loop structure (Supp Table 3 [online only]). The exception trnS1 had been observed in many Pyraloidea and other insect mitogenomes (Zhang et al. 2016, Zhang and Ye 2017, Yang et al. 2018, Zhu et al. 2018). The length and base composition of these tRNAs were similar to other pyralids (Yang et al. 2018, Zhu et al. 2018). The number of base pairs in the DHU-stem ranged from 3 to 4. All the TΨC-stem of C. exigua had 5 base pairs except 3 bp in trnT, 6 bp in trnS1, and 4 bp in trnG, trnN and trnF (Supp Table 3 [online only]). Six unmatched base pairs including five U-U and one A-C pairs were identified in the stem region of tRNAs. The mismatch of base pairs existed widely in Pyraloidea mitogenomes (Chai et al. 2012, Yang et al. 2018).

Two ribosomal genes, lrRNA and srRNA, were 1357 bp and 782 bp, located between trnL1 (CUN) and trnV, and between trnV and the A+T-rich region, respectively (Table 1). Both the AT-skew and GC-skew of tRNA genes and rRNA genes were slightly positive in the two leaffolders (Supp Table 2 [online only]). The A+T-rich region of C. exigua mitogenome extended over 340 bp with extremely high AT content (94.4%) and was located between the rrnS and trnM genes (Table 1; Supp Table 2 [online only]).

Comparative Analysis of C. exigua and C. medinalis Mitogenomes

In order to reveal the difference between the two leaffolders, we compared the mitogenome of C. exigua with that of three individuals of C. medinalis. In general, the analysis results showed that differences between mitochondrial genomes exist not only between the two leaffolders but also within species of C. medinalis.

The PCGs

Except for cox1, nad3, and nad6 gene, the start codons of other genes were the same in the mitogenomes of the two leaffolders (Table 1). Including C. medinalis, most of the pyralis had ATT start condon for both nad3 and nad6 genes (Yang et al. 2018, Zhu et al. 2018). Besides C. exigua, the ATC start condon of nad3 could be found in Dichocrocis punctiferalis, Pseudargyria interruptella, and Scirpophaga incertulas, and the ATC start condon of nad6 was only found in Spoladea recurvalis (Wu et al. 2013, He et al. 2015, Song et al. 2016). Interestingly, for the cox1 gene, we found that it has the common start codon ‘ATG,’ which is different from the previously reported that ‘CGA’ as the start codon of this gene of most Pyraloidea species. We reanalyzed the starting position of the cox1 genes of all reported Pyraloidea species and found that other four species in family Crambidae (GenBank JX144892, KJ739310, KM453724, KC493629) also have the common start codon ‘ATG.’ In addition, although the three bases ‘TAG’ before the start codon ‘CGA’ appeared in many species of Pyraloidea, the ‘TTG’ in the same position could also be found in six species (including one of three individuals of C. medinalis) of family Crambidae (GenBank KP347977, KM244688, JF339041, KJ174087, KF859965, JQ305693) and two species of Galleriinae in family Pyralidae (GenBank HQ897685, KT750964). In view of the above reasons, we suggest that when annotating mitogenome of species in Pyraloidea, it is necessary to reconsider the starting position of this gene.

The distribution of codon families in two leaffolders were very similar, and both exhibited that leucine (Leu), Isoleucine (Ile), and phenylalanine (Phe) were the three most abundant codon families, whereas codon family Cysteine (Cys), arginine (Arg), and glutamate (Glu) were the three least abundant (Fig. 2A). This situation was in accordance with other Pyraloidea species (Dai et al. 2018, Zhu et al. 2018). Interestingly, the Phe amount of Cm1 was higher than that of Ce, Cm2, and Cm3. The RSCU analysis showed that total 5 and 7 codons could not be identified in the mitogenomes of C. exigua and C. medinalis, respectively (Fig. 2B). As reported before, missing codons ranged from 1 to 8 were found in most known Pyraloidea mitogenomes (Yang et al. 2018, Zhu et al. 2018). Particularly, the codons GCG and AGG were not presented in C. medinalis mitogenome but could be found in C. exigua mitogenome (Fig. 2B). Therefore, C. medinalis had more missing codons than C. exigua, and these missing codons showed high G/C content as the previously reported (Zhu et al. 2018), which may be one of the reasons that the A+T content of PCGs of C. medinalis is slightly higher than that of C. exigua.

Fig. 2.

The codon distribution (A) and the relative synonymous codon usage (RSCU) (B) of C. exigua (Ce) and C. medinalis (Cm) mitochondrial genomes. CDspT = condons per thousand codons. Codon families are plotted on the x-axis. Codons represented above the bar are not found in the mitogenomes. Cm1 indicates GenBank JN246082; Cm2 indicates GenBank JQ305693; Cm3 indicates GenBank JQ647917.

The genetic distance analysis of mitochondrial genes of C. exigua and C. medinalis revealed that regardless of 13 PCGs, 2 rRNA genes, and 22 tRNA genes, the difference between species was greater than that of within species. In both interspecific and intraspecific comparisons, the rrnL gene showed greater variability than rrnS gene and tRNA genes. Among 13 PCGs, the genetic distance of each gene between species was greater than 0.03, ranging from 0.0317 (nad4l) to 0.1973 (atp8), while the genetic distance within C. medinalis ranged from 0 (cox3, nad3) to 0.1139 (atp8) (Table 2). Notably, mitochondrial genes (amino acid sequences) of Cm2 and Cm3 were more homologous than that of Cm1 and Cm2 (8 identical PCGs versus 3 identical PCGs). Since Cm1 and Cm2 were both collected from China, while Cm3 was collected from Korea, it is hard to explain why Cm2 and Cm3 are more similar in mitochondrial genes. It can only speculate that there may be significantly different individuals in the natural population of C. medinalis.

Table 2.

Intraspecific and interspecific genetic distances of mitochondrial genes in Cnaphalocrocis exigua (Ce) and C. medinalis (Cm) (%)

Sequence type	Compared group	nad2	cox1	cox2	atp8	atp6	cox3	nad3	nad5	nad4	nad4l	nad6	cytb	nad1	rrnS	rrnL	22tRNA
Nucleotide	Ce-Cm1	0.0982	0.0969	0.0696	0.1973	0.0817	0.0836	0.0562	0.0658	0.0714	0.0462	0.0789	0.0695	0.0980	0.0422	0.0949	0.0408
	Ce-Cm2	0.0748	0.0628	0.0664	0.1181	0.0817	0.0836	0.0562	0.0633	0.0624	0.0352	0.0852	0.0667	0.0694	0.0423	0.0758	0.0421
	Ce-Cm3	0.0748	0.0628	0.0664	0.1105	0.0834	0.0836	0.0562	0.0633	0.0624	0.0317	0.0810	0.0619	0.0694	0.0423	0.0799	0.0421
	Cm1-Cm2	0.0543	0.0403	0.0059	0.1066	0.0000	0.0000	0.0000	0.0035	0.0212	0.0281	0.0094	0.0276	0.0497	0.0000	0.0650	0.0179
	Cm1-Cm3	0.0543	0.0403	0.0059	0.1139	0.0015	0.0000	0.0000	0.0035	0.0212	0.0245	0.0056	0.0230	0.0497	0.0000	0.0666	0.0186
	Cm2-Cm3	0.0020	0.0000	0.0000	0.0062	0.0015	0.0000	0.0000	0.0000	0.0015	0.0034	0.0075	0.0044	0.0000	0.0000	0.0051	0.0020
Amino acid	Ce-Cm1	0.1365	0.0490	0.0352	0.3654	0.0714	0.0496	0.0684	0.0588	0.0336	0.0515	0.1067	0.0419	0.0769	—	—	—
	Ce-Cm2	0.0890	0.0098	0.0308	0.2500	0.0714	0.0496	0.0684	0.0519	0.0336	0.0412	0.1236	0.0393	0.0481	—	—	—
	Ce-Cm3	0.0890	0.0098	0.0308	0.2308	0.0759	0.0496	0.0684	0.0519	0.0538	0.0309	0.1124	0.0340	0.0481	—	—	—
	Cm1-Cm2	0.0772	0.0392	0.0044	0.2037	0.0000	0.0000	0.0000	0.0104	0.0247	0.0515	0.0281	0.0288	0.0449	—	—	—
	Cm1-Cm3	0.0772	0.0392	0.0044	0.2222	0.0045	0.0000	0.0000	0.0104	0.0247	0.0412	0.0169	0.0236	0.0449	—	—	—
	Cm2-Cm3	0.0000	0.0000	0.0000	0.0185	0.0045	0.0000	0.0000	0.0000	0.0000	0.0103	0.0225	0.0052	0.0000	—	—	—

Cm1 indicates GenBank JN246082; Cm2 indicates GenBank JQ305693; Cm3 indicates GenBank JQ647917.

We also did evolutionary force analysis of each PCG of C. exigua and C. medinalis through calculating the rate of nonsynonymous (Ka), the rate of synonymous (Ks), and the ratio of Ka/Ks, respectively. In the comparison of the two species, the values of Ka/Ks of PCGs were all less than 1, which indicated that these genes were evolving under purifying selection. However, the ratio of Ka/Ks of atp8 gene was much higher than other PCGs (Supp Table 4 [online only]). Notably, the atp8 gene sequence of Cm1 showed too different from that of Cm2 and Cm3, there was only 1 amino acid distinction between Cm2 and Cm3, but exist 13 and 14 amino acid variations between Cm1–Cm2 pair and Cm1–Cm3 pair, respectively. Furthermore, the ratios of Ka/Ks for atp8 in Cm1–Cm2 pair (5.58) and Cm1–Cm3 pair (5.943) were even much higher than that of Ce–Cm pairs (0.993 for Ce–Cm1 pair, 0.814 for Ce–Cm2 pair, and 0.745 for Ce–Cm3 pair), and it implied that this gene might had suffered some kind of positive selection in partial populations of C. medinalis.

Since the atp8 gene exhibits significantly different characteristics, we further compared the atp8 gene of 38 species of Pyraloidea. As the result, the atp8 gene showed polymorphic in size not only within Cnaphalocrocis, but also within superfamily Pyraloidea (ranged from 156 to 168 bp) (Supp Fig. 1 [online only]). The A+T content of atp8 gene of Cm1 (93.9%), Cm2 (93.2%), and Cm3 (93.2%) was higher than C. exigua (90.4%) and the average level (91.0%) of 38 Pyraloidea species. Furthermore, Ile and Asn was the most used amino acid in atp8 of C. exigua and C. medinalis mitogenomes, respectively. Interestingly, the amino acid composition analysis show that the cysteine (Cys, C) and Glutamic acid (Glu, E) were specifically existed in atp8 of family Crambidae (Supp Fig. 1 [online only]).

As an inner membrane polypeptide of the F0 component, atp8 is responsible for the correct assembly of ATP synthase holoenzyme (Tzagoloff et al. 2004). Previous studies had found that genetic diversity of atp8 is associated with high-altitude adaptation in yak (Wang et al. 2017). At present, it is unclear whether the atp8 gene plays a key role in insect migration habits, and different individuals of C. medinalis and C. exigua have different energy production capabilities. Given that the atp8 gene plays an important role in energy metabolism, we believe that further investigation and research are still needed in this regard.

The tRNA Genes

In general, some tRNA genes showed highly conservative properties between these two species and between different C. medinalis individuals, while some tRNA genes showed certain variability. For instance, the trnM and trnW genes were not only conserved in different individuals of C. medinalis but also conserved in both leaffolders. Of the three individuals of C. medinalis, only 10 tRNA genes were identical and the conservative property of 22 tRNAs in the Cm2–Cm3 pair was higher than that of other pairs, which corresponds to Cm2 and Cm3 having 19 identical tRNA genes (Supp Table 3 [online only]). There were 11 and 12 differential tRNAs with base variation between Cm1 and Cm2, Cm1, and Cm3, respectively, and the differential bases were mainly found in the loop structures of the dihydrouridine (DHU) arm and the TΨC arm of tRNA. This also indicated that the variation of tRNA in the same species was mainly derived from the variation of the bases on the loop. The variation of tRNA in C. medinalis also indicated that there may be a large intraspecific genetic divergence in the natural populations of this species.

It was worth mentioning that the anticondon of trnS1 (AGN) were ‘TCT’ and ‘GCT’ in C. exigua and C. medinalis mitogenomes, respectively. According to the known mitogenomes, most pyralids contained a ‘GCT’ anticodon of trnS1, except for Paracymoriza distinctalis (TCT) and Eudonia angustea (ACT) (Timmermans et al. 2014, Ye and You 2016, Yang et al. 2018, Zhu et al. 2018). Furthermore, 11 other known pyralids had the same trnS1 gene sequence with that of C. medinalis and 8 of them belonged to subfamily Spilomelinae of Crambidae, thus implying the trnS1 gene may be conserved in more unreported moth species.

Overlap and Intergenic Regions

The overlaps of C. exigua and C. medinalis were in seven regions (24 bp in total) and six regions (22 bp in Cm1, 23 bp in Cm2 and Cm3), respectively (Table 1). In both leaffolders, the largest overlap was between trnW and trnC genes and contained a same motif ‘AAGCCTTA,’ the second large overlap was between atp8 and atp6 genes and shared a commom motif ‘ATGATAA’.

The mitogenomes of C. exigua contained 132 bp of intergenic spacer sequences, whereas Cm1, Cm2, and Cm3 had 216 bp, 223 bp, and 239 bp of intergenic spacer sequences, respectively (Table 1). There were four major intergenic spacers in C. exigua and C. medinalis mitogenomes (Fig. 3). The spacer between trnS2 and nad1 contained a conserved motif ‘ATACTAW,’ which was consistent with other Pyralids (Yang et al. 2018, Zhu et al. 2018). Interestingly, spacers between trnY and cox1 genes, nad3 and trnA genes were much longer in C. medinalis than in C. exigua, which were featured by AT repeats (Fig. 3). It was worth mentioning that this feature rarely appeared in the superfamily Pyraloidea, especially the former only found in the C. medinalis.

Fig. 3.

Comparison of the intergenic spacer sequences in the mitogenome of C. exigua (Ce) and C. medinalis (Cm). Cm1 indicates GenBank JN246082; Cm2 indicates GenBank JQ305693; Cm3 indicates GenBank JQ647917.

The A+T-rich Region

The A+T-rich region of C. medinalis also had a slightly higher A+T content than that of C. exigua (95.9 vs 94.4%) (Supp Table 2 [online only]). There was a conserved structure that consisted of the motif ‘ATAGW’ and a poly-T downstream of the rrnS gene in most known lepidopteran mitogenomes (Chai et al. 2012, Wan et al. 2013, Zhu et al. 2018). However, we found another style in the A+T-rich region of C. exigua at this site: the tetranucleotide ‘CTAG’ followed by poly-T sets. A different pattern occurred in C. medinalis, where the sequence ‘ATAG’ followed by a poly-T stretch (Fig. 4). Furthermore, a microsatellite (AT)9 was observed in C. exigua, while a microsatellite-like (AT)12 region preceded by motif ‘ATTTA’ was found in the 3′ end of C. medinalis (Cm2, Cm3) control region. A former study on 187 specimens of C. medinalis has found that the A+T-rich region of C. medinalis was 339–348 bp in length, 95.1–96.0.% in A+T content, and characterized by a conserved motif ‘ATAG,’ a poly-T stretch (10–16 bp), a microsatellite-like AT repeat (10–14 bp), and a 5-bp long-motif ‘ATTTA’ (Wan et al. 2011). Since C. exigua and C. medinalis show distinctions in the structure of A+T-rich region, the function of these conserved elements needs to be further studied.

Fig. 4.

Structures of the A+T-rich regions of C. exigua (Ce) and C. medinalis (Cm) and the other Pyraloidea species.

As a whole, it was important to note that these two species of Cnaphalocrocis were very close morphologically and hard to distinguish in the field. Therefore, it was surprising to find out that their mitochondrial genomes were so divergent.

Phylogenetic trees of Pyraloidea were reconstructed based on 13 PCGs of mitogenomes from 38 species with two species of Tortricoidea as outgroup. The results strongly supported the close sister relationship between C. exigua and C. medinalis, and both of them belonged to the subfamily Spilomelinae (Crambidae). To some extent, this results also provided an evidence to explain why these two major rice leaffolders were not only similar in morphological feature but also similar in damage characteristics of rice.

Both BI and ML phylogenetic analyses strongly supported the division of the superfamily Pyraloidea into two families (Fig. 5), Crambidae and Pyralidae. In the family Pyralidae, relationships between the four subfamilies received strongly supported and were consistent with previous research results (Regier et al. 2012, Yang et al. 2018, Zhu et al. 2018). In the family Crambidae, the sister subfamilies Pyraustinae and Spilomelinae was well supported. However, there are still some doubts about the genetic relationship of other subfamilies in this families. Considering that some subfamilies do not have any species with mitogenome being sequenced, and some subfamilies have only sequenced one or two species, we suggest that in order to improve or even solve these problems, more mitochondrial data are needed.

Fig. 5.

Phylogenetic tree of Pyraloidea was constructed using Bayesian inference and Maximum Likelihood analysis based on the 13 PCGs data. Values at each node indicated Bayesian posterior probabilities and the bootstrap percentages of maximum likelihood method.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.