2urn:lsid:arphahub.com:pub:45048D35-BB1D-5CE8-9668-537E44BD4C7Eurn:lsid:zoobank.org:pub:91BD42D4-90F1-4B45-9350-EEF175B1727AZooKeysZK1313-29891313-2970Pensoft Publishers10.3897/zookeys.637.99099909Research ArticleAnimaliaArthropodaDiplopodaInvertebrataMyriapodaPolydesmidaEvolutionary biologyGeneticsMolecular geneticsMolecular systematicsChinaComplete mitochondrial genomes of two flat-backed millipedes by next-generation sequencing (Diplopoda, Polydesmida)DongYandongyan_bio@126.com1ZhuLixin1BaiYu1OuYongyue1WangChangbao1College of Biology and Food Engineering, Chuzhou University, Chuzhou 239000, ChinaMolecular Biology Laboratory, College of Biology and food technology, 1528 Fengle road, Chuzhou University, Chuzhou, Anhui, China, 239000ChuzhouChina
Corresponding author: Yan Dong (dongyan_bio@126.com)
Academic editor: Pavel Stoev
2016281120166371205724FFEC-FF8F-FFEF-FFEC-3215FFA81E09A53FB0CA-30EC-4B71-AFD6-D00AE30218F71789401707201617112016Yan Dong, Lixin Zhu, Yu Bai, Yongyue Ou, Changbao WangThis is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.http://zoobank.org/A53FB0CA-30EC-4B71-AFD6-D00AE30218F7
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, Asiomorphacoarctata 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.
Dong Y, Zhu L, Bai Y, Ou Y, Wang C (2016) Complete mitochondrial genomes of two flat-backed millipedes by next-generation sequencing (Diplopoda, Polydesmida). ZooKeys 637: 1–20. https://doi.org/10.3897/zookeys.637.9909
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 Cermatobiuslongicornis and the millipede Prionobelum sp. is identical to that of Limuluspolyphemus (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 Strigamiamaritima (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 LimulusPolyphemus 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: Appalachioriafalcifera (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 Asiomorphacoarctata 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.
MethodsTaxon sampling, DNA extraction and PCR
Specimens of Asiomorphacoarctata 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 (Narceusannularus, Thyropygus sp., Antrokoreanagracilipes, Appalachioriafalcifera, Abacionmagnum, Brachycybelecontii, 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.
ResultsMitochondrial genome organization
The complete mt genome of Asiomorphacoarctata (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.
Mitochondrial genomes of the two millipedes sequenced in this study. AAsiomorphacoarctataBXystodesmus 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; cox1–3 for cytochrome oxidase subunits 1–3; cob for cytochrome b, nad1–6 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.
https://binary.pensoft.net/fig/112690
Organization of the mitochondrial genome of Asiomorphacoarctata.
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. Antrokoreanagracilipes) include two major non-coding regions and others contain a single non-coding region, such as Narceusannularus, Prionobelum sp. and Appalachioriafalcifera (Figure 5). Of the genomes sequenced here, Asiomorphacoarctata 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).
Sequences of the non-coding region in Asiomorphacoarctata, primary structures of tandemly repeated regions (11.4 × 38 bp).
https://binary.pensoft.net/fig/112691Transfer RNA
There are 22 potential tRNA genes in Asiomorphacoarctata 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.
Putative secondary structures of the 22 tRNA genes of Asiomorphacoarctata. Watson-Crick base-pairing is indicated by solid lines, and G–T pairs are indicated with plus signs.
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.
https://binary.pensoft.net/fig/112693Gene 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 Limuluspolyphemus (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 Appalachioriafalcifera 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.
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.
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.
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.
https://binary.pensoft.net/fig/112695
Within the Diplopoda, the Polydesmida clade, including Appalachioriafalcifera, Asiomorphacoarctata, and Xystodesmus sp. is well supported. A clade consisting of Polydesmida plus Colobognatha (Brachycybelecontii) 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; Antrokoreanagracilipes, GenBank accession NC_010221; Appalachioriafalcifera, 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 Appalachioriafalcifera. 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 Appalachioriafalciferamt 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).
ReferencesAbascalFPosadaDZardoyaR (2007) MtArt: a new model of amino acid replacement for Arthropoda.24: 1–5. doi: 10.1093/molbev/msl136AguadoMTGrandeCGerthMBleidornCNoreñaC (2016) Characterization of the complete mitochondrial genomes from Polycladida (Platyhelminthes) using next-generation sequencing.575: 199–205. doi: 10.1016/j.gene.2015.08.054AxP (1999) Gustav Fischer Verlag, Stuttgart, 384 pp.BäckerHFanenbruckMWägeleJW (2008) A forgotten homology supporting the monophyly of Tracheata: the subcoxa of insects and myriapods.247: 185–207. doi: 10.1016/j.jcz.2007.11.002BerntMDonathAJühlingFExternbrinkFFlorentzCFritzschGPützJMiddendorfMStadlerPF (2013) MITOS: improved de novo metazoan mitochondrial genome annotation.69(2): 313–319. doi: 10.1016/j.ympev.2012.08.023BitschCBitschJ (2004) Phylogenetic relationships of basal hexapods among mandibulate arthropods: a cladistic analysis based on comparative morphological characters.33: 511–550. doi: 10.1111/j.0300-3256.2004.00162.xBlankeaAWesenerT (2014) Revival of forgotten characters and modern imaging techniques help to produce a robust phylogeny of the Diplopoda (Arthropoda, Myriapoda).43(1): 63–75. doi: 10.1016/j.asd.2013.10.003BooreJL (1999) Animal mitochondrial genomes.27: 1767–1780. doi: 10.1093/nar/27.8.1767BooreJLLavrovDBrownWM (1998) Gene translocation links insects and crustaceans.392(6677): 667–668. doi: 10.1038/33577BoudreauxHB (1979) Significance of inter segmental tendon system in arthropod phylogeny and a monophyletic classification of Arthropoda. In: GuptaAP (Ed.) . Van Nostrand Reinhold Company, New York, 551–586.BrewerMSSierwaldPBondJE (2012) Millipede Taxonomy after 250 Years: Classification and Taxonomic Practices in a Mega-Diverse yet Understudied Arthropod Group.7(5): e37240. doi: 10.1371/journal.pone.0037240BrewerMSSwaffordLSpruillCLBondJE (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.8(7): e68005. doi: 10.1371/journal.pone.0068005BuddGETelfordMJ (2009) The origin and evolution of arthropods.457: 812–817. doi: 10.1038/nature07890CameronSL (2014) Insect mitochondrial genomics: implications for evolution and phylogeny.59: 95–117. doi: 10.1146/annurev-ento-011613-162007CarapelliASoto-AdamesFNSimonCFratiFNardiFDallaiR (2004) Secondary structure, high variability and conserved motifs for domain III of 12S rRNA in the Arthropleona (Hexapoda; Collembola).13: 659–670. doi: 10.1111/j.0962-1075.2004.00528.xCastresanaJ (2000) Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis.17(4): 540–552. doi: 10.1093/oxfordjournals.molbev.a026334ChenWJBuYCarapelliADallaiRLiSYinWYLuanYX (2011) The mitochondrial genome of Sinentomonerythranum (Arthropoda: Hexapoda: Protura): an example of highly dviergent evolution.11: 246. doi: 10.1186/1471-2148-11-246ChenWJKochMMallattJMLuanYX (2014) Comparative analysis of mitochondrial genomes in Diplura (Hexapoda, Arthropoda): Taxon sampling is crucial for phylogenetic inferences.6(1): 105–120. doi: 10.1093/gbe/evt207ChilanaPSharmaARaiA (2012) Insect genomic resources: status, availability and future.102: 571–580.CoatesBS (2014) Assembly and annotation of full mitochondrial genomes for the corn rootworm species, Diabroticavirgiferavirgifera and Diabroticabarberi (Insecta: Coleoptera: Chrysomelidae), using next generation sequence data.542: 190–197. doi: 10.1016/j.gene.2014.03.035DohleW (1980) Sind die Myriapoden eine monophyletische Gruppe? Eine Diskussionder Verwandtschaftsbeziehungen der Antennaten. Abhandl. Verhandl. Naturwiss.23: 45–104.DohleW (1998) Myriapod–insect relationships as opposed to an insect–crustacean sister group relationship. In: ForteyRAThomasRH (Eds) . Chapman and Hall, London, 305–315. doi: 10.1007/978-94-011-4904-4_23DongYSunHYGuoHPanDQianCYHaoSJ (2012a) The complete mitochondrial genome of Pauropuslongiramus (Myriapoda: Pauropoda): implications on early diversification of the myriapods revealed from comparative analysis.505(1): 57–65. doi: 10.1016/j.gene.2012.05.049DongYXuJJHaoSJSunHY (2012b) The complete mitochondrial genome of the giant pill millipede,Sphaerotheriidae sp. (Myriapoda: Diplopoda: phaerotheriida).23(5): 333–335. doi: 10.3109/19401736.2012.683184EdgecombeGD (2004) Morphological data, extant Myriapoda, and the myriapod stem-group.73: 207–252.EdgecombeGD (2010) Arthropod phylogeny: an overview from the perspectives of morphology, molecular data and the fossil record.39(2–3): 74–87. doi: 10.1016/j.asd.2009.10.002EdgecombeGD (2011) Phylogenetic relationships of Myriapoda. In: MinelliA (Ed.) . Brill, Leiden, 1–20. doi: 10.1163/9789004188266_002EdgecombeGD (2015) Diplopoda - Phylogenetic relationships. In: MinelliA (Ed.) , Volume 2, Chapter 2. Brill, Leiden, 353–362. doi: 10.1163/9789004188273_016EnghoffHGolovatchSShortMStoevPWesenerT (2015) Diplopoda – Taxonomic overview. In: MinelliA (Ed.) , Volume 2, Chapter 2. Brill, Leiden, 363–453. doi: 10.1163/9789004188273_017EnghoffHDohleWBlowerJG (1993) Anamorphosis in millipedes (Diplopoda): the present state of knowledge with some developmental and phylogenetic considerations.109: 103–234. doi: 10.1111/j.1096-3642.1993.tb00305.xFernándezREdgecombeGDGiribetG (2016) Exploring phylogenetic relationships within Myriapoda and the effects of matrix composition and occupancy on phylogenomic reconstruction.65(5): 871–89. doi: 10.1093/sysbio/syw041GaiYMaHSunXMaJLiCYangQ (2013) The complete mitochondrial genome of Cermatobiuslongicornis (Chilopoda: Lithobiomorpha: Henicopidae).24(4): 331–332. doi: 10.3109/19401736.2012.760078GaiYSongDSunHYangQZhouK (2008) The complete mitochondrial genome of Symphylella sp. (Myriapoda: Symphyla): Extensive gene order rearrangement and evidence in favor of Progoneata.49(2): 574–85. doi: 10.1016/j.asd.2009.10.002GaiYHSongDXSunHYZhouKY (2006) Myriapod monophyly and relationships among myriapod classes based on nearly complete 28S and 18S rDNA sequences.23: 1101–1108. doi: 10.2108/zsj.23.1101GissiCIannelliFPesoleG (2008) Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of congeneric species.101: 301–320. doi: 10.1038/hdy.2008.62HoffmanRGolovatchSAdisJde MoraisJ (2002) Diplopoda. In: AdisJ (Ed.) . Pensoft, Sofia-Moscow, 505–533.KirknessEFHaasBJSunWBraigHRPerottiMAClarkJMLeeSHRobertsonHMKennedyRCElhaikEGerlachDKriventsevaEVElsikCGGraurDHillCAVeenstraJAWalenzBTubíoJMRibeiroJMRozasJJohnstonJSReeseJTPopadicATojoMRaoultDReedDLTomoyasuYKrausEMittapalliOMargamVMLiHMMeyerJMJohnsonRMRomero-SeversonJVanzeeJPAlvarez-PonceDVieiraFGAguadéMGuirao-RicoSAnzolaJMYoonKSStrycharzJPUngerMFChristleySLoboNFSeufferheldMJWangNDaschGAStruchinerCJMadeyGHannickLIBidwellSJoardarVCalerEShaoRBarkerSCCameronSBruggnerRVRegierAJohnsonJViswanathanLUtterbackTRSuttonGGLawsonDWaterhouseRMVenterJCStrausbergRLBerenbaumMRCollinsFHZdobnovEMPittendrighBR (2010) Genome sequences of the human body louse and its primary endosymbiont provide insights into the permanent parasitic lifestyle.107(27): 12168–12173. doi: 10.1073/pnas.1003379107KnausBJCronnRListonAPilgrimKSchwartzMK (2011) Mitochondrial genome sequences illuminate maternal lineages of conservation concern in a rare carnivore.11: 10. doi: 10.1186/1472-6785-11-10LaslettDCanbackB (2008) ARWEN: a program to detect tRNA genes in metazoan mitochondrial nucleotide sequences.24: 172–175. doi: 10.1093/bioinformatics/btm573LavrovDBooreJBrownW (2000a) The complete mitochondrial DNA sequence of the Horseshoe Crab Limuluspolyphemus.17(5): 813–824. doi: 10.1093/oxfordjournals.molbev.a026360LavrovDBrownWBooreJ (2000b) A novel type of RNA editing occurs in the mitochondrial tRNAs of the centipede Lithobiusforficatus.97: 13738–13742. doi: 10.1073/pnas.250402997LoweTMEddySR (1997) tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence.25: 955–964. doi: 10.1093/nar/25.5.0955MaPFGuoZHLiDZ (2012) Rapid sequencing of the bamboo mitochondrial genome using Illumina technology and parallel episodic evolution of organelle genomes in grasses.7: e30297. doi: 10.1371/journal.pone.0030297McknightMLShafferHB (1997) Large, rapidly evolving intergenic spacers in the mitochondrial DNA of the Salamander family Ambystomatidae (Amphibia: Caudata).14: 1167–1176. doi: 10.1093/oxfordjournals.molbev.a025726MeusemannKvon ReumontBMSimonSRoedingFStraussSKückPEbersbergerIWalzlMPassGBreuersSAchterVvon HaeselerABurmesterTHadrysHWägeleJWMisofB (2010) A Phylogenomic Approach to Resolve the Arthropod Tree of Life.27: 2451–2464. doi: 10.1093/molbev/msq1MoraCTittensorDPAdlSSimpsonAGBWormB (2011) How Many Species Are There on Earth and in the Ocean?9(8): e1001127. doi: 10.1371/journal.pbio.1001127OjalaDMontoyaJAttardiG (1981) tRNA punctuation model of RNA processing in human mitochondria.290: 470–474. doi: 10.1038/290470a0PocockRI (1893) On the classification of the tracheate Arthropoda.16: 271–275.PosadaDBuckleyTR (2004) Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests.53: 793–808. doi: 10.1080/10635150490522304PosadaDCrandallKA (1998) MODELTEST: testing the model of DNA substitution.14: 817–818. doi: 10.1093/bioinformatics/14.9.817RegierJCShultzJWGanleyARHusseyAShiDBallBZwickAStajichJECummingsMPMartinJWCunninghamCW (2008) Resolving arthropod phylogeny: exploring phylogenetic signal within 41 kb of protein-coding nuclear gene sequence.57: 920–938. doi: 10.1080/10635150802570791RegierJCShultzJWZwickAHusseyABallBWetzerRMartinJWCunninghamCW (2010) Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences.463: 1079–1083. doi: 10.1038/nature08742RehmPMeusemannKBornerJMisofBBurmesterT (2014) Phylogenetic position of Myriapoda revealed by 454 transcriptome sequencing.77C: 25–33. doi: 10.1016/j.ympev.2014.04.007RobertsonHELaprazFRhodesACTelfordMJ (2015) The Complete Mitochondrial Genome of the Geophilomorph Centipede Strigamiamaritima.10(3): e0121369. doi: 10.1371/journal.pone.0121369RonquistFHuelsenbeckJP. (2003) MRBAYES 3: Bayesian phylogenetic inference under mixed models.19: 1572–1574. doi: 10.1093/bioinformatics/btg180Rota-StabelliOCampbellLBrinkmannHEdgecombeGDLonghornSJPetersonKJPisaniDPhilippeHTelfordMJ (2011) A congruent solution to arthropod phylogeny: phylogenomics, microRNAs and morphology support monophyletic Mandibulata.278: 298–306. doi: 10.1098/rspb.2010.0590ShearWA (2008) Spinnerets in the milliped order Polydesmida, and the phylogenetic significance of spinnerets in millipeds (Diplopoda).1: 123–146. doi: 10.1163/187525408X395904ShellyRM (2007) Taxonomy of extant Diplopoda (Millipeds) in the modern era: Perspectives for future advancements and observations on the global diplopod community (Arthropoda: Diplopoda).1668: 343−362.SierwaldPBondJE (2007) Current status of the myriapod class diplopoda (Millipedes): Taxonomic diversity and phylogeny.52: 401–420. doi: 10.1146/annurev.ento.52.111805.090210SierwaldPShearWAShelleyRMBondJE (2003) Millipede phylogeny revisited in the light of the enigmatic order Siphoniulida.41: 87–99. doi: 10.1046/j.1439-0469.2003.00202.xSnodgrassRE (1952) Cornell University Press, New York, 363 pp.SongFLiHShaoRShiABaiXZhengXHeissECaiW (2016) Rearrangement of mitochondrial tRNA genes in flat bugs (Hemiptera: Aradidae).6: 25725. doi: 10.1038/srep25725TalaveraGCastresanaJ (2007) Improvement of Phylogenies after Removing Divergent and Ambiguously Aligned Blocks from Protein Sequence Alignments.56: 564–577. doi: 10.1080/10635150701472164TamuraKDudleyJNeiMKumarS (2007) MEGA4: molecular evolutionarygenetics analysis (MEGA) software version 4.0.24: 1596–1599. doi: 10.1093/molbev/msm092ThompsonJDGibsonTJPlewniakFJeanmouginFHigginsDG (1997) The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools.25: 4876–4882. doi: 10.1093/nar/25.24.4876WooHJLeeYSParkSJLimJTJangKHChoiEHChoiYGHwangUW (2007) Complete mitochondrial genome of a troglobite millipede Antrokoreanagracilipes (Diplopoda, Juliformia, Julida), and juliformian phylogeny.23: 182–191.YamauchiMMMiyaMNishidaM (2003) Complete mitochondrial DNA sequence of the swimming crab, Portunustrituberculatus (Crustacea: Decapoda: Brachyura).311: 129–135. doi: 10.1016/S0378-1119(03)00582-1Supplementary materials10.3897/zookeys.637.9909.suppl190384484E5410B-A47B-5B1C-AC49-1F3D25031689
Supplementary tables
molecular data
https://binary.pensoft.net/file/112696This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.Yan Dong, Lixin Zhu, Yu Bai, Yongyue Ou, Changbao Wang