with mitochondrial genome phylogeny

11 downloads 0 Views 2MB Size Report
Apr 7, 2016 - mitochondrial (mt) genome of Atylotus miser as well as the nearly complete ... subunits I-III; CytB, a gene for apocytochrome b; ND1–ND6 and.
Gene 586 (2016) 184–196

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Research paper

The complete mitochondrial genome of the Atylotus miser (Diptera: Tabanomorpha: Tabanidae), with mitochondrial genome phylogeny of lower Brachycera (Orthorrhapha) Kai Wang a,1, Xuankun Li a,1, Shuangmei Ding a, Ning Wang a,b, Meng Mao c, Mengqing Wang a,d,⁎, Ding Yang a,⁎⁎ a

Department of Entomology, China Agricultural University, Beijing, China Institute of Grassland Research, Chinese Academy of Agricultural Sciences, Hohhot, China c School of Biological Sciences, University of Wollongong, Wollongong, Australia d Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China b

a r t i c l e

i n f o

Article history: Received 22 December 2015 Received in revised form 22 March 2016 Accepted 5 April 2016 Available online 7 April 2016 Keywords: Mitochondrial genome Phylogeny Diptera Brachycera

a b s t r a c t Brachycera is a clade with over 80,000 described species and originated from the Mesozoic, and its larvae employ comprehensive feeding strategies. The phylogeny of the lower Brachycera has been studied intensively over the past decades. In order to supplement the lack of genetic data in this important group, we sequenced the complete mitochondrial (mt) genome of Atylotus miser as well as the nearly complete mt genomes of another 11 orthorrhaphous flies. The mt genome of A. miser is 15,858 bp, which is typical of Diptera, with 13 proteincoding genes, 22 tRNA genes, 2 rRNA genes and a 993 bp control region. The rest of the orthorrhaphous mt genomes in our study have the similar structure with A. miser. Additionally, we conducted a phylogenetic analysis of 20 mt genomes using Maximum-likelihood and Bayesian methods in order to reconstruct the evolutionary relationship of Orthorrhapha. The results show that all infraorders of Brachycera are monophyletic, and a relationship of Tabanomorpha + ((Xylophagomorpha + Stratiomyomorpha) + Muscomorpha) has been proposed. Within Xylophagomorpha, Nemestrinoidae forms the sister group of Xylophagidae. © 2016 Published by Elsevier B.V.

1. Introduction Mitochondrial (mt) genomes sequencing and analyses have developed rapidly in the past decades(Gibson et al., 2005; Rota-Stabelli et al., 2010; Bernt et al., 2013; Cameron, 2014b; Liu et al., 2014). Themt genome has been used for species identification (Ye et al., 2010; Krzywinski et al., 2011), molecular evolution, phylogenetic inference (Junqueira et al., 2004; Oliveira et al., 2008; Cameron et al., 2007; Moreno et al., 2010), population structure and phylogeography (Mousson et al., 2005; ScheVer and Lewis, 2006; Chatterjee et al., 2005; Nardi et al., 2010; Beckenbach, 2012), as well as genome structure and rearrangement (Matsumoto et al., 2009; Beckenbach and Joy, 2009; Behura et al., 2011; Yang et al., 2011; Torres et al., 2009; Nelson et al., 2011, 2012), because of its abundance in animal tissues, the small

Abbreviations: ATP6 and ATP8, genes for the ATPase subunits 6 and 8; COI–COIII, genes for cytochrome c oxidase subunits I-III; CytB, a gene for apocytochrome b; ND1–ND6 and ND4L, genes for NADH dehydrogenase subunits 1–6 and 4L;; rrnL, lrRNA, large (16S) rRNA subunit (gene); rrnS, srRNA, small (12S) rRNA subunit (gene). ⁎ Correspondence to: M. Wang, Department of Entomology, China Agricultural University, Beijing, China. ⁎⁎ Corresponding author. E-mail addresses: [email protected] (M. Wang), [email protected], [email protected] (D. Yang). 1 These authors contributed equally to this work.

http://dx.doi.org/10.1016/j.gene.2016.04.013 0378-1119/© 2016 Published by Elsevier B.V.

genome size, faster rate of evolution, low level or absence of sequence recombination, and evolutionary conserved gene products (Lin and Danforth, 2004; Simon et al., 2006; Gissi et al., 2008). Besides that, mt genomes can provide more genome-level characteristics than the shorter sequences of individual genes (Dowton et al., 2002; Boore et al., 2005; Masta and Boore, 2008; Boore, 2006a). As a major suborder of Diptera, Brachycera is a monophyletic group range over Nemestrinidae, Rhagionidae, Tabanidae, Xylomyidae, Stratiomyidae, Xylophagidae, Empididae, Dolichopodidae, Asilidae and Bombylidae etc. Brachycera is a clade with a Mesozoic origin and have over 80,000 described species. Well-preserved tabanids, nemestrinids, bombyliids and mydids have just been recovered from the Upper Jurassic (Ren, 1998). A number of species-rich families of the lower Brachycera diversified in the mid-Cretaceous, coincident with the radiation of angiosperms (Grimaldi, 1999). Most brachyceran larvae inhabit moist terrestrial habitats, and their adults are more stout bodied and compact than those of the lower Diptera. Orthorrhapha includes all Brachycera excluding Cyclorrhapha. However, as the group is paraphyletic, Orthorrhapha cannot be recognized formally, instead, it is a grouping of convenience only. The relationships within Brachycera have been the subject of numerous qualitative analyses using morphological characters. These analyses are often based on characters from one or a few character systems and the results on relationships have been controversial. More than five

K. Wang et al. / Gene 586 (2016) 184–196

different relationships among Brachycera have been proposed. Griffiths divided Brachycera into six major subgroups based on morphological characters, under the names Stratiomyomorpha, Vermileonomorpha, Tabanomorpha, Nemestrinoidea, Pleroneura, and Eremoneura. He also suggested Stratiomyomorpha as the sister-group of all the remaining Brachycera (Griffiths, 1994). Woodley divided Brachycera into four infraorders: Xylophagomorpha, Stratiomyomorpha, Tabanomorpha and Muscomorpha, and Muscomorpha was composed of Nemestrinoidea, Asiloidea, Empidoidea and Muscoidea (Fig. 1a) (Woodley, 1989). Yeates mostly agreed with Woodley and further suggested the relationships within Brachycera as (Stratiomyomorpha + (Xylophagomorpha + Tabanomorpha)) + (Nemestrinidae + (Acroceridae + (Asiloidea + (Empidoidea + Cyclorrhapha)))), but found Nemestrinoidea was a paraphyletic group (Fig. 1b) (Yeates, 2002). A recent morphological study supported a grade lower Brachycera as (Xylophagomorpha + (Stratiomyomorpha + (Tabanomorpha + (Acroceridae + (Asiloidea + (Empidoidea + Cyclorrhapha))))). (Fig. 1d) (Lambkin et al., 2013). Weigmann et al. detected a monophyletic Orthorrhapha based on morphological and molecular data (nuclear and mitochondrial genomes data), which was composed of (((Nemestrinidae + Xylophagomorpha) + Tabanomorpha) + ((Acroceridae + Stratiomyomorpha) + Asiloidea)), and the relationships of Brachycera were proposed as (Orthorrhapha + (Empidoidea + Cyclorrhapha)) (Fig. 1c) (Wiegmann et al., 2011). We used mt genome data to reconstruct the phylogenetic relationship of Orthorrhapha. Before this study, only two complete orthorrhaphous mt genome sequences (Cydistomyia duplonotata, NC_ 008756; Trichophthalma punctata, NC_008755) were available in GenBank (as of March 2015). Here, we sequenced the complete mt genome of Atylotus miser, as well as another 11 nearly-complete orthorrhaphous mt genomes (Table 1). We annotated these genomes using procedures and quality control methods proposed by Cameron (2014a). These new mt genomes will facilitate our understanding of the taxonomic position and evolutionary relationships of Orthorrhapha,

185

Table 1 Summary of mt genome sequences from Brachycera. Family

Species

GenBank accession Length Remark no. (bp)

Xylophagidae Xylomyidae Stratiomyidae Rhagionidae Athericidae Tabanidae Tabanidae Tabanidae Nemestrinidae Empididae Empididae Dolichopodidae Asilidae Asilidae Tipulidae Chironomidae Syrphidae Muscidae Drosophilidae Agromyzidae

Dialysis sp. Xylomya moiwana Ptecticus aurifer Rhagio sp. Suragina sp. Atylotus miser Cydistomyia duplonotata Chrysops silvifacies Trichophthalma punctata Rhamphomyia insignis Heleodromia immaculata Dolichopus bigeniculatus Leptogaster longicauda Satanas sp. Tipula abdominalis Chironomus tepperi Simosyrphus grandicornis Haematobia irritans Drosophila yakuba Liriomyza bryoniae

KT225293 KT225302 KT225297 KT225298 KT225301 KT225291 NC_008756 KT225292 NC_008755 KT225299 KT225295 KT225294 KT225296 KT225300 JN861743 NC_016167 NC_008754 NC_007102 NC_001322 NC_016713

14,465 14,681 14,754 14,481 14,599 15,858 16,247 14,544 16,396 14,636 14,515 14,466 14,407 14,415 14,566 15,652 16,141 16,078 16,019 16,183

–,* –,* –,* –,* –,* * –,* –,* –,* –,* –,* –,* *

Note: “–” incomplete data; “*” data we sequenced.

and help to select optimized primer for atypical regions in further molecular research of related taxa. Phylogenetic trees based on the mt genome data from Orthorrhapha were inferred by both Maximumlikelihood and Bayesian methods. 2. Materials and methods 2.1. Ethics statement No specific permits were required for the insects collected for this study. These specimens were collected by sweeping. The field studies

Fig. 1. Alternative hypotheses on relationships of major Brachycera. (a). Woodley (1989). (b). Yeates (2002). (c). Wiegmann et al. (2011). (d). Lambkin et al. (2013).

186

K. Wang et al. / Gene 586 (2016) 184–196

did not involve endangered or protected species. The species herein studied are not included in the “List of Protected Animals in China”. 2.2. Sampling and DNA extraction The specimen of A. miser used for DNA extraction was collected by Jingying Yang and Lei Zhang from Yunwu Mountain, Guyuan, Ningxia, China, in 3rd August 2013. Other specimens used for DNA extraction were also collected from China (S1 Table). After collection, all specimens were initially preserved in 95% ethanol in the field, and then transferred to −20 °C for the long-term storage upon the arrival at China Agricultural University. These specimens were examined and identified by the corresponding author Ding Yang with ZEISS Stemi 2000-c microscope. Whole genomic DNA was extracted from the thoracic muscle tissues using TIANamp Genomic DNA Kit (TIANGEN, Beijing, China). 2.3. PCR amplification and sequencing The mt genome of A. miser was amplified in 31 overlapping PCR fragments using NEB Long Taq DNA polymerase (New England Biolabs, Ipswich, MA). Firstly, twenty-five fragments were amplified using the universal primers (S2 Table). Then, six specifically designed primers based on the known sequences were used for the secondary PCRs (S2 Table). The same fragments amplification and sequencing method were applied to other specimens, except for the species-specific primers were shown in S2 Table. PCR cycling consisted these conditions: 95 °C for 30s, 40 cycles at 95 °C for 10 s, 42–55 °C (depending on the primers pair used) for 50 s, 65 °C for 1 kb/min (depending on the size of target amplicon) (S2 Table), and 65 °C for 10 min. The quality of PCR products was assessed through electrophoresis in a 1% agarose gel and stained with Gold View (ACME). All amplicons were sequenced in both directions using the BigDye Terminator Sequencing Kit (Applied Bio Systems, USA) and the ABI 3730XL Genetic Analyzer (PE Applied Biosystems, San Francisco, CA, USA) using amplification primers and internal primers developed by primer walking. 2.4. Bioinformatic analysis Sequences were proof-read and aligned into contigs in BioEdit version 7.0.5.3 (Hall, 1999). After fully sequencing the mt genome, it was annotated using both automated annotation methods and by hand. Quality control of the hand alignments was performed by comparing with homologous sequences from previously sequenced C. duplonotata (NC_008756) mt genomes to identify several tRNAs apparently, and to verify PCGs and rRNAs annotations (X.K. Li et al., 2015). Nucleotide substitution rates, base composition and codon usage were analyzed with MEGA 5.0 (Tamura et al., 2011). Nucleotide compositional skew was measured using the following formula: AT-skew = (A − T)/(A + T) (Perna and Kocher, 1995). 2.5. Phylogenetic analysis A total of 20 species of dipteran insects were used in phylogenetic analysis, including two outgroup species from Chironomidae and Tipulidae. Details of the species used in this study were listed in Table 1. Sequences of 13 PCGs, two rRNAs and 19 tRNAs were used in phylogenetic analysis (Because of the lack of tRNAIle,tRNAGln and tRNAMet in some species, we just use remaining 19 tRNAs). Each PCG was aligned individually based on Hand annotation method followed the procedures proposed by Cameron (2014a). The sequences of tRNAs and rRNAs were aligned respectively using MEGA 5.0 (Tamura et al., 2011), ambiguous positions in the alignment of RNAs were filtered by hand. Individual genes were concatenated using SequenceMatrix v1.7.8 (Vaidya et al., 2010). We assembled four datasets for

phylogenetic analysis: 1) sequences of 13 PCGs (PCG123) with 11,163 residues, 2) sequences of 13 PCGs, two rRNAs and 19 tRNAs (PCG123RNA) with 14,453 residues, 3) sequences of 13 PCGs exclude the third codon sites (PCG12) with 7442 residues and 4) sequences of 13 PCGs exclude the third codon sites, two rRNAs and 19 tRNAs (PCG12RNA) with 10,732 residues. We used PartitionFinder v1.1.1 (Lanfear et al., 2012) to select the optimal partition strategy and substitution models for each partition. We created an input configuration file with 63/29 (with/without 3rd codon positions) pre-defined partitions of the dataset, and used the “greedy” algorithm with branch lengths estimated as “unlinked” and Bayesian information criterion (BIC) to search for the best-fit scheme (S3 Table) (X.K. Li et al., 2015). Subsequently, we performed maximum likelihood (ML) and Bayesian inference (BI) using the best-fit partitioning schemes recommended by PartitionFinder (S3 Table). For ML analysis, we used RAxML 8.0.0 (Stamatakis, 2006) with 1000 bootstrap replicates and using the rapid bootstrap feature (random seed value 12,345) (Stamatakis et al., 2008). The Bayesian analysis was conducted with MrBayes 3.2.2 (Ronquist and Huelsenbeck, 2003). Two simultaneous runs of 2 million generations were conducted for the dataset, each with one cold and three heated chains. Samples were drawn every 1000 Markov chain Monte Carlo (MCMC) steps, with the first 25% discarded as burn-in. When the average standard deviation of split frequencies was below 0.01, we considered that the stationarity was reached and stopped run (X.K. Li et al., 2015). 3. Results and discussion 3.1. Genome organization and structure The complete mt genome of A. miser is a typical circular, doublestranded molecule with 15,858 bp in size (Fig. 2). The mt genome of A. miser contains 37 genes, including 13 PCGs, 22 tRNA genes, two rRNA genes and a large control region, which are usually present in bilaterian animals (Cameron, 2014a). The gene order is consistent with the inferred ancestral insect mt genome pattern, which is conserved among all brachyceran flies sequenced to date. Twenty-three genes are encoded on the majority strand (J-strand), while the minority strand (N-strand) encodes the remaining 14 genes. There are eight gene boundaries where sequence overlaps between neighboring genes, ranging from 1 to 8 bp in size, which together total 30 bp. The length of intergenic sequences (excluding the control region) is 1–16 bp found at 16 gene boundaries for a total of 90 bp (Table 2). The longest overlapping sequence is located between tRNATrp and tRNACys. While the longest intergenic sequence is located between tRNAGlu and tRNAPhe. Among all sequenced mt genomes of orthorrhaphous flies, the stable overlapping sequences are located between tRNATrp and tRNACys, ATP8 and ATP6, ND4 and ND4L. There are also 4 conserved intergenic sequences found in all twelve mt genomes in this study: tRNAGlu-tRNAPhe, ND5-tRNAHis, tRNASer(UCN)-ND1 and ND1- tRNALeu(CUN). The overlapping sequence between tRNATrp and tRNACys is usually 8 bp long (AAGCCT TA), which exists in all sequenced orthorrhaphous flies, but we found two exceptions in our sequenced data, Rhamphomyia insignis (24 bp intergenic sequences) and Suragina sp. (8 bp intergenic sequences). Generally, there is a 7 bp overlapping sequence between ATP8 and ATP6 in insects' mt genome (Lavrov et al., 2000). However, among all sequenced orthorrhaphous flies, ATP8–ATP6 consistently overlaps by 7 bp with a-1 frame shift (A TGA TAA → ATG ATA A) as found in all sequenced cyclorrhaphan flies (X.K. Li et al., 2015). ND4-ND4L always overlapped by 7 bp sequences (TTAACAT). tRNAGlu-tRNAPhe usually had a intergenic base about 20 bp. Interestingly, we found that Xylomya moiwana has 224 bp intergenic fragment here while Rhagio sp. does not have this intergenic base. Similarly, ND5-tRNAHis, tRNASer(UCN)-ND1 and ND1-tRNALeu(CUN) usually have intergenic base of about 15 bp, 15 bp, and 10 bp long respectively. And the region between tRNASer(UCN)

K. Wang et al. / Gene 586 (2016) 184–196

187

Fig. 2. Mitochondrial map of A. miser. Circular maps were drawn with CGView (Grant and Stothard, 2008). Arrows indicated the orientation of gene transcription. The tRNAs are denoted by the color blocks and are labelled according to the IUPACIUB single-letter amino acid codes (L1: CUN; L2: UUR; S1: AGN; S2: UCN). The GC content was plotted using a black sliding window, as the deviation from the average GC content of the entire sequence. GC-skew was plotted as the deviation from the average GC-skew of the entire sequence. The inner cycle indicated the location of genes in the mt genome.

and ND1 was supposed to contain transcription termination signal site (mtTERM) in the transcription process of mt genome (Cameron and Whiting, 2008; Sheffield et al., 2008).

3.2. Base composition and codon usage As with other insects, the nucleotide composition of the A. miser mt genome is biased toward A and T. The overall A + T content is 77.7%, nearly equal to that of C. duplonotata and slightly higher than that of T. punctata. A detailed comparison of nucleotide composition, ATskew, and GC-skew is given in Table 3. A. miser and C. duplonotata have weakly negative AT-skew, while all three species have positive GC-skew, which are consistent with the strand skew biases found in most metazoan mt genomes (weakly positive AT-skew and strongly negative GC-skew for the J-strand), except for T. punctata. In insects, the degree of AT-skew is related to gene direction, replication and

codon positions, whereas the degree of GC-skew is affected by reversals in replication orientation (Wei et al., 2010). The A + T bias is also reflected by relative codon usage of the PCGs. In the mt genome of A. miser, A + T rich codons, such as TTA (Leu), ATT (Ile), TTT (Phe), ATA (Met), AAT (Asn), are more frequently used than G + C rich codons. Among both strands, NNA-form codons are predominantly used while NNC-typecodons are rarely used. On J-strand, NNA is preponderant and NNG is least used with the reverse found in N-strand encoded PCGs (NNU commonest, NNC least common) (Fig. 3, S4 Table). Among all sequenced orthorrhaphous PCGs, the nucleotide composition is also biased toward A and T, which range from 63.69% (in Satanas sp.) to 77.78% (in Heleodromia immaculata) (S5 Table). A comparative analysis of nucleotide composition (A + T%, G + C%) versus skew (AT- and GC-skew) across the orthorrhaphous PCGs is shown in Fig. 4. All orthorrhaphous PCGs have almost equal G and C content (Fig. 4a). However, all Orthorrhapha have a negative GC-skew in J-

188

K. Wang et al. / Gene 586 (2016) 184–196

Table 2 Organization of the mt genome of A. miser. Gene

tRNAIle tRNAGln tRNAMet ND2 tRNATrp tRNACys tRNATyr COI tRNALeu(UUR) COII tRNALys tRNAAsp ATP8 ATP6 COIII tRNAGly ND3 tRNAAla tRNAArg tRNAAsn tRNASer(AGN) tRNAGlu tRNAPhe ND5 tRNAHis ND4 ND4L tRNAThr tRNAPro ND6 CytB tRNASer(UCN) ND1 tRNALeu(CUN) rrnL tRNAVal rrnS Control-region

Strand

Position

J N J J J N N J J J J J J J J J J J J J J J N N N N N J N J J J N N N N N

Size

From

To

1 67 144 214 1244 1305 1378 1443 2980 3047 3735 3809 3876 4031 4708 5499 5565 5922 5988 6054 6121 6187 6270 6338 8073 8140 9472 9771 9836 9904 10,433 11,568 11,651 12,601 12,666 13,999 14,071 14,866

66 135 213 1245 1312 1372 1444 2979 3045 3734 3805 3875 4037 4708 5496 5564 5918 5988 6051 6119 6185 6253 6337 8057 8139 9478 9768 9835 9901 10,428 11,569 11,635 12,590 12,665 13,998 14,070 14,865 15,858

Codon Start

66 69 70 1032 69 68 67 1537 66 688 71 67 162 678 789 66 354 67 64 66 65 67 68 1720 67 1339 297 65 66 525 1137 68 940 65 1333 72 795 993

Anticodon

Intergenic nucleotides

Stop GAT TTG CAT

ATT

0 8 0 −2 −8 5 −2 0 1 0 3 0 −7 −1 2 0 3 −1 2 1 1 16 0 15 0 −7 2 0 2 4 −2 15 10 0 0 0 0

TAA TCA GCA GTA

TCG

T

ATG

T

TAA CTT GTC ATT ATG ATG

TAA TAA TAA

ATT

TAA

TCC TGC TCG GTT GCT TTC GAA ATT

T

ATG ATG

T TAA

GTG

TGT TGG ATT ATG

TAA TAG

ATA

T

TGA TAG TAC

Note: Intergenic nucleotide, minus indicates overlapping between genes. tRNAX: where X is the abbreviation of the corresponding amino acid.

strand PCGs (Fig. 4b) and a strongly positive GC-skew in N-strand PCGs (Fig. 4c). As for AT-skew, all Orthorrhapha have consistently negative AT-skew except that Satanas sp. has weakly positive AT-skew in Jstrand PCGs (Fig. 4d, e, f), but the level of nucleotide skew in N-strand PCGs is higher than that of all PCGs and J-strand PCGs. 3.3. Protein-coding genes Thirteen protein-coding genes were identified in the mt genome of A. miser, with characteristics similar to those of other dipteran

species. The average A + T content across all protein-coding genes is 76.1%, slightly higher than that of C. duplonotata (75.5%) and T. punctata (72.0%). Table 3 shows the AT-skews and GC-skews of the three codon positions for A. miser, C. duplonotata and T. punctata. Among the three species, the A + T content of the third codon positions (94.4%, 93.1%, and 83.5% for A. miser, C. duplonotata and T. punctata, respectively) is higher than those of the first (68.5%, 67.8% and 66.6%) and second codon positions (65.4%, 65.7% and 65.8%). There is moderately a negative AT-skew for the PCGs as a whole (− 0.16) driven by strongly negative AT-skew at second

Table 3 Comparison of mitochondrial nucleotide composition in three species. A+T%

G+C%

AT-skew

GC-skew

Region Whole mitogenome Protein-coding genes 1st codon position 2nd codon position 3rd codon position tRNA genes lrRNA srRNA Control region

A.mis

C.dup

T.pun

A.mis

C.dup

T.pun

A.mis

C.dup

T.pun

A.mis

C.dup

T.pun

77.7 76.1 68.5 65.4 94.4 77.0 81.9 79.2 87.8

77.9 75.5 67.8 65.7 93.1 76.9 82.3 78.8 92.6

74.0 72.0 66.6 65.8 83.5 74.8 78.5 77.5 81.2

22.3 23.9 31.5 34.6 5.6 23.0 18.1 20.8 12.2

22.1 24.5 32.2 34.3 6.9 23.1 17.7 21.2 7.4

26.0 28.0 33.4 34.2 16.5 25.2 21.5 22.5 18.8

0.00 −0.16 −0.11 −0.39 −0.04 0.01 −0.02 −0.02 −0.03

0.00 −0.15 −0.10 −0.39 −0.02 0.00 −0.01 −0.01 0.03

0.09 −0.14 −0.05 −0.40 −0.02 0.02 −0.11 −0.12 −0.12

−0.17 0.04 0.30 −0.15 −0.19 −0.05 0.33 0.31 −0.31

−0.18 0.02 0.26 −0.16 −0.16 0.16 0.32 0.31 0.33

−0.24 −0.02 0.18 −0.13 −0.23 0.17 0.36 0.30 0.18

Note: A.mis indicates Atylotus miser, C.dup indicates Cydistomyia duplonotata and T.pun indicates Trichophthalma punctata. The A + T and G + C biases of protein-coding genes were calculated by AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C), respectively.

K. Wang et al. / Gene 586 (2016) 184–196

189

Fig. 3. Relative synonymous codon usage (RSCU) in the A. miser mt genome. Codon families are provided on the X-axis. Stop codon is not given. Red codon, codon not present in the chain/ genome.

codon positions (− 0.39), with weak skew at first, and third codon positions (− 0.11 and − 0.04 respectively). The absence of significant GC-skew across the PCGs as a whole (0.04) masks strong skews at each codon position with first codon positions strongly positive (0.30) balanced by strongly negative skews at second and third codon positions (− 0.15 and − 0.19 respectively) (Table 3). ATN, GTG (Val) and TTG (LeuUUR) are accepted as the canonical start codons for invertebrate PCGs (Cameron et al., 2012), while TCG (SerUCN) has been proposed to be an additional start codon commonly found in flies (Cameron et al., 2007). All 13 PCGs in A. miser mt genome use canonical start codons. ND2, ATP8, ND3, ND5 and ND6 start with ATT (Ile); ATA (Met) were found in ND1; CO1 starts with TCG (SerUCN); and ATG (Met) is used in the remaining 6 PCGs as start codons (Table 2).

Among the orthorrhaphous flies sequenced to date, ATG (Met) is the most frequently used start codon followed by ATT (Ile). ATG (Met) is almost exclusively used in the ATP6, CO2, CO3, CYTB, ND4 and ND4L genes of almost all orthorrhaphous species. T. punctata uses ATA (Met) for CO3 and ND4L of Satanas sp (Asilidae). starts with non-canonical start codon, AAA (Lys). ATT (Ile) is used for ND2, ATP8, ND3, ND5 and ND6 in most species. Furthermore, ATA (Met) and ATC (Ile) are also found in ND2, ATP8, ND3, ND5, ND6, CO3 and ND1 in minority of orthorrhaphous species. TTG (LeuUUR) is used in ND1 in most species. In addition, another three non-canonical start codons CCG (Pro), TTA (LeuUUR) and GTT (Val) have been proposed. CCG (Pro) was found as the start codon for CO1 in Rhagio sp. (Rhagionidae), Suragina sp. (Athericidae) and T. punctata. TTA (LeuUUR) was shown as the start codon for ND1 in Dialysis sp.

190 K. Wang et al. / Gene 586 (2016) 184–196 Fig. 4. AT% vs AT-skew and GC% vs GC-skew in orthorrhaphous PCGs. Measured in bp percentage (X-axis) and level of nucleotide skew (Y-axis). (a) and (d). Relationship between nucleotide content and skew in all PCGs. (b) and (e). Relationship between nucleotide content and skew in PCGs encoded on majority strand. (c) and (f). Relationship between nucleotide content and skew in PCGs encoded on minority strand.

K. Wang et al. / Gene 586 (2016) 184–196

(Xylophagidae), Suragina sp. (Athericidae), R. insignis (Empididae), and Satanas sp. (Asilidae), and it was also found in ND5 of Dialysis sp. GTT (Val) was proposed as the start codon for ND1 in Dolichopus bigeniculatus (Dolichopodidae) and Leptogaster longicauda (Asilidae). The non-canonical start codons are not newfangled in insects. For example, AAA (Lys), CCG (Pro) and TTA (LeuUUR) as non-canonical start codons were also found in other insects' PCGs (Fenn et al., 2008; Lee et al., 2006; Nardi et al., 2001; Shao and Barker, 2003). However, GTT (Val) has not been proposed as the start codon yet. Therefore, more mt genome data of insects is needed to verify its correctness. The stop codons most commonly used in A. miser are TAA (ATP6, ATP8, CO3, ND2, ND4L, ND6, ND3) and TAG (CYTB) found in 8 of the 13 PCGs. The remaining five PCGs (ND1, ND4, ND5, CO1, CO2) employ the partial stop codon T, which has been found in many insect mt genomes and is completed to a full TAA stop codon via posttranscriptional polyadenylation (Ojala et al., 1981). In all reported orthorrhaphous flies till now, TAA is the most common stop codon which is used in every PCG in at least one species, most notably in the ATP6, ATP8, CO3, ND4L and ND6 genes from all species sequenced. TAG has been found in CYTB, ND3, ND4 and ND5. The incomplete stop codons TA and T were found in CO1,CO2, CYTB, ND1, ND2, ND3, ND4 and ND5 as well as being commonly found in CO1 and ND1 (Fig. 5, S6 Table).

3.4. Transfer RNAs The typical complement of 22 typical tRNAs found in the arthropod mt genomes were identified in the A. miser mt genome, ranging in size from 64 bp (tRNAArg) to 72 bp (tRNAVal) and 1478 bp in total. The overall

191

A + T content of tRNAs is 77.0%. Fourteen genes are encoded on the Jstrand and the remains are encoded on the N-strand. Most tRNAs could be folded into the typical clover-leaf structure (Fig. 6), whereas tRNASer(AGN) is an exception lacking a DHU arm as has been observed for this gene in other metazoan mt genomes (Wolstenholme, 1992). A total of 14 G–U pairs and 5 mismatched base U\\U pairs were found in A. miser mt tRNA secondary structures, no A–A pairs or A–C pairs were found. The G–U pairs are located in the AA arm (4 bp), DHU arm (7 bp), AC arm (3 bp), respectively. The mismatched base U\\U pairs are located in the TψC arm (2 bp) and AA arm (3 bp).

3.5. Ribosomal RNAs The rRNAs of A. miser are 1333 bp for lrRNA and 795 bp for srRNA in length. Their A + T content are 81.9% and 79.2%, respectively. Because rRNA genes lack functional annotation features, analogous to the start and stop codons of PCGs, it is impossible to determine the boundaries from DNA sequence alone (Boore, 2001, 2006a). Hence, they were assumed to extend to the boundaries of flanking genes. As in other dipteran species, the lrRNA gene is flanked by tRNA Leu (CUN) and tRNAVal, while the srRNA gene is flanked by tRNAVal and the control region. We inferred secondary structures of lrRNA and srRNA of A. miser using the published rRNA secondary structures of a sepsid fly Nemopoda mamaevi, as a model (X.K. Li et al., 2015). The lrRNA has 42 helices in five structural domains (I-II, IV-VI, domain III is absent as in other insects), which is similar to the typical structure of arthropods (Cannone et al., 2002). The secondary structure of srRNA in A. miser includes three domains and 25 helices, again similar to other Diptera (Yang et al., 2011) (Figs. 7, 8).

Fig. 5. Usage of start and stop codons in Orthorrhapha mt genomes. (A) Start codons usage of 13 PCGs in Orthorrhapha. (B) Stop codons usage of 13 PCGs in Orthorrhapha.

192

K. Wang et al. / Gene 586 (2016) 184–196

Fig. 6. Putative secondary structures of tRNAs found in the mt genome of A. miser. All tRNAs can be folded into the usual clover-leaf secondary structure. The tRNAs are labelled with the abbreviations of their corresponding amino acids. Inferred Watson–Crick bonds are illustrated by lines, whereas GU bonds are illustrated by dots.

3.6. The control region The length of the control region of A. miser is 993 bp, shorter than that of C. duplonotata (1376 bp) and T. punctata (1597 bp). It has an A + T content of 87.8%, which is lower than that of C. duplonotata (92.6%) and higher than that of T. punctata (81.2%). All three species have weak AT-skews. A. miser has a strongly negative GC-skew,

while C. duplonotata and T. punctata's have strongly positive GCskew. 3.7. Phylogeny Four datasets were used in the phylogenetic analyses. There are 11,163 residues in the PCG123 matrix (containing nucleotides of 13

K. Wang et al. / Gene 586 (2016) 184–196

193

Fig. 7. Predicted secondary structure of the lrRNA gene in A. miser. Inferred Watson–Crick bonds are illustrated by lines, GU bonds by dots.

PCGs), 14,453 residues in the PCG123RNA matrix (containing nucleotides of 13 PCGs, two rRNAs and 19 tRNAs), 7442 residues in the PCG12 matrix (nucleotides of 13 PCGs excluding the third codon sites) and 10,732 residues in the PCG12RNA matrix (nucleotides of 13 PCGs excluding the third codon sites, two rRNAs and 19 tRNAs). The phylogenetic trees conducted with both Bayesian and ML analyses have very similar topologies across four datasets (Fig. 9). The monophyly of Stratiomyomorpha was consistently supported (posterior probability = 1.00, ML bootstrap = 100), as was the monophyly of Tabanomorpha (posterior probability = 1, ML bootstrap = 99/99/98/98), Xylophagomorpha (posterior probability = 1/0.88/1/0.94, ML bootstrap = 66/3/0/2), Muscomorpha (posterior probability = 1/1/1/0.98, ML bootstrap = 49/27/0/7), Asiloidea (posterior probability = 1, ML bootstrap = 99/100/0/100), and Empidoidea (posterior probability = 1/1/0.63/1, ML bootstrap = 84/92/26/100). These results support Yeates et al.'s and Wiegmann et al.'s findings regarding basal branching events within Orthorrhapha (Yeates, 2002; Wiegmann et al., 2011; Yeates and Wiegmann, 2007). A novel relationship between these four infraorders was consistently supported by our analyses: Tabanomorpha + ((Xylophagomorpha + Stratiomyomorpha) + Muscomorpha). However, this result was neither supported by Yeates et al., who hold the relationship that (Stratiomyomorpha + (Xylophagomorpha + Tabanomorpha)) + Muscomorpha (Yeates and Wiegmann, 2007; Yeates, 2002), nor proposed by Sinclair et al. which is Xylophagomorpha + Tabanomorpha + (Stratiomyomorpha + Muscomorpha) (Sinclair et al., 1994). Tabanidae was sister group of Athericidae (posterior probability = 1, ML bootstrap = 82/97/92/93), and these two families formed a clade as the sister-group of Rhagionidae (posterior probability = 1, ML bootstrap = 99/99/98/98). These relationships were supported by

Yeates et al.'s, Wiegmann et al.'s and Woodley's research (Woodley, 1989; Yeates, 2002; Wiegmann et al., 2011). Stratiomyidae was sister group of Xylomyidae (posterior probability = 1, ML bootstrap = 100) which is consistent with the previous study of Woodley, Yeates et al. and Wiegmann et al. (Woodley, 1989; Yeates, 2002; Wiegmann et al., 2011). Dolichopodidae was sister group of Empididae (posterior probability = 1/0.99/0.63/0, ML bootstrap = 84/92/26/0) and these two families formed superfamily Empidoidea, which had a sister-group relationship with Cyclorrhapha (posterior probability = 1/0.99/0.65/0.81, ML bootstrap = 94/99/41/75). These results support Yeates et al.'s finding (Yeates, 2002; Yeates and Wiegmann, 2007). Our result of the relationship within Muscomorpha was Asiloidae + (Empidoidea + Cyclorrhapha) (posterior probability = 1/0.99/0.65/ 0.81, ML bootstrap = 94/99/41/75), which is not conform to Wiegmann et al. (Wiegmann et al., 2011), but supported by previous studies based on morphological characters (Yeates and Wiegmann, 2007; Lambkin et al., 2013). The position of Nemestrinidae has been a debated question. Yeates et al. and Woodley placed it into Muscomorpha (Yeates and Wiegmann, 2007; Woodley, 1989), while Griffiths and Nagatomi suggested that it belongs to Tabanomorpha (Griffiths, 1994; Nagatomi, 1992). However, we detected a weakly supported sister-group relationship between Nemestrinidae and Xylophagidae (posterior probability = 1/0.88/ 1/0.94, ML bootstrap = 66/3/0/2), and therefore suggested that Nemestrinidae belongs to Xylophagomorpha. This result was also supported by Wiegmann et al. (Wiegmann et al., 2011). Apparently, considerable additional data of lower Brachycera are required to reliably resolve relationships between the four infraorders of Brachycera (Xylophagomorpha, Stratiomyomorpha, Tabanomorpha, Muscomorpha) and the position of Nemestrinidae.

194

K. Wang et al. / Gene 586 (2016) 184–196

Fig. 8. Predicted secondary structure of the srRNA gene in A. miser. Roman numerals denote the conserved domain structure. Inferred Watson–Crick bonds are illustrated by lines, GU bonds by dots.

4. Conclusions The mt genome of A. miser is 15,858 bp, which is typical of Diptera, with 13 protein-coding genes, 22 tRNA genes, two rRNA genes and a 993 bp

control region. The rest of the orthorrhaphous mt genomes in our study have similar structure with A. miser. Additionally, we conducted a phylogenetic analysis of 20 mt genomes using Maximum-likelihood and Bayesian methods in order to reconstruct the evolutionary relationship

K. Wang et al. / Gene 586 (2016) 184–196

195

Fig. 9. Phylogenetic tree of Brachycera families based on mt genome data. Cladogram of relationships resulting from Bayesian analyses with datasets PCG12, PCG123 and PCG12RNA, ML analyses with PCG12RNA and PCG123.With Tipula abdominalis (Tipulidae) and Chironomus tepperi (Chironomidae) as outgroups. Squares at the nodes are Bayesian posterior probabilities for 1, 2, 3 and 4, ML bootstrap values for 5, 6, 7 and 8. Dataset of PCG12, 1 and 5, PCG12RNA, 2 and 6, PCG123, 3 and 7, PCG123RNA, 4 and 8. Black indicates posterior probabilities = 1.00 or ML bootstrap = 100, gray indicates posterior probabilities ≥0.90 or ML bootstrap ≥70, white indicates posterior probabilities b0.90 or ML bootstrap b70, ‘ns’ indicates not support, star indicates posterior probabilities = 1.00 or ML bootstrap = 100 in eight trees.

of Orthorrhapha. The results show that all infraorders of Brachycera are monophyletic, and a relationship of Tabanomorpha + ((Xylophagomorpha + Stratiomyomorpha) + Muscomorpha) has been proposed. Within Xylophagomorpha, Nemestrinoidae forms the sister group of Xylophagidae. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2016.04.013. Competing interests The authors have declared that no competing interest exists. Acknowledgments DY was funded by the National Natural Science Foundation of China (31272354, 3141101151). NW was funded by the National Natural Science Foundation of China (41301049, 4151101023). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References Beckenbach, A.T., 2012. Mitochondrial genome sequences of Nematocera (Lower Diptera): evidence of rearrangement following a complete genome duplication in a winter crane fly. Genome Biol. Evol. 4 (2), 89–101. Beckenbach, A.T., Joy, J.B., 2009. Evolution of the mitochondrial genomes of gall midges (Diptera: Cecidomyiidae): rearrangement and severe truncation of tRNA genes. Genome Biol. Evol. 1, 278–287. Behura, S.K., Lobo, N.F., Haas, B., deBruyn, D., DD, Lovin, et al., 2011. Complete sequences of mitochondria genomes of Aedes aegypti and Culex quinquefasciatus and comparative analysis of mitochondrial DNA fragments inserted in the nuclear genomes. Insect Biochem. Mol. Biol. 41, 770–777.

Bernt, M., Bleidorn, C., Braband, A., Dambach, J., Donath, A., et al., 2013. A comprehensive analysis of bilaterian mitochondrial genomes and phylogeny. Mol. Phylogenet. Evol. 69, 252–364. Boore, J.L., 2001. Complete mitochondrial genome sequence of the polychaete annelid. Mol. Biol. Evol. 18, 1413–1416. Boore, J.L., 2006a. Requirements and standards for organelle genome databases. OMICS 10, 119–126. Boore, J.L., 2006b. The use of genome-level characters for phylogenetic reconstruction. Trends Ecol. Evol. 21, 439–446. Boore, J.L., Macey, J.R., Medina, M., 2005. Sequencing and comparing whole mitochondrial genomes of animals. Methods Enzymol. 395, 311–348. Cameron, S.L., 2014a. How to sequence and annotate insect mitochondrial genomes for systematic and comparative genomics research. Syst. Entomol. 39 (3), 400–411. Cameron, S.L., 2014b. Insect mitochondrial genomics: implications for evolution and phylogeny. Annu. Rev. Entomol. 59, 95–117. Cameron, S.L., Whiting, M.F., 2008. The complete mitochondrial genome of the tobacco hornworm, Manduca sexta (Insecta: Lepidoptera: Sphingidae), and an examination of mitochondrial gene variability within butterflies and moths. Gene 408, 112–123. Cameron, S.L., Lambkin, C.L., Barker, S.C., Whiting, M.F., 2007. A mitochondrial genome phylogeny of Diptera: whole genome sequence data accurately resolve relationships over broad timescales with high precision. Syst. Entomol. 32 (1), 40–59. Cameron, S.L., Lo, N., Bourguignon, T., Svenson, G.J., Evans, T.A., 2012. A mitochondrial genome phylogeny of termites (Blattodea: Termitoidae): robust support for interfamilial relationships and molecular synapomorphies define major clades. Mol. Phylogenet. Evol. 65, 163–173. Cannone, J.J., Subramanian, S., Schnare, M.N., Collett, J.R., D'Souza, L.M., et al., 2002. The comparative RNA web (CRW) site: an online database of comparative sequence and structure information for ribosomal, intron, and other RNAs. BMC Bioinformatics 3, 2. Chatterjee, S., Taraphdar, T., Mohandas, T.P., 2005. Molecular analysis of divergence in tachinid uzi (Exorista sorbillans) populations in India. Genetica 125, 1–15. Dowton, M., Castro, L., Austin, A.D., 2002. Mitochondrial gene rearrangements as phylogenetic characters in the invertebrates: the examination of genome ‘morphology’. Invertebr. Syst. 16, 345–356. Fenn, J.D., Song, H., Cameron, S.L., Whiting, M.F., 2008. A preliminary mitochondrial genome phylogeny of Orthoptera (Insecta) and approaches to maximizing phylogenetic signal found within mitochondrial genome data. Mol. Phylogenet. Evol. 49 (1), 59–68.

196

K. Wang et al. / Gene 586 (2016) 184–196

Gibson, A., Gowri-Shankar, V., Higgs, P.G., Rattray, M.A., 2005. Comprehensive analysis of mammalian mitochondrial genome base composition and improved phylogenetic methods. Mol. Biol. Evol. 22, 251–264. Gissi, C., Iannelli, F., Pesole, G., 2008. Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of conge-neric species. Heredity 101, 301–320. Grant, J.R., Stothard, P., 2008. The CGView server: a comparative genomics tool for circular genomes. Nucleic Acids Res. 36, W181–W184. Griffiths, G.C.D., 1994. Relationships among the major subgroups of Brachycera (Diptera): a critical review. Can. Entomol. 126, 861–880. Grimaldi, D., 1999. The co-radiations of pollinating insects and angiosperms in the Cretaceous. Ann. Mo. Bot. Gard. 86, 373–406. Hall, T.A., 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98. Junqueira, A.C., Lessinger, A.C., Torres, T.T., Da Silva, F.R., Vettore, A.L., et al., 2004. The mitochondrial genome of the blowfly Chrysomya chloropyga (Diptera: Calliphoridae). Gene 339, 7–15. Krzywinski, J., Li, C., Morris, M., Conn, J.E., Lima, J.B., et al., 2011. Analysis of the evolutionary forces shaping mitochondrial genomes of a Neotropical malaria vector complex. Mol. Phylogenet. Evol. 58 (3), 469–477. Lambkin, C.L., Sinclair, B.J., Pape, T., Courtney, G.W., Skevington, J.H., et al., 2013. The phylogenetic relationships among infraorders and superfamilies of Diptera based on morphological evidence. Syst. Entomol. 38, 164–179. Lanfear, R., Calcott, B., Ho, S.Y.W., Guindon, S., 2012. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analysis. Mol. Biol. Evol. 29, 1695–1701. Lavrov, D.V., Boore, J.L., Brown, W.M., 2000. The complete mitochondrial DNA sequence of the horseshoe crab Limulus polyphemus. Mol. Biol. Evol. 17, 813–824. Lee, E.S., Shin, K.S., Kim, M.S., Park, H., Cho, S., et al., 2006. The mitochondrial genome of the smaller tea tortrix Adoxophyes honmai (Lepidoptera: Tortricidae). Gene 373, 52–57. Li, H., Gao, J.Y., Cai, W.Z., 2015a. Complete mitochondrial genome of the assassin bug Oncocephalus breviscutum (Hemiptera: Reduviidae). Mitochondrial DNA http://dx. doi.org/10.3109/19401736.2013.840602. Li, X.K., Ding, S.M., Cameron, S.L., Kang, Z.H., Wang, Y.Y., et al., 2015b. The mitochondrial genome of sepsid fly Nemopoda mamaevi Ozerov with re-annotation and comparative protein-coding genes analysis of cyclorrhaphan flies. PLoS One 10 (3), e0123594. Lin, C.P., Danforth, B.N., 2004. How do insect nuclear and mitochondrial gene substitution patterns differ? Insights from Bayesian analyses of combined datasets. Mol. Phylogenet. Evol. 30, 686–702. Liu, Y., Cox, C.J., Wang, W., Goffinet, B., 2014. Mitochondrial phylogenomics of early land plants: mitigating the effects of saturation, compositional heterogeneity, and codonusage bias. Syst. Biol. 63, 862–878. Masta, S.E., Boore, J.L., 2008. Parallel evolution of truncated transfer RNA genes in arachnid mitochondrial genomes. Mol. Biol. Evol. 25, 949–959. Matsumoto, Y., Yanase, T., Tsuda, T., Noda, H., 2009. Species-specific mitochondrial gene rearrangements in biting midges and vector species identification. Med. Vet. Entomol. 23, 47–55. Moreno, M., Marinotti, O., Krzywinski, J., Tadei, W.P., James, A.A., et al., 2010. Complete mtDNA genomes of Anopheles darlingi and an approach to anopheline divergence time. Malar. J. 9, 127. Mousson, L., Dauga, C., Garrigues, T., Schaffner, F., Vazeille, M., et al., 2005. Phylogeography of Aedes (Stegomyia) aegypti (L.) and Aedes (Stegomyia) albopictus (Skuse) (Diptera: Culicidae) based on mitochondrial DNA variations. Genet. Res. 86, 1–11. Nagatomi, A., 1992. Notes on the phylogeny of various taxa of the orthorrhaphous Brachycera (Insecta: Diptera). Zool. Sci. 9, 843–857. Nardi, F., Carapelli, A., Fanciulli, P.P., Dallai, R., Frati, F., et al., 2001. The complete mitochondrial DNA sequence of the basal Hexapod Tetrodontophora bielanensis: evidence for heteroplasmy and tRNA translocations. Mol. Biol. Evol. 18 (7), 1293–1304. Nardi, F., Carapelli, A., Boore, J.L., Roderick, G.K., Dallai, R., et al., 2010. Domestication of olive fly through a multi-regional host shift to cultivated olives: comparative dating using complete mitochondrial genomes. Mol. Phylogenet. Evol. 57, 678–686. Nelson, L.A., Cameron, S.L., Yeates, D.K., 2011. The complete mitochondrial genome of the gall-forming fly, Fergusonina taylori Nelson and Yeates (Diptera: Fergusoninidae). Mitochondrial DNA 22 (3), 197–199.

Nelson, L.A., Cameron, S.L., Yeates, D.K., 2012. The complete mitochondrial genome of the flesh fly, Sarcophaga impatiens Walker (Diptera: Sarcophagidae). Mitochondrial DNA 23 (2), 42–43. Ojala, D., Montoya, J., Attardi, G., 1981. tRNA punctuation model of RNA processing in human mitochondria. Nature 290, 470–474. Oliveira, M.T., Barau, J.G., Martins-Junqueira, A.C., Feija, P.C., da Rosa, A.C., et al., 2008. Structure and evolution of the mitochondrial genomes of Haematobia irritans and Stomoxis calcitrans: the Muscidae (Diptera: Calyptratae) perspective. Mol. Phylogenet. Evol. 48, 850–857. Perna, N.T., Kocher, T.D., 1995. Patterns of nucleotide composition at fourfold degenerate sites of animal mitochondrial genomes. J. Mol. Evol. 41, 353–358. Ren, D., 1998. Flower-associated Brachycera flies as fossil evidence for Jurassic angiosperm origins. Science 280, 85–88. Ronquist, F., Huelsenbeck, J.P., 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. Rota-Stabelli, O., Kayal, E., Gleeson, D., Daub, J., Boore, J.L., et al., 2010. Ecdysozoan mitogenomics: evidence for a common origin of the legged invertebrates, the Panarthropoda. Genome Biol. Evol. 2, 425–440. ScheVer, S.J., Lewis, M.L., 2006. Mitochondrial phylogeography of the vegetable pest Liriomyza trifolii (Diptera: Agromyzidae): diverged clades and invasive populations. Ann. Entomol. Soc. Am. 99, 991–998. Shao, R., Barker, S.C., 2003. 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. Mol. Biol. Evol. 20 (3), 362–370. Sheffield, N.C., Song, H., Cameron, S.L., Whiting, M.F., 2008. A comparative analysis of mitochondrial genomes in Coleoptera (Arthropoda: Insecta) and genome descriptions of six new beetles. Mol. Biol. Evol. 25 (11), 2499–2509. Simon, C., Buckley, T.R., Frati, F., Stewart, J.B., Beckenbach, A.T., 2006. In-corporating molecular evolution into phylogenetic analysis, and a new compilation of conserved polymerase chain reaction primers for animal mitochondrial DNA. Annu. Rev. Ecol. Evol. Syst. 37, 545–579. Sinclair, B.J., Cumming, J.M., Wood, D.M., 1994. Homology and phylogenetic implications of male genitalia in Diptera-lower Brachycera. Entomol. Scand. 24, 407–432. Stamatakis, A., 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analysis with thousands of taxa and mixed models. Bioinformatics 22, 2688–2690. Stamatakis, A., Hoover, P., Rougemont, J., 2008. A rapid bootstrap algorithm for the RAxML web servers. Syst. Biol. 57, 758–771. Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., et al., 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. Torres, T.T., Dolezal, M., Schlotterer, C., Ottenwalder, B., 2009. Expression profiling of Drosophila mitochondrial genes via deep mRNA sequencing. Nucleic Acids Res. 37 (22), 7509–7518. Vaidya, G., Lohman, D.J., Meier, R., 2010. SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information. Cladistics 27, 171–180. Wei, S.J., Shi, M., Chen, X.X., Sharkey, M.J., van Achterberg, C., et al., 2010. New views on strand asymmetry in insect mitochondrial genomes. PLoS One 5, e12708. Wiegmann, B.M., Trautwein, M.D., Winkler, I.S., Barr, N.B., Kim, J.W., et al., 2011. Episodic radiations in the fly tree of life. PNAS 108 (14), 5690–5695. Wolstenholme, D.R., 1992. Animal mitochondrial DNA: structure and evolution. Int. Rev. Cytol. 141, 173–216. Woodley, N.E., 1989. Phylogeny and classification of the ‘Orthorraphous’ Brachycera. Manual of Nearctic Diptera Vol. 3 pp. 1371–1395. Yang, F., Du, Y.Z., Wang, L.P., Cao, J.M., Yu, W.W., 2011. The complete mitochondrial genome of the leafminer Liriomyza sativae (Diptera: Agromyzidae): great difference in the A + T-rich region compared to Liriomyza trifolii. Gene 485, 7–15. Ye, J., Fang, R., Yi, J.P., Zhou, G.L., Zheng, J.Z., 2010. The complete sequence determination and analysis of four species of Bactrocera mitochondrial genome. Plant Quar. 24 (3), 11–14. Yeates, D.K., 2002. Relationships of extant lower Brachycera (Diptera): a quantitative synthesis of morphological characters. Zool. Scr. 31, 105–121. Yeates, D.K., Wiegmann, B.M., 2007. Phylogeny and systematic of Diptera: two decades of progress and prospects. Zootaxa 1668, 565–590.