Complete mitochondrial genomes of two flat-backed millipedes by next-generation sequencing (Diplopoda, Polydesmida)

Yan Dong; Lixin Zhu; Yu Bai; Yongyue Ou; Changbao Wang

doi:10.3897/zookeys.637.9909

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×10⁶ generations and sampled every 100 generations to yield a posterior probability distribution of 2×10⁴ 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; 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.

Table 1.

Download as

CSV

XLSX

Organization of the mitochondrial genome of Asiomorpha coarctata.

Feature	From	To	Length (nt)	Codons		Spacer/Overlap(-)
Feature	From	To	Length (nt)	Start Stop		Spacer/Overlap(-)
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

Table 2.

Download as

CSV

XLSX

Organization of the mitochondrial genome of Xystodesmus sp.

Feature	From	To	Length (nt)	Codons		Spacer/Overlap(-)
Feature	From	To	Length (nt)	Start Stop		Spacer/Overlap(-)
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 falciferamt 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

Posada D, Crandall KA (1998) MODELTEST: testing the model of DNA substitution. Bioinformatics 14: 817–818. doi: 10.1093/bioinformatics/14.9.817

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

Supplementary material

Supplementary material 1

Supplementary tables

Yan Dong, Lixin Zhu, Yu Bai, Yongyue Ou, Changbao Wang

Data type: molecular data

This 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.

Download file (24.65 kb)