Comparative mitogenomic analysis of mirid bugs (Hemiptera: Miridae) and evaluation of potential DNA barcoding markers

Sample and DNA extraction

Adult specimens of five mirid bugs were collected from alfalfa field in Shishe Town, Xifeng District, Qingyang City, Gansu Province, China, in July 2013. Samples and voucher specimens have been deposited in the State Key Laboratory of Grassland Agro-Ecosystems, College of Pastoral Agricultural Science and Technology, Lanzhou University, Lanzhou, China. All specimens were initially preserved in 100% ethanol in the field, and transferred to −20 °C until used for DNA extraction. The total genomic DNA was extracted from thorax muscle of a single specimen using the Omega Insect DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s protocols.

PCR amplification and sequencing

For each mirid species, mitogenome was amplified with 10–13 overlapping fragments by using universal insect mitogenome primers (Simon et al., 2006) and species-specific primers designed from sequenced fragments. All primers used in this study are provided in Table S1 . PCR and sequence reactions were conducted following our previous studies (Yuan et al., 2015a; Yuan et al., 2015b).

Genome annotation and sequence analysis

Sequence files were assembled into contigs with BioEdit 7.0.9.0 (Hall, 1999). Each of the five mirid mitogenomes newly sequenced in the present study was annotated following our previous studies (Yuan et al., 2015a; Yuan et al., 2015b). Nucleotide composition and codon usage were analyzed with MEGA 6.06 (Tamura et al., 2013). For each protein-coding gene (PCG) of all 15 mirid mitogenomes, the number of synonymous substitutions per synonymous site (Ks) and the number of nonsynonymous substitutions per nonsynonymous site (Ka) were calculated with MEGA 6.06 (Tamura et al., 2013). We also calculated the genetic distances for 13 PCGs and two ribosomal RNA genes (rRNAs) with MEGA 6.06 (Tamura et al., 2013) under the Kimura-2-parameter model (K2P). Strand asymmetry was calculated using the formulas: AT-skew = [A−T]/[A+T] and GC-skew = [G−C]/[G+C] (Perna & Kocher, 1995). To determine whether the Ka/Ks values of cox1-barcode sequences were relevant with the scope of sequences used, we downloaded all 7,759 cox1 sequences of Miridae available in GenBank (March 1, 2017). After removing sequences shorter than 658 bp, a total of 2,326 sequences in 144 genera were obtained (Table S7). Except for forty genera with only one sequence, the remaining 104 genera were used to calculate the values of Ka, Ks and Ka/Ks.

Phylogenetic analysis

Complete or nearly complete mitogenome sequences of eleven mirid bugs (15 samples, Table 1) and two outgroups from Pentatomomorpha (Corizus tetraspilus and Eurydema gebleri) (Yuan et al., 2015a; Yuan et al., 2015b) were used to perform phylogenetic analyses. The complete sequences of 13 PCGs (excluding stop codons), two rRNAs and 22 transfer RNA genes (tRNAs) were used for phylogenetic analyses. Each PCG was aligned individually based on codon-based multiple alignments by using the MAFFT algorithm within the TranslatorX (Abascal, Zardoya & Telford, 2010) online platform. Sequences of each rRNA gene were individually aligned using the MAFFT v7.0 online server with G-INS-i strategy (Katoh & Standley, 2013). Alignments of individual genes were then concatenated as a combined matrix with DAMBE 5.3.74 (Xia, 2013). Two datasets were assembled for phylogenetic analyses: (1) nucleotide sequences of 13 PCGs (P123) with 11,133 residues and (2) nucleotide sequences of 13 PCGs, two rRNAs and 22 tRNAs (P123RT) with 14,905 residues. To evaluate phylogenetic potential of single mitochondrial genes, i.e., each of 13 PCGs, rrnL and rrnS, as well as the combined 22 tRNAs, were also used in phylogenetic analyses.

The optimal partitioning schemes and corresponding nucleotide substitution models for each datasets were determined by PartitionFinder v1.1.1 (Lanfear et al., 2012). We created input configuration files that contained pre-defined data blocks by genes and codons, e.g., 39 partitions for P123, 42 partitions for P123RT and 3 partitions for each of 13 PCGs. The “greedy” algorithm with branch lengths estimated as “unlinked” and Bayesian information criterion (BIC) were used to search for the best-fit scheme (Table S2). The best-fit partitioning schemes selected by PartitionFinder were used in all subsequent phylogenetic analyses. We used jModelTest 2.1.7 (Posada, 2008) to determine the best evolutionary model for rrnL, rrnS and the combined 22 tRNAs.

Phylogenetic analyses were conducted with Bayesian inference (BI) and maximum likelihood (ML) methods available on the CIPRES Science Gateway v3.3 (Miller, Pfeiffer & Schwartz, 2010). Bayesian analyses were performed with MrBayes 3.2.3 (Ronquist et al., 2012) on Extreme Science and Engineering Discovery Environment (XSEDE 8.0.24). Two independent runs with four chains (three heated and one cold) each were conducted simultaneously for 1 × 10⁶ generations. Each run was sampled every 100 generations. Stationarity is assumed to be reached when ESS (estimated sample size) value is above 100 and PSRF (potential scale reduction factor) approach 1.0 as suggested in MrBayes 3.2.3 manual (Ronquist et al., 2012). The first 25% samples were discarded as burn-in, and the remaining trees were used to calculate posterior probabilities (PP) in a 50% majority-rule consensus tree. ML analyses were carried out using RAxML-HPC2 (Stamatakis, 2014) on XSEDE 8.0.24 with the GTRGAMMA model, and the node reliability was assessed by 1,000 bootstraps (BS).

General features of mirid mitogenomes

In the present study, we sequenced and annotated the mitogenomes of five mirid bugs: two were completely sequenced, whereas three were nearly complete mitogenomes (lacking sequences of three tRNAs and the putative control region) (Table 1, Table S3). The mitogenome sequences of five mirids have been deposited in GenBank of NCBI underaccession numbers: KU234536 –KU234540. The two completely sequenced mitogenomes contained 37 typical mitochondrial genes (i.e., 13 PCGs, 22 tRNA genes and two rRNAs) and a large non-coding region (putative control region) (Table S3). The order and orientation of the mitochondrial genes were identical to that of the putative ancestral insect mitogenome (Boore, 1999; Cameron, 2014). Gene overlaps and spacers were presented in several conserved positions in the mirid mitogenomes, e.g., trnS2/nad1 (7 bp), trnW/trnC (−8 bp), atp8/atp6 (−7 bp) and nad4/nad4L (−7 bp).

The tRNAs in the three Adelphocoris species could be folded into a classical clover-leaf secondary structure (Fig. 1). However, trnS1 (AGN) in Apo. lucorum and L. pratenszs species lacked the DHU stem-loop structures, as previously observed in many other true bugs (Wang et al., 2014b; Yuan et al., 2015a; Yuan et al., 2015b). All 22 tRNAs in Apolygus and Lygus species used the standard anticodon, whereas two tRNAs in the three Adelphocoris species were exceptions: trnS1 was predicted to have anticodon UCU, whereas trnK had the anticodon UUU (Table S3). The sequences and structures of anticodon arms and aminoacyl acceptor stems were well conserved within Miridae, whereas most of the variations (nucleotide substitutions and indels) were found in the DHU loops, pseudouridine (TψC) arms and variable loops (Fig. 1).

All of the mirid mitogenomes showed similar nucleotide composition in the J-strand: high A+T content, positive AT- and negative GC-skews (Fig. S1), as is usually observed in insect mitogenomes (Hassanin, Leger & Deutsch, 2005). For each part of mitogenomes, the A+T content, AT- and GC-skews showed low variability among different taxonomic levels (i.e., family, subfamily, genus and species). For AT-skew, a negative value was found in PCGs, rrnL, rrnS and the 2nd codon position of PCGs, whereas the 1st and 3rd codon positions of PCGs had positive AT-skew values. Except for Apo. lucorum, a positive AT-skew value also was found for tRNAs in all mirids. For GC-skew, the 1st codon position of PCGs, rrnL, rrnS and tRNAs were markedly positive, whereas the 2nd and 3rd codon positions of PCGs showed negative values. The pattern of codon usage in all analyzed mirid mitogenomes was consistent with previous findings in insects (e.g., Wang et al., 2015; Yuan et al., 2015a), namely that A+T-rich codons were preferably used (Table S4).

Gene variability in mirid mitogenomes

Among 13 PCGs, four genes (cox1, cox3, nad1 and nad3) had no length variability in the 15 examined mirid mitogenomes (Table S5). Four genes showed size variation only in one species (i.e., atp6, cox2 and nad4 in N. tenuis, cob in L. lineolaris). N. tenuis belonging to Bryocorinae showed the most length differences with species from Mirinae. The length of nad5 was most variable, but conserved within each genus. No length variation was found in the same species, whereas the most variations were found among genera. For the two genera including more than two species, no size variation was present in Adelphocoris, and only cob in L. lineolaris showed size difference with the others in Lygus. For rrnL and rrnS, intra-generic and -specific size differences were slight, whereas large differences were found among genera.

The numbers of variable and informative sites varied among genes and taxa (Table 2). The smallest gene atp8 showed the highest informative sites in Miridae, Mirinae and Mirini, followed by nad2 and nad6, while nad4L in Adelphocoris and nad3 in Lygus. The cox1 gene showed least informative sites in Miridae, Mirinae and Mirini, while nad6 in Adelphocoris and nad4L in Lygus. Generally, rrnL and rrnS had the lower informative sites than most PCGs; rrnL was slightly higher than rrnS, except for Adelphocoris. Compared to the whole cox1, cox1-barcode sequences contained relatively small informative sites in most taxonomic levels (except for Lygus). Except for cox1, other two longest genes (nad4 and nad5) contained moderate informative sites in most taxonomic levels (except for Lygus).

As expected, the largest genetic distances were found among subfamilies, followed by Tribes and genera, whereas the smallest among species (Fig. 2). Some genes (e.g., atp6, cox2, nad1, nad3 and nad6) showed no intraspecific variations at nucleotide and/or amino acid levels, indicating that these genes were highly conserved. The values of K2P and Ka in 13 PCGs of Lygus were larger than that of Adelphocoris, indicating that Lygus had the higher substitution rates. For the K2P distance, nad6 was the highest in Miridae and Lygus, whereas atp8 in Mirinae and Mirini. The K2P distance of nad4L was the highest within Adelphocoris, and relatively high K2P distance was also found in Miridae, Mirinae, Mirini and Lygus. The two rRNAs showed similar K2P distance with those of cox1-3 and cob in Miridae, Mirinae and Mirini, but had the lowest substitution rates in Lygus. For Miridae, Mirinae and Mirini, atp8 showed the highest Ka value, whereas nad4L in Adelphocoris and nad6 in Lygus showed the highest Ka values. Two genes (nad4 and nad5) showed relatively high substitution rates at the nucleotide and amino acid levels. Cox1 showed the lowest substitution rate in Miridae, Mirinae, Mirini, Adelphocoris and Lygus. The Ka/Ks values for all PCGs were far lower than 1 (<0.59) (Fig. 2), suggesting that these genes were evolving under purifying selection. However, we found that the Ka/Ks values for cox1-barcode sequences of the 15 mirids were always larger than 1 (1.34 in Ade. suturalis to 15.20 in Mirini; Table S6). In addition, analyses for cox1-barcode sequences of Miridae from GenBank showed that forty-four genera had no Ka and/or Ks values (Table S7), whereas the values of Ka/Ks for the remaining 61 genera were larger than 1 (1.5–263.0, Table S7), indicating that these sequences may be under positive selection or selection relaxation during the evolutionary process of Miridae. These results suggested that different genes had different substitution rates, whereas the same genes in different taxonomic levels showed large differences.

Mirid phylogeny based on combined mitochondrial genes

Phylogenetic relationships of 15 mirid mitogenome sequences were inferred using BI and ML methods based on two mitogenomic datasets (P123 and P123RT). The results showed that the two methods with the same dataset resulted in identical tree topology, whereas slight difference was found in the relationships of three species between the two datasets (Fig. 3). For the P123 dataset, Ade. nigritylus was sister to Ade. suturalis (PP = 0.65, BS = 84). For the P123RT dataset, Ade. nigritylus had a closer relationship with Ade. fasciaticolli and Ade. lineolatus (PP = 0.64, BS = 30), which was consistent with previous mirid phylogeny based on mitogenomic data (Wang et al., 2014b). However, low supports were present in both analyses (Fig. 3), indicating that the phylogenetic position of Ade. nigritylus was unstable. All analyses consistently supported the relationship of Nesidiocoris + (Trigonotylus + (Adelphocoris + (Apolygus + Lygus))), as previous mitogenomic analyses (Wang et al., 2014b). For Lygus, the two datasets consistently recovered a phylogeny of (rugulipennis + (lineolaris + (hesperus + pratenszs))). The monophyly of Adelphocoris was strongly supported by all analyses with high supports (PP = 1.0, BS = 100), as was the monophyly of Lygus (PP = 1.0, BS = 100).

Mirid phylogeny based on single mitochondrial gene

We performed phylogenetic analyses using BI and ML methods with single mitochondrial genes, including each of 13 PCGs, rrnL, rrnS and the combined 22 tRNAs (Figs. S2 and S3). The results showed that the tree topologies were variable among different datasets, indicating incongruent phylogenetic signals among genes, as reported in previous similar studies (Duchêne et al., 2011; Havird & Santos, 2014; Nadimi, Daubois & Hijri, 2016). However, all analyses consistently supported the monophyly of each of Adelphocoris (PP = 0.87–1.0, BS = 82–100) and Lygus (PP = 0.74–1.0, BS = 59–100), and several relationships within each of these two genera were recovered by different individual genes (Table 3, Figs. S2 and S3). Four datasets (nad4, nad5, rrnL and tRNAs) with the two analytical methods supported the phylogeny: Nesidiocoris + (Trigonotylus + (Adelphocoris + (Apolygus + Lygus))), as recovered by the P123 and P123RT datasets. This relationship among the five genera was also supported by BI analyses with two other genes (cox1 and nad2). These four datasets supported the Lygus phylogeny of (rugulipennis + (lineolaris + (hesperus + pratenszs))), as four other genes (cob, cox3, nad1 and nad6). For Adelphocoris, BI and ML analyses of four genes (cox2, cox3, nad5 and rrnL) supported (fasciaticollis + lineolatus) + nigritylus (PP = 0.85–1.0, BS = 34–98), whereas other four genes (atp6, cox1, nad1 and nad4) recovered the sister relationship of nigritylus and suturalis (PP = 0.63–0.9, BS = 53–97). It was notable that the cox1-barcode sequences did not support the sister relationship of Apolygus + Lygus consistently recovered by many independent datasets (Fig. 3 and Figs. S2–S4).