The Complete Mitochondrial Genome of Spirobolus bungii (Diplopoda, Spirobolidae): The First Sequence for the Genus Spirobolus

Millipedes (Diplopoda) comprise one of the most important groups of large soil arthropods in terrestrial ecosystems; however, their phylogenetic relationships are poorly understood. Herein, the mitochondrial genome (mitogenome) of Spirobolus bungii was sequenced and annotated, which was 14,879 bp in size and included 37 typical mitochondrial genes (13 protein-coding genes (PCGs), two ribosomal RNA genes (rRNAs), and 22 transfer RNA genes (tRNAs)). Most of the 13 PCGs had ATN (AT/A/T/G) as the start codon except for COX1, which used CGA, and most PCGs ended with the T end codon. By comparing the gene arrangements of the mitogenomes among Diplopoda species, rearrangement occurred between and within orders. In contrast to Narceus annularus, the mitogenome genes of S. bungii had consistent orders but were transcribed in completely opposite directions, which was a novel finding in Spirobolidae. Moreover, the phylogenetic relationships within Diplopoda, which were based on the sequences of 13 PCGs, showed that S. bungii was clustered with N. annularus, followed by Abacion magmun. This indicated that there might be a close relationship between Callipodida and Spirobolida. These results could contribute to further studies on the genetics and evolutionary processes of S. bungii and other Diplopoda species.


Introduction
Millipedes Spirobolus bungii (S. bungii) belongs to the Spirobolidae family of the Diplopoda class [1]. Diplopoda comprise one of the most important groups of large soil arthropods in the terrestrial ecosystems [2], with key decomposition and nutrient cycling functions in forests [3]. They also serve as model organisms for addressing myriad evolutionary, ecological, and biological concepts and questions [4]. Diplopoda are found worldwide and reside within forests, meadows, mountains, caves, farmlands, urban green spaces, and residential areas [1]. While there have been interesting studies on millipedes in recent years, they remain a largely unexplored group, with only 12,000 of the predicted 60,000 [5] to 80,000 [6] species that are currently described. To date, there are very few studies on Diplopoda and even fewer for species in China [7,8]. Furthermore, phylogenetic studies based on morphological characteristics between diplopod taxa are rare [9,10].
Molecular data have become increasingly important in recent years. In animals, the typical mitochondrial genome (mitogenome) is a circular double-stranded DNA molecule, which encodes 13 protein-coding genes (PCGs) for the enzymes required for oxidative phosphorylation, two ribosomal RNA genes (rRNAs), and 22 transfer RNA genes (tRNAs) necessary for the translation of the proteins encoded by the mitogenome [11,12]. Compared with individual genes, the mitogenome remains a promising tool for inferring phylogenetic relationships due to its high information content. Recently, some mitogenomes in Diplopoda were published and applied to explore phylogenetic relationships [13][14][15][16]. However, only a few mitogenomes have been published for Spirobolida [17]. Further, the arrangement of genes in mitogenomes is remarkably variable across Diplopoda [13,17,18].
In this study, for the first time, the S. bungii mitogenome was assembled and characterized. The structural organization, nucleotide composition, codon usage, and AT/GC-skew were analyzed. Additionally, we conducted phylogenetic analyses based on 13 PCGs available elsewhere for the purpose of investigating the phylogenetic position of S. bungii within Diplopoda, which might further elucidate the genetics and evolutionary processes of S. bungii and other Diplopoda species.

Sample Collection and DNA Extraction
The specimens used in this study were collected from the Purple Mountain (30 • 01 N, 118 • 48 E) in 2019, where an existing deciduous broadleaved mixed forest is dominated by oaks (e.g., Quercus varialis BL, Q. accutissima Carruth), in Nanjing, Jiangsu Province, China. Following morphological identification, the samples were stored at −20°C in the Ecology Laboratory of Nanjing Forestry University (Accession No: NFU20191103). The total genomic DNA was prepared from a small portion of body segments of a single individual using the SDS-protease K-alcohol phenyl-trichlormethane method. The remaining tissue was stored at −20 • C in 90% ethanol to preserve the specimens.

Mitogenome Sequencing, Assembly, and Annotation
The complete genomic library of S. bungii was established using an Illumina HiSe-qNano DNA Sample Prep Kit (Illumina, San Diego, CA, USA), whereas the sequencing was performed using next-generation sequencing (NGS) via Illumina Hiseq2000 (Illumina, USA). To generate clean data, low-quality sequences were removed. About 40 million reads with a GC content of 43.65% were assembled to obtain a complete mitogenome using SPAdes v3.11.1 [19]. Thus, the complete mitochondrial genome sequence was used to predict the transcriptional direction of each gene component using the Improved de novo Metazoan Mitochondrial Genome Annotation (MITOS) platform [20]. The annotated mitochondrial genome sequence of S. bungii was submitted to GenBank (Accession: NC_056899.1).

Sequence Analysis
The mitochondrial ring structure was plotted by comparative genomics (CG) View Server [21], and 22 tRNA clover two-dimensional structures were predicted using tRNAscan-Se [22]. The composition skew was calculated according to the following formulae: AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) [23]. Next, a visual graph of the composition skew was created using the ggplot2 packages in R v.4.2.0. Moreover, the R script for the relative synonymous codon usage (RSCU) frequency graph was generated from PhyloSuite [24], which was then run in R v.4.2.0.

Phylogenetic Analysis
To clarify the phylogenetic position of S. bungii, the available complete mitogenomes were obtained from GenBank and were comprised of nine orders and 27 species (Table 1). Stylochyrus rarior (GenBank accession: CQ927176.2) from order Mesostigmata was used as the outgroup. A total of 27 species, including S. bungii, were employed to develop phylogenetic trees based on 13 PCGs. All operations were performed with the PhyloSuite software package [24]. The sequences were aligned in batches using MAFFT software [25]. Ambiguously aligned areas were removed using Gblocks [26]. ModelFinder was utilized to partition the codons and identify the best substitution model for phylogenetic analyses [27]. Phylogenetic trees were constructed with Bayesian inference (BI) and maximum likelihood (ML). The ML phylogenies were inferred using IQ-TREE [28] under the model automatically selected by IQ-TREE ('Auto' option in IQ-TREE) for 5000 ultrafast [29] bootstraps, as well as the Shimodaira-Hasegawa-like approximate likelihood-ratio test [30]. BI analysis was performed using MrBayes v.3.2.6 [31] with four chains (one cold chain and three hot chains).
Two independent runs of 2,000,000 generations were conducted with sampling every 100 generations. The first 25% of trees were discarded as burn-in.
Genes 2022, 13, x FOR PEER REVIEW 5 of 15 between −6 and 40 bp. Among these intergenic spacers, the longest was 17 bp (found between trnQ and trnT) ( Table 2).   The complete mitochondrial genome was 14,879 bp in size, and its overall base composition was 26.60% for A, 32.62% for T, 28.44% for G, and 12.34% for C, with a GC content of 40.78% (Table 3), which was slightly higher than other Diplopoda species (Table A1) [15,32]. The AT-skew of S. bungii was negative, while the GC-skew was positive, which was opposed to Narceus. annularus in the same family Spirobolidae ( Figure 2). Further, the GC-skews of all Polydesmida species were positive, while the AT-skews for all of this order were negative, which was completely opposed to the Spirostretida order ( Figure 2).

The PCGs
The total length of the PCGs was 10,977 bp, which was consistent with other Diplopoda species (Table A1). The base composition of the PCGs was A = 24.53%, T = 32.22%, G = 20.46%, and C = 22.78% (Table 3). In contrast to the whole mitochondrial genome, the AT-and GC-skews were both negative, which were the same as the almost Spirostreptida species (Figure 2).
The RSCU of the S. bungii mitogenome is presented in Figure 3, which indicates that Leu, Val, and Gly were the three most frequently utilized amino acids, and Cys had the lowest concentration ( Figure 3B). Nine of the twenty-two amino acids (i.e., Pro, Thr, Leu1, Arg, Ala, Ser1, Ser2, Val, and Gly) had four codons, while the others had two ( Figure 3A).

The PCGs
The total length of the PCGs was 10,977 bp, which was consistent with other Diplopoda species (Table A1). The base composition of the PCGs was A = 24.53%, T = 32.22%, G = 20.46%, and C = 22.78% (Table 3). In contrast to the whole mitochondrial genome, the AT-and GC-skews were both negative, which were the same as the almost Spirostreptida species (Figure 2).
The RSCU of the S. bungii mitogenome is presented in Figure 3, which indicates that Leu, Val, and Gly were the three most frequently utilized amino acids, and Cys had the lowest concentration ( Figure 3B). Nine of the twenty-two amino acids (i.e., Pro, Thr, Leu1, Arg, Ala, Ser1, Ser2, Val, and Gly) had four codons, while the others had two ( Figure 3A).

Transfer RNAs and Ribosomal RNAs
The typical sets of the 22 tRNAs were identified with sizes ranging from 57 bp (trnS) to 68 bp (trnQ) ( Table 2). Moreover, the total length of the tRNAs was 1375 bp, with an A+T content of 65.02%, an AT-skew of 0.056, and a GC-skew of 0.089. Among all secondary structures of the 22 tRNA genes from the S. bungii mitochondrial genome, except for trnS1, all had a typical cloverleaf structure (Figure 4), as observed in other Diplopoda mitogenomes [10].

Transfer RNAs and Ribosomal RNAs
The typical sets of the 22 tRNAs were identified with sizes ranging from 57 bp (trnS) to 68 bp (trnQ) ( Table 2). Moreover, the total length of the tRNAs was 1375 bp, with an A+T content of 65.02%, an AT-skew of 0.056, and a GC-skew of 0.089. Among all secondary structures of the 22 tRNA genes from the S. bungii mitochondrial genome, except for trnS1, all had a typical cloverleaf structure (Figure 4), as observed in other Diplopoda mitogenomes [10]. For S. bungii, the rrnL gene (length: 1270 bp) was encoded between trnV and trnL1, and the rrnS gene was 803 bp long. The total size of the two rRNAs was 2073 bp, with an A+T content of 66.14%, an AT-skew of −0.102, and a GC-skew of 0.345, which were higher than the other regions (Table A1). The rRNA AT-skews of all these species were positive, while the GC-skews were negative except for Anaulaciulus gracilipes (Figure 2). Genes 2022, 13, x FOR PEER REVIEW 9 of 15

Phylogenetic Analysis
Based on ML and BI analyses of nucleotide data of the 13 PCGs, we reconstructed the phylogenetic relationships of 26 species of Diplopoda, with S. rarior (Arachnida) as an outgroup. The two trees were similar to each other, with strongly supported branches ( Figure 5). For the BI tree, Callipodida was clustered with Sphaerotheriida and Glomeridesmida, while it did not cluster with any species for the ML tree. However, S. bungii was most closely related to N. annularus, and the relationships between Callipodida, Spirobolida, Julida, and Spirostreptida were stable, which was congruent with a previous study of mitochondrial genomes [32].
For S. bungii, the rrnL gene (length: 1270 bp) was encoded between trnV an and the rrnS gene was 803 bp long. The total size of the two rRNAs was 2073 bp, A+T content of 66.14%, an AT-skew of −0.102, and a GC-skew of 0.345, which were than the other regions (Table A1). The rRNA AT-skews of all these species were p while the GC-skews were negative except for Anaulaciulus gracilipes (Figure 2).

Phylogenetic Analysis
Based on ML and BI analyses of nucleotide data of the 13 PCGs, we reconstru phylogenetic relationships of 26 species of Diplopoda, with S. rarior (Arachnida outgroup. The two trees were similar to each other, with strongly supported b ( Figure 5). For the BI tree, Callipodida was clustered with Sphaerotheriida and G desmida, while it did not cluster with any species for the ML tree. However, S. bun most closely related to N. annularus, and the relationships between Callipodi robolida, Julida, and Spirostreptida were stable, which was congruent with a p study of mitochondrial genomes [32].

Gene Arrangement among Diplopoda Classes
By comparing the gene arrangements of the mitogenomes between Diplopoda species, rearrangement occurred between and within orders ( Figure 6). The positions of trnT for Julida differed from those of Spirobolida and Spirostreptida, which had similar gene arrangement patterns ( Figure 6). Within Julida, the positions of trnC and trnW were inversed (Figure 6), which were found in fireflies [37]. An interesting phenomenon occurred where the gene orders of the mitogenomes between S. bungii and N. annularus were consistent, while they were transcribed in completely opposite directions ( Figure 6). This was also found in the Glomeridesmidae family of the Glomeridesmida order ( Figure 6). Further, the positions of trnP in Antrokoreana gracilipes and Anaulacilus koreanus, which belonged to Julida order, were consistent but transcribed in opposite directions ( Figure 6).
By comparing the gene arrangements of the mitogenomes between Diplopoda species, rearrangement occurred between and within orders ( Figure 6). The positions of trnT for Julida differed from those of Spirobolida and Spirostreptida, which had similar gene arrangement patterns ( Figure 6). Within Julida, the positions of trnC and trnW were inversed (Figure 6), which were found in fireflies [37]. An interesting phenomenon occurred where the gene orders of the mitogenomes between S. bungii and N. annularus were consistent, while they were transcribed in completely opposite directions ( Figure 6). This was also found in the Glomeridesmidae family of the Glomeridesmida order ( Figure  6). Further, the positions of trnP in Antrokoreana gracilipes and Anaulacilus koreanus, which belonged to Julida order, were consistent but transcribed in opposite directions ( Figure 6).

Conclusions
The mitogenome of S. bungii was determined to be 14,879 bp in length, with a GC content of 40.78%. Additionally, based on a mitogenomic analysis of S. bungii, we found an intriguing phenomenon, where the AT-and GC-skews of the S. bungii mitogenome were opposed to most Diplopoda, while those of the 13 PCGs were consistent, except for Polydesmida. Consequently, this mitogenome, particularly the 13 PCGs, will assist with elucidating the genetic diversity, evolutionary origins, and genetic relationships of Diplopoda. The arrangement of genes in mitogenomes was remarkably variable across Diplopoda. Conversely, the mitogenome genes had consistent orders; however, for the Glomeridesmida and Spirobolida orders, they were transcribed in opposite directions. This indicated that the phenomenon was prevalent in Diplopoda, which will warrant additional investigations in the future. Furthermore, these results provide valuable data for the future resolution of phylogenetic relationships in this tribe.

Conclusions
The mitogenome of S. bungii was determined to be 14,879 bp in length, with a GC content of 40.78%. Additionally, based on a mitogenomic analysis of S. bungii, we found an intriguing phenomenon, where the AT-and GC-skews of the S. bungii mitogenome were opposed to most Diplopoda, while those of the 13 PCGs were consistent, except for Polydesmida. Consequently, this mitogenome, particularly the 13 PCGs, will assist with elucidating the genetic diversity, evolutionary origins, and genetic relationships of Diplopoda. The arrangement of genes in mitogenomes was remarkably variable across Diplopoda. Conversely, the mitogenome genes had consistent orders; however, for the Glomeridesmida and Spirobolida orders, they were transcribed in opposite directions. This indicated that the phenomenon was prevalent in Diplopoda, which will warrant additional investigations in the future. Furthermore, these results provide valuable data for the future resolution of phylogenetic relationships in this tribe.

Data Availability Statement:
The data that support the findings of this study are openly available in the US National Center for Biotechnology Information (NCBI database) (available online: https: //www.ncbi.nlm.nih.gov/nuccore/NC_056899.1, accessed on 22 June 2022).

Conflicts of Interest:
The authors declare no conflict of interest.