Research Article
Print
Research Article
Complete mitochondrial genomes of two flat-backed millipedes by next-generation sequencing (Diplopoda, Polydesmida)
expand article infoYan Dong, Lixin Zhu, Yu Bai, Yongyue Ou, Changbao Wang
‡ Molecular Biology Laboratory, College of Biology and food technology, 1528 Fengle road, Chuzhou University, Chuzhou, Anhui, China, 239000, Chuzhou, China
Open Access

Abstract

A lack of mitochondrial genome data from myriapods is hampering progress across genetic, systematic, phylogenetic and evolutionary studies. Here, the complete mitochondrial genomes of two millipedes, Asiomorpha coarctata Saussure, 1860 (Diplopoda: Polydesmida: Paradoxosomatidae) and Xystodesmus sp. (Diplopoda: Polydesmida: Xystodesmidae) were assembled with high coverage using Illumina sequencing data. The mitochondrial genomes of the two newly sequenced species are circular molecules of 15,644 bp and 15,791 bp, within which the typical mitochondrial genome complement of 13 protein-coding genes, 22 tRNAs and two ribosomal RNA genes could be identified. The mitochondrial genome of A. coarctata is the first complete sequence in the family Paradoxosomatidae (Diplopoda: Polydesmida) and the gene order of the two flat-backed millipedes is novel among known myriapod mitochondrial genomes. Unique translocations have occurred, including inversion of one half of the two genomes with respect to other millipede genomes. Inversion of the entire side of a genome (trnF-nad5-trnH-nad4-nad4L, trnP, nad1-trnL2-trnL1-rrnL-trnV-rrnS, trnQ, trnC and trnY) could constitute a common event in the order Polydesmida. Last, our phylogenetic analyses recovered the monophyletic Progoneata, subphylum Myriapoda and four internal classes.

Keywords

gene order, mitochondrial genome, Myriapoda, next-generation sequence, Polydesmida

Introduction

The Myriapoda comprise more than 18,000 species worldwide and are a diverse and ecologically important group of terrestrial arthropods (Miyazawa et al. 2014). Four groups have been united as Myriapoda: the Chilopoda (centipedes) and Diplopoda (millipedes), containing the vast majority of species, and the poorly investigated Pauropoda and Symphyla. An ancient group, the myriapods inhabited Pangea, and as a result they occur on all continents today except for Antarctica (Edgecombe 2004, 2010). Diplopoda (millipedes) is the third most diverse class of Arthropoda, with more than 11,000 species described and an estimated diversity of 80,000 species (Sierwald and Bond 2007; Shelley 2007; Zhang 2011; Brewer et al. 2012; Enghoff et al. 2015). Millipedes are an important component of modern terrestrial ecosystems and play a major role in the breakdown of organic matter (Hopkin and Read 1992). However, there are few studies documenting aspects of the group’s phylogeny, evolution, behavior, physiology, and ecology (Hoffman et al. 2002).

Extant myriapods have not been the subject of extensive molecular phylogenetic study and studies that have been done tend to focus on relationships at the phylum level. Comparative morphological, molecular and higher-level systematic evidence has largely confirmed the monophyly of myriapods (Boudreaux 1979; Ax 1999; Bitsch and Bitsch 2004; Gai et al. 2006; Bäcker et al. 2008; Regier et al. 2008, 2010; Rehm et al. 2014; Fernández et al. 2016), even though this was once a controversial topic (Pocock 1893; Snodgrass 1952; Dohle 1980). The phylogeny of the millipedes is still an open topic regarding their position within the Myriapoda and earliest splits inside the diplopod lineage (Edgecombe 2011; Fernández et al. 2016).

Mitochondrial genomes are used extensively to study phylogeography and phylogenetic relationships (Boore et al. 1998; Gissi et al. 2008; Cameron 2014). In addition to the sequences of mitochondrial genes, the secondary structures of RNAs as well as the mitochondrial gene order have been explored in a phylogenetic context (Boore 1999; Carapelli et al. 2004; Chen et al. 2011; Song et al. 2016). As in most other arthropods, myriapod mitochondrial (mt) genomes are a single circular DNA molecule encoding 13 proteins, 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and one A+T-rich region involved in the control of mtDNA replication and transcription.

Gene order in the centipede Cermatobius longicornis and the millipede Prionobelum sp. is identical to that of Limulus polyphemus (Arthropoda: Xiphosura) (Lavrov et al. 2000a; Dong et al. 2012b; Brewer et al. 2013; Gai et al. 2013). However, the arrangement of genes in mt genomes is remarkably variable in Strigamia maritima (Chilopoda: Geophilomorpha) and Symphylella sp. (Symphyla: Scolopendrellidae) (Gai et al. 2008; Robertson et al. 2015). All millipedes in which the mt genome has been sequenced, except Sphaerotheriida, have nad6 + cob placements that differ from that of Limulus Polyphemus and the nad6 + cob pattern was supposed to be sound molecular evidence supporting the Helminthomorpha clade. Dong et al. (2012a) compared nine known myriapod mt genomes and posited that a translocation of trnT out of the 5' end of nad4L is a common event in derived progoneate lineages. Although taxon sampling is limited, gene synteny has supplied evolutionary evidence relating to Myriapoda phylogenetic and evolutionary history. Full mitochondrial genomes of sixteen myriapod species have hitherto been sequenced, however, this number is still far from sufficient considering the high species richness of this group (Sierwald and Bond 2007; Budd and Telford 2009; Dong et al. 2012a; Robertson et al. 2015). This lack of mitochondrial genome data is hampering phylogenetic and evolutionary studies within the subphylum Myriapoda.

Compositional heterogeneity and accelerated substitution rates have proven to be major sources of systematic bias in mtDNA based phylogeny (Rota-Stabelli et al. 2011). Avoiding inadequate outgroups, selecting conserved amino acid alignment regions and bolstering taxon sampling are keys to phylogenetic reconstruction using mt genomes (Rota-Stabelli et al. 2011; Chen et al. 2014; Robertson et al. 2015). However, complete representation of the four myriapod classes in many studies is not included (Rota-Stabelli et al. 2011; Brewer et al. 2013) and the mt genomes of class Pauropoda, the presumed sister lineage of millipedes, are available in Genbank but not included in previous studies (Brewer et al. 2013). Relationships among the Myriapoda remain unresolved in Robertson et al. (2015), including among the four classes.

Prior to our study, the mt genomes of one flat-backed millipede representing the family Xystodesmidae: Appalachioria falcifera (Keeton, 1959) was sequenced. In this species, the entire side of the mt genome is inverted and all genes are on a single strand. Whole genome shotgun reads sequenced with the Illumina sequencing platform were used to obtain two complete mt genomes from the millipedes Asiomorpha coarctata and Xystodesmus sp. These species are representatives of the families Paradoxosomatidae and Xystodesmidae and further our understanding of how gene rearrangement occurred in the Polydesmida. Phylogenetic analysis including myriapods, three other arthropod classes and outgroups (species of Onychophora and Priapulida) were also performed to explore the internal relationships within the Myriapoda using sequence alignments from mitochondrial genes.

Methods

Taxon sampling, DNA extraction and PCR

Specimens of Asiomorpha coarctata and Xystodesmus sp. were collected from Langya Mountain, Chuzhou, Anhui, China (32°16'N, 118°16'E) and stored at the Molecular Biology Laboratory, Chuzhou University, Chuzhou, Anhui, China (MBLCZU). Species identification was performed by Dong Y (first author) and Qian CY (Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China). Voucher specimens (MBLCZU000145, MBLCZU000146) were deposited at the MBLCZU. Total genomic DNA was extracted from one individual representing each species using the DNeasy tissue Kit (Qiagen China, Shanghai).

Standard PCR reactions to amplify three different fragments of mtDNA (cox1, cob and nad5) were undertaken for each sample. Primers are listed in Suppl. material 1: Table S1. Amplified PCR products were gel-purified and then analyzed by primer walking on an ABI-PRISM3730 Automated DNA Sequencer.

Genome sequencing and analyses

For Illumina sequencing, double index sequencing libraries with average insert sizes of around 300 bp were prepared. The libraries were sequenced as 250 bp paired-end runs on an Illumina Hi-Seq 2000 (about 2 Gb raw data each species). The resulting bait sequences (cox1, cob and nad5) were subsequently employed as references in the manner detailed below. De novo assemblies were conducted with Geneious v8.1 using the Map to Reference program with the following settings applied: medium-low sensitivity/fast; iterate up to five times; with a maximum of 2% mismatches, a maximum gap size of 3 bp and requiring a minimum overlap of 100 bp; do not trim. After manual inspection, the longest contigs resulting from the respective assemblies were then aligned and ensured for correct translation frames with MEGA v5.0, together with reference protein-coding gene sequences from seven millipedes (Narceus annularus, Thyropygus sp., Antrokoreana gracilipes, Appalachioria falcifera, Abacion magnum, Brachycybe lecontii, Prionobelum sp.). The contig ends of A. coarctata and Xystodesmus sp. overlapped.

Mitochondrial genome annotation and analyses

The assembled consensus sequence of each millipede mtDNA was further annotated and analyzed. Preliminary annotation using MITOS webserver provided overall information on mt genomes (Bernt et al. 2013). Protein-coding genes were annotated by identification of their open reading frames, and alignments of homologous genes of other reported myriapod mt genomes. Blast searches in the National Center for Biotechnology Information also helped to identify and annotate the PCGs. Transfer RNA genes were identified by comparing the results predicted by tRNAscan-SE Search Server v.1.21 and ARWEN based on cloverleaf secondary structure information (Lowe and Eddy 1997; Laslett and Canback 2008). Based on known gene order information, the boundaries of the 16S rRNA (rrnS) gene were assumed to be delimited by the ends of the trnV and trnL2 pair. The 12S rRNA (rrnL) gene was assumed to start from the end of trnV, and its end was roughly identified by alignment with other published millipede sequences. Nucleotide frequencies and codon usage were determined by MEGA v5.05 (Tamura et al. 2007).

Data sets and sequence alignment

Phylogenetic trees were reconstructed focused on the Arthropoda defining Onychophora and Priapulida as outgroups in the analyses, including 31 taxa (Suppl. material 1: Table S2). Amino acid sequences of 13 mitochondrial protein-coding genes were aligned separately using Clustal X 1.81 based on default settings for each gene (Thompson et al. 1997) then used as guides to align nucleotide sequences. GenBank accession numbers for taxa used in this study are given in Suppl. material 1: Table S1. Each alignment was analyzed with Gblocks under default settings to select conserved amino acid aligned regions (Castresana 2000).

Phylogenetic analyses

The optimal partition strategy and models were selected by PartitionFinder v1.1.1. We created an input configuration file that contained 13 pre-define partitions by gene. We used the ‘greedy’ algorithm with branch lengths estimated as ‘unlinked’ and Akaike Information Criterion (AIC) to search for the best-fit scheme (Suppl. material 1: Table S3). Model selection of the amino acid data set was performed with ProtTest version 2.4 (Abascal et al. 2007), and under AIC, MtRev+I+G was the best-fit model.

Maximum likelihood phylogenetic analysis searches were carried out with RAxML through a web portal (http://phylobench.vital-it.ch/raxml-bb/index.php). Bootstrap values, indicating the robustness of the internal nodes of the gene trees, were set at 100 replicates. Bayesian analyses of nucleotide and amino acid data sets were performed with MrBayes v3.1.2 (Ronquist and Huelsenbeck 2003), using the GTR+I+G model and MtRev+I+G model, respectively. Four Markov chains were run for 2×106 generations and sampled every 100 generations to yield a posterior probability distribution of 2×104 trees. The first 2000 trees were discarded as burn-in. Three replicates of these Bayesian runs were conducted, retrieving the same topology.

Results

Mitochondrial genome organization

The complete mt genome of Asiomorpha coarctata (GenBank accession KU721885) is 15,644 bp long and was assembled with coverage between 54–1062 reads; the genome of Xystodesmus sp. (GenBank accession KU721886) is 15,791 bp in length and has a coverage within 86–1392 reads. The sizes of these two mt genomes are within the range reported for those of other myriapods ranging from 14,487–16,833 bp (GenBank accessions NC_016676 and NC_021403).

The two genomes encode 13 protein-coding regions, 22 tRNA and two rRNA genes, consistent with metazoan mitochondrial DNA structure, but contain a number of unique features (Figure 1). Most protein-coding genes start with the codon ATG, with the exception of nad1, nad2, nad3, and nad4, which begins with TTG, TTG, GTG, and GTG, respectively. Seven protein-coding genes use the typical termination codons TAG (cox1, cox2, nad4L, nad4 and nad5) and TAA (atp8 and nad1), while others in the mt genome of A. coarctata use incomplete stop codons (Table 1). Several genes show complete stop codons, either TAG (as in cox1, cox2, nad3, nad4L and nad5) or TAG (as in atp8 and cytb) in Xystodesmus sp. (Table 2). Incomplete stop codons are frequently found in other myriapod mitochondrial protein-encoding genes (Dong et al. 2012a) and may be completed by polyadenylation after cleavage of the polycistronic transcript (Ojala et al. 1981). Both novel genomes have an overlapping gene region only between atp8/atp6.

Figure 1.

Mitochondrial genomes of the two millipedes sequenced in this study. A Asiomorpha coarctata B Xystodesmus sp. Circular maps were drawn with Geneious v9.1.2. Arrows indicate the orientation of gene transcription. Abbreviations of gene names are: atp6 and atp8 for ATP synthase subunits 6 and 8; cox13 for cytochrome oxidase subunits 1–3; cob for cytochrome b, nad16 and nad4L for NADH dehydrogenase subunits 1–6 and 4L; and lrRNA and srRNA for large and small rRNA subunits. tRNA genes are indicated with their one-letter corresponding amino acids. CR for control region. The GC content was plotted using a green sliding window and the AT content was blue.

Organization of the mitochondrial genome of Asiomorpha coarctata.

Feature From To Length (nt) Codons Spacer/Overlap(-)
Start Stop
cox1 1 1533 1533 ATG TAG 11
cox2 1545 2222 678 ATG TAG 8
trnK 2231 2297 67 CTT 1
trnD 2299 2363 65 GTC 0
atp8 2364 2522 159 ATG TAA -7
atp6 2516 3182 667 ATG TA 0
cox3 3183 3967 785 ATG TA 0
trnG 3968 4031 64 TCC 0
nad3 4032 4381 334 GTG TA 0
trnA 4382 4443 62 TGC 0
trnR 4444 4509 66 TCG 2
trnN 4512 4579 68 GTT 0
trnS1 4580 4638 69 GCT 0
trnE 4639 4703 65 TTC 0
nad6 4704 5176 473 ATG AGT 0
cytb 5177 6297 1121 ATG TA 0
trnS2 6298 6358 61 TGA 0
CR1 6359 7321 963 0
rrnS 7322 8097 776 0
trnV 8098 8162 65 TAC 0
rrnL 8163 9402 1240 0
trnL1 9403 9465 63 TAG 0
trnL2 9469 9533 65 TAA 3
nad1 9534 10460 TTG TAA 0
trnP 10461 10524 64 CCA 0
nad4L 10525 10806 62 ATG TAG 0
nad4 10800 12143 1344 GTG TAG 0
trnT 12144 12206 63 TGT 0
CR2 12207 12393 187 0
trnH 12394 12458 65 GTG 1
nad5 12460 14160 1701 ATG TAG 3
trnF 14164 14231 68 GAA 0
trnY 14232 14299 68 GTA 1
trnQ 14301 14371 71 TTG 0
trnC 14372 14436 65 GCA 0
trnI 14437 14503 67 GAT 0
trnM 14504 14570 67 CAT 0
nad2 14571 15575 1005 TTG GGA 0
trnW 15576 15644 69 TCA 0
Total 15644 -7/30

Organization of the mitochondrial genome of Xystodesmus sp.

Feature From To Length (nt) Codons Spacer/Overlap(-)
Start Stop
cox1 1 1533 1533 ATG TAG 6
cox2 1540 2217 678 ATG TAG 3
trnK 2221 2287 67 CTT 0
trnD 2288 2351 64 GTC 0
atp8 2352 2513 162 ATG TAA -7
atp6 2507 3175 669 ATG TA 0
cox3 3179 3963 785 ATG TA 0
trnG 3964 4029 66 TCC 0
nad3 4030 4380 351 GTG TAG 1
trnA 4382 4445 64 TGC 4
trnR 4450 4514 65 TCG 3
trnN 4518 4579 68 GTT 0
trnS1 4580 4640 61 GCT 2
trnE 4643 4710 68 TTC 2
nad6 4713 5181 469 ATG T 0
cytb 5182 6303 1122 ATG TAA 0
trnS2 6304 6360 57 TGA 0
CR1 6361 7392 963 0
rrnS 7393 8172 780 0
trnV 8173 8238 66 TAC 0
rrnL 8239 9465 1227 0
trnL1 9466 9530 65 TAG 51
trnL2 9582 9650 69 TAA 3
nad1 9651 10575 925 TTG T 0
trnP 10576 10640 65 TGG 2
nad4L 10643 10924 282 ATG TAG 0
nad4 10918 12264 1347 GTG TAA 1
trnT 12266 12327 62 TGT 0
CR2 12328 12535 208 0
trnH 12536 12599 64 GTG 1
nad5 12602 14299 1698 ATG TAG 3
trnF 14303 14369 67 GAA 0
trnY 14370 14434 65 GAT 15
trnQ 14450 14515 66 TTG 2
trnC 14518 14582 65 GCA 3
trnI 14586 14650 65 GAT 0
trnM 14651 14716 66 CAT 0
nad2 14717 15722 1006 TTG T 0
trnW 15723 15790 68 TCA 1
Total 15791 -7/103

Non-coding regions

Some millipede mt genomes that have been sequenced (e.g. Antrokoreana gracilipes) include two major non-coding regions and others contain a single non-coding region, such as Narceus annularus, Prionobelum sp. and Appalachioria falcifera (Figure 5). Of the genomes sequenced here, Asiomorpha coarctata and Xystodesmus sp. include two major non-coding regions (Table 1 and 2).

The largest non-coding region (CR1, 963 bp) is located between trnS2 and rrnS in A. coarctata (Table 1). The non-coding region in A. coarctata contains tandemly repeated regions (11.4 × 38 bp), and the repeated unit is ‘GTAATAATATAGATAGAGTAATATAACCTTATATAGGA’ (Figure 2).

Figure 2.

Sequences of the non-coding region in Asiomorpha coarctata, primary structures of tandemly repeated regions (11.4 × 38 bp).

Transfer RNA

There are 22 potential tRNA genes in Asiomorpha coarctata and Xystodesmus sp., respectively (Figure 3 and 4), as there are in most other published metazoan mtDNAs (Gissi et al. 2008). All of these tRNA genes are α-strand-encoded bearing more protein-coding sequence (Figure 1), and the newly sequenced mt genomes show dihydrouridine arm (DHU arm, D-arm) loss in trnS1 and trnS2. According to our analysis based on the ARWEN program, trnS1 lacks the D-arm in all other millipede species.

Figure 3.

Putative secondary structures of the 22 tRNA genes of Asiomorpha coarctata. Watson-Crick base-pairing is indicated by solid lines, and G–T pairs are indicated with plus signs.

Figure 4.

Putative secondary structures of the 22 tRNA genes of Xystodesmus sp. Watson-Crick base-pairing is indicated by solid lines, and G–T pairs are indicated with plus signs.

Gene order

Gene order arrangements were compared with mt genome organization in other myriapods (Figure 5), including the ancestral gene order of the mitogenome for the Myriapoda shared by Prionobelum sp. and Limulus polyphemus (Dong et al. 2012b). The overall arrangement of the genes around the A. coarctata and Xystodesmus sp. mt genomes is unique compared to other myriapod species. All coding regions are on a single strand in A. coarctata and Xystodesmus sp. which has been reported in the mt genome of the millipede Appalachioria falcifera with an entire side of the genome inverted (Brewer et al. 2013). These three flat-backed millipedes are identical but for a single tRNA translocation in A. coarctata and Xystodesmus sp. The two newly sequenced mt genomes have undergone gene and tRNA translocations compared to other myriapod sequences.

Figure 5.

Comparison of gene arrangements in mtDNA of the arthropod ground pattern. Gene segments are not drawn to scale. Genes shaded gray have different relative positions compared to the ground pattern. Underlining indicates the gene is encoded on the opposite strand, and arrows indicate translocation of trnT. CR: putative control region. Gene arrangements of two diplopods, N. annularus and Thyropygus sp. are similar and represented as one.

Phylogenetic inference

Bayesian inference and maximum likelihood phylogenetic analysis were performed using conserved blocks of amino acid and nucleotide data sets (Figure 6). The topological pattern obtained using Bayesian inference and maximum likelihood analyses based on both the amino acid and nucleotide data sets were similar.

Figure 6.

Phylogenetic tree of the Arthropoda, including Myriapoda, Hexapoda, Crustacea and Chelicerata and outgroups reconstructed based on protein-coding genes from mtDNA genomes. Each group of four numbers indicates node confidence values (from top left): Bayesian posterior probabilities in percent (BPP) in amino acid and nucleotide datasets; maximum likelihood bootstrapping values (MLBP) in amino acid and nucleotide datasets.

Within the Diplopoda, the Polydesmida clade, including Appalachioria falcifera, Asiomorpha coarctata, and Xystodesmus sp. is well supported. A clade consisting of Polydesmida plus Colobognatha (Brachycybe lecontii) is also recovered. All millipedes in this study except Prionobelum sp. recovered as a monophyletic group, Helminthomorpha. The grouping with Helminthomorpha and Pentazonia (the basal lineage Prionobelum sp.) is supported.

Four myriapod clades are resolved, with Chilopoda (BPP = 1.00 and 1.00; MLBP = 99 and 98) as the basal group. A monophyletic Progoneata (“Symphyla + Pauropoda” + Diplopoda) is recovered as the sister group of the Chilopoda (BPP = 1.00 and 0.99; MLBP = 89 and 80). Our Bayesian inference analysis recognizes a monophyletic Myriapoda with strong support, while the maximum likelihood analysis is only weakly supported.

Discussion

Next-generation sequencing such as the 454 pyrosequencing, Solexa, and SOLiD provided by Roche, Illumina and Applied Biosystems, has the ability to generate a large number of sequences within a very short time compared to Sanger’s method (Chilana et al. 2012). These methods have been fused to rapid generate sequence data and successful de novo sequence assemblies for arthropods (Kirkness et al. 2010) and more recently to the assembly of full mt genomes (Knaus et al. 2011; Ma et al. 2012; Coates 2014; Aguado et al. 2016). Next-generation sequencing results in better assembly, leading to fewer gaps, larger contigs and greater accuracy of the final consensus sequence. It is essential for accurately identifying more complex rearrangements, for example in the order Polydesmida, which has a large inversion in the mt genome in this study.

We found that the genome size was larger in A. coarctata (15,644 bp) and Xystodesmus sp. (15,791 bp) than other millipede mt genomes (14747–15282 bp; Antrokoreana gracilipes, GenBank accession NC_010221; Appalachioria falcifera, GenBank accession NC_021933). Intergenic spacer length variation may have arisen through retention of partial duplication, or incomplete multiple deletions of redundant genes (McKnight and Shaffer 1997; Yamauchi et al. 2003) under the duplication-and-deletion mechanism. For this reason, we speculate that multiple intergenic spacers distributed in the two larger mt genomes may serve as a guide in deducing derived gene arrangement.

Most of the tRNAs appear to be truncated and lack one of the arms found in the canonical tRNAs of pauropods and symphylans (Gai et al. 2008; Dong et al. 2012a). One to two nucleotide mismatches could be found in the acceptor arms of the five tRNAs in A. coarctata, and in the four tRNAs in Xystodesmus sp. Nucleotide mismatch in the arms of tRNAs may also occur in other myriapod groups (Woo et al. 2007; Lavrov et al. 2000b; Gai et al. 2008, 2013; Brewer et al. 2012; Dong et al. 2012a, 2012b).

Compared with other millipedes, the mt gene order in A. coarctata and Xystodesmus sp. is very similar to that in Appalachioria falcifera. These three species belong to the order Polydesmida (Figure 5). There are no differences in the relative position of the protein-coding genes, but the trnT gene and non-coding regions of Appalachioria falcifera mt genome are translocated with respect to the newly sequenced mt genomes here. Although only nine millipede mt genomes are compared, an extreme variety in gene arrangements is known in millipedes, and inversion of an entire side of the genome (trnF-nad5-trnH-nad4-nad4L, trnP, nad1-trnL2-trnL1-rrnL-trnV-rrnS, trnQ, trnC and trnY) could be a synapomorphy in the Polydesmida lineage.

Gene arrangement in three flat-backed millipedes is similar to that in other Helminthomorpha sequenced previously in which nad6 + cob placements occurred. We agree with Brewer et al. (2013) that the inversion of the mt genome in flat-backed millipedes is a derived event associated with losing the second non-coding region. The mitochondrial gene arrangements of the order Polydesmida and the infraclass Helminthomorpha lineages are reshuffled regularly. To better understand the evolutionary implications of gene arrangements in the Myriapoda, mt genome research with broader taxon sampling will be required.

This phylogenomic study provides a strongly supported phylogenetic framework for the monophyletic origin of the Myriapoda and the monophyly of the Progoneata and extant myriapod subgroups. Among the Myriapoda, the union of Diplopoda with Pauropoda as a monophyletic group (= Dignatha) was once widely accepted (Dohle 1980, 1998; Enghoff et al. 1993; Blanke and Wesener 2014; Edgecombe 2015), but a number of molecular analyses of nuclear and mitochondrial sequences support the combination of Pauropoda and Symphyla (Gai et al. 2006; Regier et al. 2010; Dong et al. 2012a). Our results strongly support the traditional morphology-based Progoneata (Diplopoda + Pauropoda + Symphyla) defined by the presence of gonopores behind the second pair of legs. The notion of Progoneata recovered here is consistent with that favored by morphology and molecular analyses (Gai et al. 2008; Edgecombe 2010, 2011; Regier et al. 2010; Dong et al. 2012a).

The Polydesmida and Helminthomorpha are recovered as monophyletic. Gene order in Prionobelum sp. which is the basal lineage of the millipede mt genome is assumed to represent the millipede or myriapod ground pattern (Dong et al. 2012b), and inversion of the entire side of the genome occurred as a common event in the order Polydesmida lineage proposed in this study. In our gene order comparison of these millipedes, phylogenetic results are mainly concordant with gene arrangement analyses (Figure 5). Combining the implications of phylogenetic analyses and gene arrangement could yield valuable understanding of myriapod evolutionary history.

Polydesmida has been considered the sister group of Juliformia and a ‘ring-forming’ clade from the perspective of morphology (Enghoff et al. 1993). Polydesmida also unite the Nematophora (including Stemmiulida, Callipodida and Chordeumatida) as sister group for the shared presence of preanal spinnerets (Enghoff 1984; Sierwald et al. 2003; Shear 2008; Blanke and Wesener 2014). The Polydesmida allied with Platydesmida (basal representatives of the subclass Colobognatha) as a clade in our phylogenetic analyses. The combination Polydesmida + Colobognatha has been recovered in previous analyses (Sierwald and Bond 2007; Brewer et al. 2013). More taxa must be sequenced from the Colobognatha and better analysis methods used to test the position of the Polydesmida.

Many internal relationships of the Diplopoda remain unresolved and several groups are paraphyletic. The Bayesian inference tree seems to have better support at shallower nodes with amino acid data sets, whereas the maximum likelihood tree has better support at deeper nodes. The maximum likelihood tree has very low support at most nodes and confuses the relationships of taxa that are confidently placed in monophyletic groups by other studies (Meusemann et al. 2010; Regier et al. 2010; Brewer et al. 2013). Our finding that the Bayesian methods outperformed likelihood-based approaches is consistent with the results reported by Talavera and Castresana (2007).

Acknowledgements

The authors declare no competing interests. Special thanks to Pavel Stoev and three reviewers for constructive comments that improved this paper. This research was supported by the National Natural Science Foundation of China Project (31401971), Anhui Provincial Natural Science Foundation (1408085QC67), and Natural Science Foundation of the Higher Education Institutions of Anhui Province (KJ2013Z245).

References

  • Abascal F, Posada D, Zardoya R (2007) MtArt: a new model of amino acid replacement for Arthropoda. Molecular Biology and Evolution 24: 1–5. doi: 10.1093/molbev/msl136
  • Aguado MT, Grande C, Gerth M, Bleidorn C, Noreña C (2016) Characterization of the complete mitochondrial genomes from Polycladida (Platyhelminthes) using next-generation sequencing. Gene 575: 199–205. doi: 10.1016/j.gene.2015.08.054
  • Ax P (1999) Das System der Metazoa II. Ein Lehrbuch der phylogenetisch en Systematik. Gustav Fischer Verlag, Stuttgart, 384 pp.
  • Bäcker H, Fanenbruck M, Wägele JW (2008) A forgotten homology supporting the monophyly of Tracheata: the subcoxa of insects and myriapods. Zoologischer Anzeiger 247: 185–207. doi: 10.1016/j.jcz.2007.11.002
  • Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, Pütz J, Middendorf M, Stadler PF (2013) MITOS: improved de novo metazoan mitochondrial genome annotation. Molecular Phylogenetics and Evolution 69(2): 313–319. doi: 10.1016/j.ympev.2012.08.023
  • Bitsch C, Bitsch J (2004) Phylogenetic relationships of basal hexapods among mandibulate arthropods: a cladistic analysis based on comparative morphological characters. Zoologica Scripta 33: 511–550. doi: 10.1111/j.0300-3256.2004.00162.x
  • Blankea A, Wesener T (2014) Revival of forgotten characters and modern imaging techniques help to produce a robust phylogeny of the Diplopoda (Arthropoda, Myriapoda). Arthropod Structure & Development 43(1): 63–75. doi: 10.1016/j.asd.2013.10.003
  • Boore JL (1999) Animal mitochondrial genomes. Nucleic Acids Research 27: 1767–1780. doi: 10.1093/nar/27.8.1767
  • Boore JL, Lavrov D, Brown WM (1998) Gene translocation links insects and crustaceans. Nature 392(6677): 667–668. doi: 10.1038/33577
  • Boudreaux HB (1979) Significance of inter segmental tendon system in arthropod phylogeny and a monophyletic classification of Arthropoda. In: Gupta AP (Ed.) Arthropod Phylogeny. Van Nostrand Reinhold Company, New York, 551–586.
  • Brewer MS, Sierwald P, Bond JE (2012) Millipede Taxonomy after 250 Years: Classification and Taxonomic Practices in a Mega-Diverse yet Understudied Arthropod Group. PLoS ONE 7(5): e37240. doi: 10.1371/journal.pone.0037240
  • Brewer MS, Swafford L, Spruill CL, Bond JE (2013) Arthropod Phylogenetics in Light of Three Novel Millipede (Myriapoda: Diplopoda) Mitochondrial Genomes with Comments on the Appropriateness of Mitochondrial Genome Sequence Data for Inferring Deep Level Relationships. PLoS ONE 8(7): e68005. doi: 10.1371/journal.pone.0068005
  • Budd GE, Telford MJ (2009) The origin and evolution of arthropods. Nature 457: 812–817. doi: 10.1038/nature07890
  • Cameron SL (2014) Insect mitochondrial genomics: implications for evolution and phylogeny. Annual Review of Entomology 59: 95–117. doi: 10.1146/annurev-ento-011613-162007
  • Carapelli A, Soto-Adames FN, Simon C, Frati F, Nardi F, Dallai R (2004) Secondary structure, high variability and conserved motifs for domain III of 12S rRNA in the Arthropleona (Hexapoda; Collembola). Insect Molecular Biology 13: 659–670. doi: 10.1111/j.0962-1075.2004.00528.x
  • Castresana J (2000) Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Molecular Biology and Evolution 17(4): 540–552. doi: 10.1093/oxfordjournals.molbev.a026334
  • Chen WJ, Bu Y, Carapelli A, Dallai R, Li S, Yin WY, Luan YX (2011) The mitochondrial genome of Sinentomon erythranum (Arthropoda: Hexapoda: Protura): an example of highly dviergent evolution. BMC Evolutionary Biology 11: 246. doi: 10.1186/1471-2148-11-246
  • Chen WJ, Koch M, Mallatt JM, Luan YX (2014) Comparative analysis of mitochondrial genomes in Diplura (Hexapoda, Arthropoda): Taxon sampling is crucial for phylogenetic inferences. Genome Biology and Evolution 6(1): 105–120. doi: 10.1093/gbe/evt207
  • Chilana P, Sharma A, Rai A (2012) Insect genomic resources: status, availability and future. Current Science 102: 571–580.
  • Coates BS (2014) Assembly and annotation of full mitochondrial genomes for the corn rootworm species, Diabrotica virgifera virgifera and Diabrotica barberi (Insecta: Coleoptera: Chrysomelidae), using next generation sequence data. Gene 542: 190–197. doi: 10.1016/j.gene.2014.03.035
  • Dohle W (1980) Sind die Myriapoden eine monophyletische Gruppe? Eine Diskussionder Verwandtschaftsbeziehungen der Antennaten. Abhandl. Verhandl. Naturwiss. Ver Hamburg 23: 45–104.
  • Dohle W (1998) Myriapod–insect relationships as opposed to an insect–crustacean sister group relationship. In: Fortey RA, Thomas RH (Eds) Arthropod Relationships. Chapman and Hall, London, 305–315. doi: 10.1007/978-94-011-4904-4_23
  • Dong Y, Sun HY, Guo H, Pan D, Qian CY, Hao SJ (2012a) The complete mitochondrial genome of Pauropus longiramus (Myriapoda: Pauropoda): implications on early diversification of the myriapods revealed from comparative analysis. Gene 505(1): 57–65. doi: 10.1016/j.gene.2012.05.049
  • Dong Y, Xu JJ, Hao SJ, Sun HY (2012b) The complete mitochondrial genome of the giant pill millipede,Sphaerotheriidae sp. (Myriapoda: Diplopoda: phaerotheriida). Mitochondrial DNA 23(5): 333–335. doi: 10.3109/19401736.2012.683184
  • Edgecombe GD (2004) Morphological data, extant Myriapoda, and the myriapod stem-group. Contributions to Zoology Bijdragen Tot De Dierkunde 73: 207–252.
  • Edgecombe GD (2010) Arthropod phylogeny: an overview from the perspectives of morphology, molecular data and the fossil record. Arthropod Structure & Development 39(2–3): 74–87. doi: 10.1016/j.asd.2009.10.002
  • Edgecombe GD (2011) Phylogenetic relationships of Myriapoda. In: Minelli A (Ed.) The Myriapoda. Brill, Leiden, 1–20. doi: 10.1163/9789004188266_002
  • Edgecombe GD (2015) Diplopoda - Phylogenetic relationships. In: Minelli A (Ed.) Treatise on Zoology – Anatomy, Taxonomy, Biology. The Myriapoda, Volume 2, Chapter 2. Brill, Leiden, 353–362. doi: 10.1163/9789004188273_016
  • Enghoff H, Golovatch S, Short M, Stoev P, Wesener T (2015) Diplopoda – Taxonomic overview. In: Minelli A (Ed.) Treatise on Zoology – Anatomy, Taxonomy, Biology. The Myriapoda, Volume 2, Chapter 2. Brill, Leiden, 363–453. doi: 10.1163/9789004188273_017
  • Enghoff H, Dohle W, Blower JG (1993) Anamorphosis in millipedes (Diplopoda): the present state of knowledge with some developmental and phylogenetic considerations. Zoological Journal of the Linnean Society 109: 103–234. doi: 10.1111/j.1096-3642.1993.tb00305.x
  • Fernández R, Edgecombe GD, Giribet G (2016) Exploring phylogenetic relationships within Myriapoda and the effects of matrix composition and occupancy on phylogenomic reconstruction. Systematic Biology 65(5): 871–89. doi: 10.1093/sysbio/syw041
  • Gai Y, Ma H, Sun X, Ma J, Li C, Yang Q (2013) The complete mitochondrial genome of Cermatobius longicornis (Chilopoda: Lithobiomorpha: Henicopidae). Mitochondrial DNA 24(4): 331–332. doi: 10.3109/19401736.2012.760078
  • Gai Y, Song D, Sun H, Yang Q, Zhou K (2008) The complete mitochondrial genome of Symphylella sp. (Myriapoda: Symphyla): Extensive gene order rearrangement and evidence in favor of Progoneata. Molecular Phylogenetics and Evolution 49(2): 574–85. doi: 10.1016/j.asd.2009.10.002
  • Gai YH, Song DX, Sun HY, Zhou KY (2006) Myriapod monophyly and relationships among myriapod classes based on nearly complete 28S and 18S rDNA sequences. Zoological Science 23: 1101–1108. doi: 10.2108/zsj.23.1101
  • Gissi C, Iannelli F, Pesole G (2008) Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of congeneric species. Heredity 101: 301–320. doi: 10.1038/hdy.2008.62
  • Hoffman R, Golovatch S, Adis J, de Morais J (2002) Diplopoda. In: Adis J (Ed.) Amazonian Arachnida and Myriapoda. Pensoft, Sofia-Moscow, 505–533.
  • Kirkness EF, Haas BJ, Sun W, Braig HR, Perotti MA, Clark JM, Lee SH, Robertson HM, Kennedy RC, Elhaik E, Gerlach D, Kriventseva EV, Elsik CG, Graur D, Hill CA, Veenstra JA, Walenz B, Tubío JM, Ribeiro JM, Rozas J, Johnston JS, Reese JT, Popadic A, Tojo M, Raoult D, Reed DL, Tomoyasu Y, Kraus E, Mittapalli O, Margam VM, Li HM, Meyer JM, Johnson RM, Romero-Severson J, Vanzee JP, Alvarez-Ponce D, Vieira FG, Aguadé M, Guirao-Rico S, Anzola JM, Yoon KS, Strycharz JP, Unger MF, Christley S, Lobo NF, Seufferheld MJ, Wang N, Dasch GA, Struchiner CJ, Madey G, Hannick LI, Bidwell S, Joardar V, Caler E, Shao R, Barker SC, Cameron S, Bruggner RV, Regier A, Johnson J, Viswanathan L, Utterback TR, Sutton GG, Lawson D, Waterhouse RM, Venter JC, Strausberg RL, Berenbaum MR, Collins FH, Zdobnov EM, Pittendrigh BR (2010) Genome sequences of the human body louse and its primary endosymbiont provide insights into the permanent parasitic lifestyle. Proceedings of the National Academy of Sciences of USA 107(27): 12168–12173. doi: 10.1073/pnas.1003379107
  • Knaus BJ, Cronn R, Liston A, Pilgrim K, Schwartz MK (2011) Mitochondrial genome sequences illuminate maternal lineages of conservation concern in a rare carnivore. BMC Ecology 11: 10. doi: 10.1186/1472-6785-11-10
  • Laslett D, Canback B (2008) ARWEN: a program to detect tRNA genes in metazoan mitochondrial nucleotide sequences. Bioinformatics 24: 172–175. doi: 10.1093/bioinformatics/btm573
  • Lavrov D, Boore J, Brown W (2000a) The complete mitochondrial DNA sequence of the Horseshoe Crab Limulus polyphemus. Molecular Biology and Evolution 17(5): 813–824. doi: 10.1093/oxfordjournals.molbev.a026360
  • Lavrov D, Brown W, Boore J (2000b) A novel type of RNA editing occurs in the mitochondrial tRNAs of the centipede Lithobius forficatus. Proceedings of the National Academy of Sciences of USA 97: 13738–13742. doi: 10.1073/pnas.250402997
  • Lowe TM, Eddy SR (1997) tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Research 25: 955–964. doi: 10.1093/nar/25.5.0955
  • Ma PF, Guo ZH, Li DZ (2012) Rapid sequencing of the bamboo mitochondrial genome using Illumina technology and parallel episodic evolution of organelle genomes in grasses. PLoS ONE 7: e30297. doi: 10.1371/journal.pone.0030297
  • Mcknight ML, Shaffer HB (1997) Large, rapidly evolving intergenic spacers in the mitochondrial DNA of the Salamander family Ambystomatidae (Amphibia: Caudata). Molecular Biology and Evolution 14: 1167–1176. doi: 10.1093/oxfordjournals.molbev.a025726
  • Meusemann K, von Reumont BM, Simon S, Roeding F, Strauss S, Kück P, Ebersberger I, Walzl M, Pass G, Breuers S, Achter V, von Haeseler A, Burmester T, Hadrys H, Wägele JW, Misof B (2010) A Phylogenomic Approach to Resolve the Arthropod Tree of Life. Molecular Biology and Evolution 27: 2451–2464. doi: 10.1093/molbev/msq1
  • Mora C, Tittensor DP, Adl S, Simpson AGB, Worm B (2011) How Many Species Are There on Earth and in the Ocean? PLoS Biology 9(8): e1001127. doi: 10.1371/journal.pbio.1001127
  • Ojala D, Montoya J, Attardi G (1981) tRNA punctuation model of RNA processing in human mitochondria. Nature 290: 470–474. doi: 10.1038/290470a0
  • Pocock RI (1893) On the classification of the tracheate Arthropoda. Zoologische Anzeiger 16: 271–275.
  • Posada D, Buckley TR (2004) Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests. Systematic Biology 53: 793–808. doi: 10.1080/10635150490522304
  • Regier JC, Shultz JW, Ganley AR, Hussey A, Shi D, Ball B, Zwick A, Stajich JE, Cummings MP, Martin JW, Cunningham CW (2008) Resolving arthropod phylogeny: exploring phylogenetic signal within 41 kb of protein-coding nuclear gene sequence. Systematic Biology 57: 920–938. doi: 10.1080/10635150802570791
  • Regier JC, Shultz JW, Zwick A, Hussey A, Ball B, Wetzer R, Martin JW, Cunningham CW (2010) Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences. Nature 463: 1079–1083. doi: 10.1038/nature08742
  • Rehm P, Meusemann K, Borner J, Misof B, Burmester T (2014) Phylogenetic position of Myriapoda revealed by 454 transcriptome sequencing. Molecular Phylogenetics and Evolution 77C: 25–33. doi: 10.1016/j.ympev.2014.04.007
  • Robertson HE, Lapraz F, Rhodes AC, Telford MJ (2015) The Complete Mitochondrial Genome of the Geophilomorph Centipede Strigamia maritima. PLoS ONE 10(3): e0121369. doi: 10.1371/journal.pone.0121369
  • Ronquist F, Huelsenbeck JP. (2003) MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574. doi: 10.1093/bioinformatics/btg180
  • Rota-Stabelli O, Campbell L, Brinkmann H, Edgecombe GD, Longhorn SJ, Peterson KJ, Pisani D, Philippe H, Telford MJ (2011) A congruent solution to arthropod phylogeny: phylogenomics, microRNAs and morphology support monophyletic Mandibulata. Proceedings of the Royal Society B-Biological Sciences 278: 298–306. doi: 10.1098/rspb.2010.0590
  • Shear WA (2008) Spinnerets in the milliped order Polydesmida, and the phylogenetic significance of spinnerets in millipeds (Diplopoda). International Journal of Myriapodology 1: 123–146. doi: 10.1163/187525408X395904
  • Shelly RM (2007) Taxonomy of extant Diplopoda (Millipeds) in the modern era: Perspectives for future advancements and observations on the global diplopod community (Arthropoda: Diplopoda). Zootaxa 1668: 343−362.
  • Sierwald P, Bond JE (2007) Current status of the myriapod class diplopoda (Millipedes): Taxonomic diversity and phylogeny. Annual Review of Entomology 52: 401–420. doi: 10.1146/annurev.ento.52.111805.090210
  • Sierwald P, Shear WA, Shelley RM, Bond JE (2003) Millipede phylogeny revisited in the light of the enigmatic order Siphoniulida. Journal of Zoological Systematics and Evolutionary Research 41: 87–99. doi: 10.1046/j.1439-0469.2003.00202.x
  • Snodgrass RE (1952) A Textbook of Arthropod Anatomy. Cornell University Press, New York, 363 pp.
  • Song F, Li H, Shao R, Shi A, Bai X, Zheng X, Heiss E, Cai W (2016) Rearrangement of mitochondrial tRNA genes in flat bugs (Hemiptera: Aradidae). Scientific Reports 6: 25725. doi: 10.1038/srep25725
  • Talavera G, Castresana J (2007) Improvement of Phylogenies after Removing Divergent and Ambiguously Aligned Blocks from Protein Sequence Alignments. Systematic Biology 56: 564–577. doi: 10.1080/10635150701472164
  • Tamura K, Dudley J, Nei M, Kumar S (2007) MEGA4: molecular evolutionarygenetics analysis (MEGA) software version 4.0. Molecular Biology and Evolution 24: 1596–1599. doi: 10.1093/molbev/msm092
  • Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG (1997) The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Research 25: 4876–4882. doi: 10.1093/nar/25.24.4876
  • Woo HJ, Lee YS, Park SJ, Lim JT, Jang KH, Choi EH, Choi YG, Hwang UW (2007) Complete mitochondrial genome of a troglobite millipede Antrokoreana gracilipes (Diplopoda, Juliformia, Julida), and juliformian phylogeny. Molecules and Cells 23: 182–191.
  • Yamauchi MM, Miya M, Nishida M (2003) Complete mitochondrial DNA sequence of the swimming crab, Portunus trituberculatus (Crustacea: Decapoda: Brachyura). Gene 311: 129–135. doi: 10.1016/S0378-1119(03)00582-1