Conserved gene arrangement in the mitochondrial genomes of barklouse families Stenopsocidae and Psocidae

Abstract

Substantial variation in gene organization and arrangement has been reported for sequenced mitochondrial (mt) genomes from the suborders of the insect order Psocoptera. In this study we sequenced the complete mt genome of Stenopsocus immaculatus, the first representative of the family Stenopsocidae from the suborder Psocomorpha. Relative to the ancestral pattern, rearrangements of a protein-coding gene (nad3) and five tRNA genes (trnQ, trnC, trnN, trnS1, trnE) were found. This pattern was similar to that of two barklice from the family Psocidae, with the exception of the translocation of trnS1, trnE and trnI. Based on comparisons of pairwise breakpoint distances of gene rearrangements, gene number and chromosome number, it was concluded that mt genomes of Stenopsocidae and Psocidae share a relatively conserved pattern of gene rearrangements; mt genomes within the Psocomorpha have been generally stable over long evolutionary history; and mt gene rearrangement has been substantially faster in the booklice (suborder Troctomorpha) than in the barklice (suborders Trogiomorpha and Psocomorpha). It is speculated that the change of life history and persistence of unusual reproductive systems with maternal inheritance contributed to the contrasting rates in mt genome evolution between the barklice and booklice.

Keywords

gene rearrangement / mitochondrial genome / Psocoptera / Stenopsocidae / TDRL model

1 Introduction

The animal mitochondrial (mt) genome is important in cell metabolism, apoptosis, disease and aging^[1^]. The mt genome is the most extensively studied genomic system in insects, with sequenced representatives from all insect orders available in GenBank^[2^]. Sequence data from the mt genome has been used widely for phylogenetic analyses across broad taxonomic levels^[3–8^] due to the abundance of mt genomes in tissues, their greater rate of evolution than the nuclear genome and their evolutionarily conserved transcription products^[2,9,10^]

The insect mt genome is a compact circular molecule and typically encodes 13 protein-coding genes (PCGs), two rRNA (rRNA) genes and 22 tRNA (tRNA) genes. In addition to this set of 37 genes, there are a variety of non-coding regions, the largest of which is usually termed the control region (CR) as it contains both an origin of replication and transcription^[11^]. In general, the arrangement of genes in insect mt genomes is highly conserved^[2,12,13^], with notable exceptions in the Hymenoptera, Thysanoptera and Psocodea (Phthiraptera and Psocoptera)^[14–19^].

The insect order Psocoptera contains 39 extant families and more than 5000 described species in three suborders: Trogiomorpha (barklice), Psocomorpha (barklice) and Troctomorpha (booklice)^[20^]. To date, complete mt genomes have been published for three barklice and five booklice^[18,19,21–24^]. The mt genomes of the five booklice (all from the genus Liposcelis) differed significantly in gene arrangement from the generally-accepted ancestral insect mt genome, and even from each other^[18^], e.g., no commonalities were observed between the ancestral gene arrangement and the arrangement in Liposcelis entomophila^[21^], and only one commonality (atp8-atp6) was seen between the ancestral gene arrangement and the other four Liposcelis species^[22^]. Liposcelis decolor and Liposcelis sculptilis, like most other insects, have the typical single-chromosome mt genome; the other three Liposcelis species, however, possess multipartite mt genomes comprising two chromosomes^[22,23^]. The mt genomes of the barklice were much less rearranged than those of the booklice. The genome of the lepidopsocid (suborder Trogiomorpha) had rearrangements of eight genes, including a PCG^[24^]. Two genomes of the other two barklice, Psococerastis albimaculata and Longivalvus hyalospilus (family Psocidae, suborder Psocomorpha), had the least rearrangements, with the transposition of one PCG and five tRNA genes, explained by the mechanism of tandem duplication-random loss (TDRL)^[19^]. The variation observed within the currently-available barklice and booklice mt genomes greatly limits their use in reliably inferring evolutionary patterns within the Psocoptera, and additional mt genomic data from other families are needed.

To further understand the evolution of the mt genome organization in the Psocoptera, we sequenced the complete mt genome of Stenopsocus immaculatus, the first representative of the family Stenopsocidae. The conserved pattern of gene arrangement in the mt genomes of two barklice families, Stenopsocidae and Psocidae suggested that mt genomes within the suborder Psocomorpha have been stable over long evolutionary timescales after some early gene rearrangements.

2 Materials and methods

2.1 Samples collection and DNA extraction

Specimens of the narrow barklouse, S. immaculatus, were collected from Ashoro Washippu, Hokkaido, Japan in 2013 and kept in 100% ethanol before being transferred to -20°C for long-term preservation at the Entomological Museum of China Agricultural University, Beijing, China. Genomic DNA was extracted from six adult whole bodies excluding the abdomen with a DNeasy DNA Extraction kit (Qiagen, Valencia, CA, USA).

2.2 Mitochondrial genome sequencing, assembly and annotation

Two fragments of mtDNA (cox1 and rrnS) were amplified using PCR and Sanger sequencing following Li et al.^[25^] (Table S1). A library was prepared from the genomic DNA with an insert size of 450 bp and was sequenced on the Illumina Hiseq 2500 platform at Berry Genomics, Beijing. 4 Gb of clean data (250 bp paired-end reads) were used to assemble the mt genome using the method map to reference in Geneious v9.0.4^[26^] with the sequences of cox1 and rrnS as the reference. The assembly parameters were: minimum overlap identity 98%, no gaps, maximum mismatches per read 2%, maximum ambiguity 2, and minimum overlap 100 bp. PCGs and rRNAs were identified using BLAST searches in GenBank and subsequently by alignment with sequences of other barklice^[19,24^]. tRNA genes were identified by tRNAscan-SE Search Server v1.21^[27^] and checked manually. tRNA genes that could not be determined by tRNAscan-SE were determined in the unannotated regions by sequence similarity to tRNAs of other barklice. The annotated entire mt genome sequences of S. immaculatus have been deposited in GenBank under accession number KX187004.

2.3 Nucleotide composition, substitution rate and gene rearrangement analyses

The nucleotide composition of the mt genomes was calculated using MEGA 6.0^[28^]. The nonsynonymous (Ka) and synonymous (Ks) substitution rate was calculated using DnaSP v.5.0^[29^]. To compare the rates in mt gene rearrangement among major lineages of Psocoptera, pairwise breakpoint distances were calculated with the CREx web server^[30^].

2.4 Sequence alignment and phylogenetic analyses

Four barklice and five booklice from Psocoptera were included in phylogenetic analyses (Table 1). The true bug, Halyomorpha halys, was used as an outgroup. Sequences of 12 PCGs (without nad4L) along with two rRNA genes were used for the phylogenetic analyses. Each PCG was aligned separately based on multiple codon alignment with the MAFFT algorithm in the TranslatorX online platform^[32^]. The two rRNA genes were aligned with the MAFFT 7.0 online server^[33^].

The individual gene alignments were assembled into a concatenated data set which was partitioned by genes using both Bayesian inference (BI) and maximum likelihood (ML) analyses with MrBayes 3.23^[34^] and RAxML-HPC2 8.1.11^[35^]. The best-fitting model (GTR+ I+ G) for the nucleotide data set was determined with jModelTest 0.1.1^[36^]. For the ML analysis, the reliability of the inferred topology was assessed by performing 1000 rapid bootstrap replicates. For Bayesian analysis, two simultaneous runs of 10 million generations were performed for the data set and trees were sampled every 1000 generations, with the first 25% discarded as burn-in. The analysis was considered to have reached stationarity when the average standard deviation of split frequencies decreased to 0.01^[37^].

Tab.1 Species phylogenetically analyzed in this study

Order/Suborder	Family	Species	Accession number	Reference
Psocoptera
Psocomorpha	Psocidae	Psococerastis albimaculata	JQ910986	[19]
		Longivalvus hyalospilus	JQ910989	[19]
	Stenopsocidae	Stenopsocus immaculatus	KX187004	Present study
Trogiomorpha	Lepidopsocidae	Unidentified species	NC_004816	[24]
Troctomorpha	Liposcelidae	L. decolor	JX870621	[23]
		L. bostrychophila	JN645275-76	[22]
		L. paeta	NC_025505-06	[21]
		L. entomophila	NC_025503-04	[21]
		L. sculptilis	KX171073	[18]
Hemiptera
Heteroptera	Pentatomidae	Halyomorpha halys	NC_013272	[31]

3 Results and discussion

3.1 The mitochondrial genome of the narrow barklouse, Stenopsocus immaculatus

The mt genome of S. immaculatus was determined to be a double-strand circular DNA molecule 16991 bp in length, including the 37 typical coding genes (13 PCGs, 22 tRNAs and two rRNAs) and a 1720 bp CR (Fig. 1). Twenty-three genes are encoded on the major strand (J-strand) and other 14 genes are encoded on the minor strand (N-strand). In addition to the CR, another 14 non-coding regions were observed with the longest one (332 bp) flanked by trnI and trnM (Table S2). There were also 13 overlaps identified, ranging from 1 to 7 bp in length.

The nucleotide composition of the S. immaculatus mt genome had a bias toward A and T, with an A+ T content of 78.3%; the highest A+ T content was found in the third codon position of PCGs (88.9%) and the CR (84.8%) (Table 2). The whole mt genome was found to have a slight T skew and C skew, and the PCGs had an obvious T skew and moderate G skew. The codon usage of PCGs also contributed to the genome-wide bias toward AT and there was a strong bias to A or T at the third codon position of PCGs (Table S3). All of the PCGs had the standard ATN start codon. Eleven PCGs had the complete stop codon, either TAG or TAA, whereas cox3 and nad5 have an incomplete stop codon T, which is presumed to be subsequently completed by the post-transcriptional polyadenylation process^[38^].

The typical 22 tRNAs found in insect mt genomes were found in S. immaculatus, with lengths varying from 61 to 69 bp. Most of the tRNAs could be folded as classic clover-leaf structures, with the exception of trnS1, in which its dihydrouridine arm simply forms a loop, instead of a stem-loop (Fig. S1). Based on the secondary structure, a total of 22 unmatched base pairs were found in the tRNAs. Eighteen of them were G-U pairs, which form a weak bond. The remaining four pairs included one A-G, one C-U, and two A-C mismatches.

Fig.1 The mitochondrial genome of the narrow barklouse, Stenopsocus immaculatus. Arrows indicate the orientation of gene transcription. PCGs are shown as blue arrows, rRNA genes as purple arrows, tRNA genes as brown arrows and the control region as gray rectangle. Abbreviations of gene names are: atp6 and atp8 for ATP synthase subunits 6 and 8, cox1–3 for cytochrome oxidase subunits 1–3, cytb for Cytochrome b, nad1–6 and nad4L for NADH dehydrogenase subunits 1–6 and 4L, rrnL and rrnS for large and small rRNA subunits. tRNA genes are shown with their one-letter corresponding amino acids; the two tRNA genes for leucine and serine have different anticodons: L1 (TAG), L2 (TAA), S1 (TCT) and S2 (TGA). The GC content is plotted using a black sliding window, as the deviation from the average GC content of the entire sequence. GC-skew is plotted as the deviation from the average GC-skew of the entire sequence. The inner cycle indicates the location of genes in the mt genome.

Full size|PPT slide

The putative CR (1720 bp) was flanked by rrnS and trnI and is highly AT-rich (84.8%). Two tandem mirror repeats, unit I (683 bp) and unit II (686 bp), were identified in this region (Fig. 2). The sequence of mirror-repeat unit II was reversed with respect to unit I and showed high sequence similarity (99.4%) to unit I, with only four point mutations. Mirror repeats are known to form triplex H-DNA which are intrinsically mutagenic in mammalian cells^[39^]. They are known to induce large scale mutations such as deletions and/or rearrangements at a higher frequency^[40^]. Most of the repetitions found in the CRs of insect mt genomes, to date, are short repeat units similar to minisatellites with high copy numbers^[41^]. The occurrence of a large mirror repeats, with few mutations, found in the CR of the S. immaculatus mt genome is a rare event deserving the attention of future studies.

Fig.2 The control region of the narrow barklouse Stenopsocus immaculatus. I and II indicate two mirror-repeat units.

Full size|PPT slide

Tab.2 Nucleotide composition of the Stenopsocus immaculatus mitochondrial genome

Feature	T (U)	C	A	G	A%+ T%	AT-skew	GC-skew
Whole genome	39.9	11.8	38.4	9.8	78.3	- 0.019	- 0.093
Control region	42.8	7.5	41.9	7.7	84.8	- 0.011	0.015
Protein-coding genes	43.6	11.5	32.1	12.8	75.7	- 0.151	0.053
First codon position	35.0	10.7	35.7	18.4	70.9	0.006	0.263
Second codon position	46.0	18.3	21.3	14.3	67.4	- 0.367	- 0.123
Third codon position	49.0	5.5	39.5	5.6	88.9	- 0.112	0.017
tRNA genes	39.7	8.3	40.1	11.9	79.8	0.005	0.175
rRNA genes	43.5	6.7	39.0	10.8	82.5	- 0.055	0.234

Note: AT-skew= (A% - T%)/(A% + T%); GC-skew= (G% - C%)/(G% + C%).

3.2 Gene rearrangements in the mitochondrial genome of Stenopsocus immaculatus

Compared with the generally-accepted ancestral mt gene arrangement of insects, rearrangements of a PCG (nad3) and five tRNAs (trnM, trnC, trnN, trnS1 and trnE) were found in two regions of the S. immaculatus mt genome: (1) between cox3 and trnH, and (2) between CR and trnY (Fig. 3). These two regions are considered as the active regions for mt gene rearrangement in hymenopterans^[15^], flat bugs^[3^] and barklice from the family Psocidae^[19^].

Fig.3 Comparison of mitochondrial gene arrangement between Psocomorpha (Stenopsocidae and Psocidae) and the hypothetical ancestor of insects. Abbreviations of gene names follow Fig. 1. Genes are transcribed from left to right except those underlined, which have the opposite transcriptional orientation. Orange frames show two active regions of gene rearrangements.

Full size|PPT slide

Three TDRL events could explain the gene rearrangements observed in the mt genome of S. immaculatus. In the first hypothesized TDRL event, the gene cluster of trnI-trnQ-trnM-nad2-trnW-trnC-trnY was duplicated in tandem and one copy of each duplicated gene was randomly deleted, leading to the new gene arrangement trnI-trnM-trnC-trnQ-nad2-trnW-trnY. The gene rearrangement in the second active region could be derived from another hypothesized TDRL event (Fig. 4). The mt gene rearrangements in S. immaculatus is slightly more complex than those in two barklice form the Psocidae (suborder Psocomorpha) with a third hypothesized TDRL event changing trnS1-trnE to trnE-trnS1.

Fig.4 Inferred TDRL events that account for the mitochondrial gene rearrangements in the narrow barklouse Stenopsocus immaculatus. (a) Genes between CR and cox1; (b) genes between cox3 and nad4. Genes with crosses below were eliminated. Two longer non-coding sequences are highlighted in orange.

Full size|PPT slide

3.3 Evolution of mitochondrial genome organization in the Psocoptera

To investigate the evolution of mt genome organization in the Psocoptera, we inferred phylogenetic relationships between the four barklice and five booklice, representing all three suborders, for which complete mt genome sequences are available. From ML and BI analyses trees were inferred with an identical topology with strong support for most of the clades (Fig. 5). Within the suborder Psocomorpha, S. immaculatus (Stenopsocidae) was recovered as the sister to a clade comprising P. albimaculata and L. hyalospilus (Psocidae). The close relationship between these two Psocomorpha families was also supported by their highly similar gene arrangement, sharing the same gene boundaries, with the exception of the translocation of trnS1 and trnE and trnI. This result indicated a relatively conserved pattern of gene rearrangements in the mt genomes of Stenopsocidae and Psocidae. A phylogenetic study had previously suggested that these two families (Stenopsocidae and Psocidae) represent the most widely divergent clades within the suborder Psocomorpha^[42^]. This suggests that mt genomes within Psocomorpha have been stable over a long period of evolutionary history following some early gene rearrangements.

Fig.5 Phylogenetic relationships among major lineages of the Psocoptera inferred from mitochondrial genome sequences. Numbers close to the branching points are ML bootstrap support values (right) and Bayesian posterior probabilities (left) in percentages. Breakpoint distances are relative to the ancestral arranging pattern. Gene number in the bracket calculates gene without repetitions.

Full size|PPT slide

To compare the evolutionary rates of mt gene rearrangement among major lineages of the Psocoptera, the pairwise breakpoint distances of gene rearrangement, gene number and chromosome number, and nonsynonymous (Ka) and synonymous (Ks) substitution rates were calculated, and are given adjacent to the phylogenetic tree (Fig. 5). The mt genomes of the Psocoptera showed accelerated evolutionary rates with breakpoint distances of gene rearrangements, relative to the ancestral insect genome, ranging from 12 to 36. There were substantial differences found between the three suborders of the Psocoptera. The five booklice from the Troctomorpha have had much faster evolutionary rates than the three barklice from the other two suborders, with notably higher breakpoint distances and nucleotide substitution rates. The five booklice also showed more variation in terms of missing/copied genes and genome fragmentation (i.e., two chromosomes) (Fig. 5).

Differences in the mt genomes among the three suborders may be attributable to the obviously different ecologies of barklice and booklice. Barklice are entirely free-living insects, whereas booklice are closely associated with the nests of birds and mammals as well as human households. Our previous study indicated that ecological changes in booklice and parasitic lice (Phthiraptera) appears to be associated with an increased rate of mt gene rearrangement^[19^]. This finding is supported by a recent study which associated an extraordinarily divergent mitochondrial karyotype with maternally transmitted sex ratio distortion in Liposcelis bostrychophila^[43^]. When considering that parthenogenetic reproduction is common in booklice, we speculate that the change of life history and persistence of unusual reproductive systems with maternal inheritance may have contributed to the contrasting rates of evolution in mt genomes between the barklice and booklice.

4 Conclusions

We present the complete mt genome of  Stenopsocus immaculatus, the first representative of the family Stenopsocidae from the suborder Psocomorpha. Comparative mt genome analyses reveal conserved pattern of gene rearrangements in Stenopsocidae and Psocidae, and mt gene rearrangement has been substantially faster in the booklice (suborder Troctomorpha) than in the barklice (suborders Trogiomorpha and Psocomorpha). The change of life history and persistence of unusual reproductive systems with maternal inheritance is likely to contribute to the contrasting rates in mt genome evolution between the barklice and booklice.

References

[1] Boore J L. Animal mitochondrial genomes. Nucleic Acids Research, 1999, 27(8): 1767–1780

[2] Cameron S L. Insect mitochondrial genomics: implications for evolution and phylogeny. Annual Review of Entomology, 2014, 59(1): 95–117

[3] Song F, Li H, Shao R, Shi A, Bai X, Zheng X, Heiss E, Cai W. Rearrangement of mitochondrial tRNA genes in flat bugs (Hemiptera: Aradidae). Scientific Reports, 2016, 6(1): 25725

[4] Gillett C P, Crampton-Platt A, Timmermans M J, Jordal B H, Emerson B C, Vogler A P. Bulk de novo mitogenome assembly from pooled total DNA elucidates the phylogeny of weevils (Coleoptera: Curculionoidea). Molecular Biology and Evolution, 2014, 31(8): 2223–2237

[5] Li H, Shao R, Song N, Song F, Jiang P, Li Z, Cai W. Higher-level phylogeny of paraneopteran insects inferred from mitochondrial genome sequences. Scientific Reports, 2015, 5(1): 8527

[6] Song F, Li H, Jiang P, Zhou X, Liu J, Sun C, Vogler A P, Cai W. Capturing the phylogeny of Holometabola with mitochondrial genome data and Bayesian site-heterogeneous mixture models. Genome Biology and Evolution, 2016, 8(5): 1411–1426

[7] Tang M, Hardman C J, Ji Y, Meng G, Liu S, Tan M, Yang S, Moss E D, Wang J, Yang C, Bruce C, Nevard T, Potts S G, Zhou X, Yu D W. High-throughput monitoring of wild bee diversity and abundance via mitogenomics. Methods in Ecology and Evolution, 2015, 6(9): 1034–1043

[8] Ma C, Yang P, Jiang F, Chapuis M P, Shali Y, Sword G A, Kang L. Mitochondrial genomes reveal the global phylogeography and dispersal routes of the migratory locust. Molecular Ecology, 2012, 21(17): 4344–4358

[9] Simon C, Buckley T R, Frati F, Stewart J B, Beckenbach A T. Incorporating molecular evolution into phylogenetic analysis, and a new compilation of conserved polymerase chain reaction primers for animal mitochondrial DNA. Annual Review of Ecology Evolution and Systematics, 2006, 37(1): 545–579

[10] Crampton-Platt A, Yu D W, Zhou X, Vogler A P. Mitochondrial metagenomics: letting the genes out of the bottle. GigaScience, 2016, 5(1): 15

[11] Wolstenholme D R. Genetic novelties in mitochondrial genomes of multicellular animals. Current Opinion in Genetics & Development, 1992, 2(6): 918–925

[12] McBride H M, Neuspiel M, Wasiak S. Mitochondria: more than just a powerhouse. Current Biology, 2006, 16(14): R551–R560

[13] Clary D O, Wolstenholme D R. The mitochondrial DNA molecular of Drosophila yakuba: nucleotide sequence, gene organization, and genetic code. Journal of Molecular Evolution, 1985, 22(3): 252–271

[14] Shao R, Campbell N J, Barker S C. Numerous gene rearrangements in the mitochondrial genome of the wallaby louse, Heterodoxus macropus (Phthiraptera). Molecular Biology and Evolution, 2001, 18(5): 858–865

[15] Dowton M, Austin A D. Evolutionary dynamics of a mitochondrial rearrangement “hot spot” in the Hymenoptera. Molecular Biology and Evolution, 1999, 16(2): 298–309

[16] Shao R, Barker S C. The highly rearranged mitochondrial genome of the plague thrips, Thrips imaginis (Insecta: Thysanoptera): convergence of two novel gene boundaries and an extraordinary arrangement of rRNA genes. Molecular Biology and Evolution, 2003, 20(3): 362–370

[17] Shao R, Li H, Barker S C, Song S. The mitochondrial genome of the guanaco louse, Microthoracius praelongiceps: insights into the ancestral mitochondrial karyotype of sucking lice (Anoplura, Insecta). Genome Biology and Evolution, 2017, 9(2): 431–445

[18] Shi Y, Chu Q, Wei D D, Qiu Y J, Shang F, Dou W, Wang J J. The mitochondrial genome of booklouse, Liposcelis sculptilis (Psocoptera: Liposcelididae) and the evolutionary timescale of Liposcelis. Scientific Reports, 2016, 6(1): 30660

[19] Li H, Shao R, Song F, Zhou X, Yang Q, Li Z, Cai W. Mitochondrial genomes of two Barklice, Psococerastis albimaculata and Longivalvus hyalospilus (Psocoptera: Psocomorpha): contrasting rates in mitochondrial gene rearrangement between major lineages of Psocodea. PLoS One, 2013, 8(4): e61685

[20] Lienhard C, Smithers C N. Psocoptera: world catalogue and bibliography. Geneva: Muséum d’Histoire Naturelle Press, 2002.

[21] Chen S C, Wei D D, Shao R, Shi J X, Dou W, Wang J J. Evolution of multipartite mitochondrial genomes in the booklice of the genus Liposcelis (Psocoptera). BMC Genomics, 2014, 15(1): 861

[22] Wei D D, Shao R, Yuan M L, Dou W, Barker S C, Wang J J. The multipartite mitochondrial genome of Liposcelis bostrychophila: insights into the evolution of mitochondrial genomes in bilateral animals. PLoS One, 2012, 7(3): e33973

[23] Chen S C, Wei D D, Shao R, Dou W, Wang J J. The complete mitochondrial genome of the booklouse, Liposcelis decolor: insights into gene arrangement and genome organization within the genus Liposcelis. PLoS One, 2014, 9(3): e91902

[24] Shao R, Dowton M, Murrell A, Barker S C. Rates of gene rearrangement and nucleotide substitution are correlated in the mitochondrial genomes of insects. Molecular Biology and Evolution, 2003, 20(10): 1612–1619

[25] Li H, Liu H, Shi A, Štys P, Zhou X, Cai W. The complete mitochondrial genome and novel gene arrangement of the unique-headed bug Stenopirates sp. (Hemiptera: Enicocephalidae). PLoS One, 2012, 7(1): e29419

[26] Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics, 2012, 28(12): 1647–1649

[27] Lowe T M, Eddy S R. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Research, 1997, 25(5): 955–964

[28] Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Molecular Biology and Evolution, 2013, 30(12): 2725–2729

[29] Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics, 2009, 25(11): 1451–1452

[30] Bernt M, Merkle D, Ramsch K, Fritzsch G, Perseke M, Bernhard D, Schlegel M, Stadler P F, Middendorf M. CREx: inferring genomic rearrangements based on common intervals. Bioinformatics, 2007, 23(21): 2957–2958

[31] Lee W, Kang J, Jung C, Hoelmer K, Lee S H, Lee S. Complete mitochondrial genome of brown marmorated stink bug Halyomorpha halys (Hemiptera: Pentatomidae), and phylogenetic relationships of hemipteran suborders. Molecules and Cells, 2009, 28(3): 155–165

[32] Abascal F, Zardoya R, Telford M J. TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations. Nucleic Acids Research, 2010, 38(Web Server issue): W7–W13 

[33] Katoh K, Standley D M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular Biology and Evolution, 2013, 30(4): 772–780

[34] Ronquist F, Teslenko M, van der Mark P, Ayres D L, Darling A, Höhna S, Larget B, Liu L, Suchard M A, Huelsenbeck J P. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology, 2012, 61(3): 539–542

[35] Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics, 2006, 22(21): 2688–2690

[36] Posada D. jModelTest: phylogenetic model averaging. Molecular Biology and Evolution, 2008, 25(7): 1253–1256

[37] Huelsenbeck J P, Ronquist F, Nielsen R, Bollback J P. Bayesian inference of phylogeny and its impact on evolutionary biology. Science, 2001, 294(5550): 2310–2314

[38] Ojala D, Montoya J, Attardi G. tRNA punctuation model of RNA processing in human mitochondria. Nature, 1981, 290(5806): 470–474

[39] Wang G, Vasquez K M. Naturally occurring H-DNA-forming sequences are mutagenic in mammalian cells. Proceedings of the National Academy of Sciences of the United States of America, 2004, 101(37): 13448–13453

[40] Jain A, Wang G, Vasquez K M. DNA triple helices: biological consequences and therapeutic potential. Biochimie, 2008, 90(8): 1117–1130

[41] Zhang D X, Hewitt G M. Insect mitochondrial control region: a review of its structure, evolution and usefulness in evolutionary studies.  Biochemical  Systematics and  Ecology, 1997,  25(2): 99–120

[42] Yoshizawa K. Phylogeny and higher classification of suborder Psocomorpha (Insecta: Psocodea: ‘Psocoptera’). Zoological Journal of the Linnean Society, 2002, 136(3): 371–400

[43] Perlman S J, Hodson C N, Hamilton P T, Opit G P, Gowen B E. Maternal transmission, sex ratio distortion, and mitochondria. Proceedings of the National Academy of Sciences of the United States of America, 2015, 112(33): 10162–10168

Supplementary materials

The online version of this article at http://dx.doi.org/10.15302/J-FASE-2017158 contains supplementary materials (Tables S1–S3; Fig. S1).

Acknowledgements

This work was supported by grants from the National Natural Science Foundation of China (31401991, 31420103902, 31372229), the Beijing Natural Science Foundation (6152016, 6144027) and the Chinese Universities Scientific Fund (2017QC100, 2017QC066, 2017ZB002).

Compliance with ethics guidelines

Xiaochen Liu, Hu Li, Yao Cai, Fan Song, John-James Wilson, and Wanzhi Cai declare that they have no conflict of interest or financial conflicts to disclose.

This article does not contain any studies with human or animal subjects performed by any of the authors.

RIGHTS & PERMISSIONS

The Author(s) 2017. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)