Ordinal-level phylogenomics of the arthropod class Diplopoda (millipedes) based on an analysis of 221 nuclear protein-coding loci generated using next-generation sequence analyses

Abstract

Background: The ancient and diverse, yet understudied arthropod class Diplopoda, the millipedes, has a muddled taxonomic history. Despite having a cosmopolitan distribution and a number of unique and interesting characteristics, the group has received relatively little attention; interest in millipede systematics is low compared to taxa of comparable diversity. The existing classification of the group comprises 16 orders. Past attempts to reconstruct millipede phylogenies have suffered from a paucity of characters and included too few taxa to confidently resolve relationships and make formal nomenclatural changes. Herein, we reconstruct an ordinal-level phylogeny for the class Diplopoda using the largest character set ever assembled for the group.

Methods: Transcriptomic sequences were obtained from exemplar taxa representing much of the diversity of millipede orders using second-generation (i.e., next-generation or high-throughput) sequencing. These data were subject to rigorous orthology selection and phylogenetic dataset optimization and then used to reconstruct phylogenies employing Bayesian inference and maximum likelihood optimality criteria. Ancestral reconstructions of sperm transfer appendage development (gonopods), presence of lateral defense secretion pores (ozopores), and presence of spinnerets were considered. The timings of major millipede lineage divergence points were estimated.

Results: The resulting phylogeny differed from the existing classifications in a number of fundamental ways. Our phylogeny includes a grouping that has never been described (Juliformia+Merocheta+Stemmiulida), and the ancestral reconstructions suggest caution with respect to using spinnerets as a unifying characteristic for the Nematophora. Our results are shown to have significantly stronger support than previous hypotheses given our data. Our efforts represent the first step toward obtaining a well-supported and robust phylogeny of the Diplopoda that can be used to answer many questions concerning the evolution of this ancient and diverse animal group.

Figure 1. Phylogeny recovered from Bayesian inference conducted by the program Phylobayes.

The traditional clades above the ordinal level are indicated by boxes. Support values are posterior probabilities/maximum likelihood boostrap values from the RAxML analysis.

Figure 2. Topology comparisons between our Phylobayes results and past studies.

The phylogeny obtained from phylobayes was compared to past hypotheses that were well-resolved regarding our included taxa. Three topologies met our criteria of taxon inclusion and resolution: Verhoeff [54], Enghoff [10,55], and Sierwald and Bond [2]. The full trees from each study are shown, the data source used to build the tree is indicated, and branches leading to taxa included in our analysis are shown in red. P-values obtained from each test are provided. Significant results (p < 0.05) are bolded and indicate the likelihood value of our Phylobayes results were significantly better than the alternative. Our RAxML results were also compared to the Phylobayes tree.

Figure 3. Ancestral character state reconstructions of features relating to millipede gonopods, ozopores, and spinnerets.

The phylogeny obtained from phylobayes imported into Mesquite and states were reconstructed under a parsimony model. Squares correspond to gonopods (red = none, green = ninth and tenth legs pairs, blue = eighth and ninth leg pairs, black = ambiguous). Circles correspond to the presence of millipede ozopores (red = none and blue = present). Pentagons correspond to millipede spinnerets (red = none and blue = present).

Figure 4. Chronograms representing the estimated divergence times for the lineages included in this study.

All analyses were conducted in Phylobayes using the Bayesian inference topology. A single fossil constraint was used, Pneumodesmus newmani - ~428 MYA, and the root prior was set at 510 MYA (as estimated by Rehm et al. [61]) with an exponential distribution. The orange circles indicate the fossil constraints, and the green circle indicates maximum age constraint placed on the Diplopoda (the emergence of land plants). A) divergence times estimated using the “UGAM” model; B) divergence times estimated using the “log normal” model.