A remarkable new species of the millipede genus Trachyjulus Peters, 1864 (Diplopoda, Spirostreptida, Cambalopsidae) from Thailand, based both on morphological and molecular evidence

Abstract A new, giant species of Trachyjulus from a cave in southern Thailand is described, illustrated, and compared to morphologically closely related taxa. This new species, T. magnussp. nov., is much larger than all other congeners and looks especially similar to the grossly sympatric T. unciger Golovatch, Geoffroy, Mauriès & VandenSpiegel, 2012, which is widespread in southern Thailand. Phylogenetic trees, both rooted and unrooted, based on a concatenated dataset of the COI and 28S genes of nine species of Cambalopsidae (Trachyjulus, Glyphiulus, and Plusioglyphiulus), strongly support the monophyly of Trachyjulus and a clear-cut divergence between T. magnussp. nov. and T. unciger in revealing very high average p-distances of the COI gene (20.80–23.62%).


Introduction
In South and Southeast Asia, as well as China, the juliform millipede family Cambalopsidae Cook, 1895 is among the most diverse, common, and often highly abundant groups that clearly dominate cave millipede faunas (Golovatch 2015). Four genera are actually involved.
By far the largest genus is Glyphiulus Gervais, 1847 with its 60+ species ranging across China and Southeast Asia to Borneo in the east (Golovatch et al. 2007a(Golovatch et al. , 2007b(Golovatch et al. , 2011b(Golovatch et al. , 2011c(Golovatch et al. , 2012bJiang et al. 2017Jiang et al. , 2018Likhitrakarn et al. 2017;Liu and Wynne 2019;Golovatch and Liu 2020). The genus Plusioglyphiulus Silvestri, 1923 encompasses 28 described species ranging from northern Thailand and Laos in the west, through Myanmar and Malaysia, to Borneo in the east and southeast (Golovatch et al. , 2011aLikhitrakarn et al. 2018). Interestingly, the famous Burmese amber, 99-100 Mya, appears to contain a typical Plusioglyphiulus yet to be described (Wesener in litt.). This is evidence both of the very old age of this genus and its long presence in situ (Likhitrakarn et al. 2018).
The genus Hypocambala Silvestri, 1895 is the smallest, but particularly widespread, presently containing 14 species in Southeast Asia, as well as scattered across several islands of the Pacific and Indian oceans (Golovatch et al. 2011d).
During recent field surveys in southern Thailand, a new, unusually large Trachyjulus species was taken from a cave. From the first glance, it seemed to be particularly similar to the grossly sympatric T. unciger, but both are readily distinguished by body size and several other characters, including gonopodal structures. To better understand the species delimitations and their variations, we compare this new species to topotypes of T. unciger (Pung-Chang (= Tham Nam) Cave, Phang-Nga Province, Thailand) not only based on their morphological characters, but also on molecular evidence. In addition, molecular-based phylogenetic relationships within the genus Trachyjulus are revealed and discussed for the first time using mitochondrial cytochrome c oxidase subunit I (COI) and nuclear gene 28S rRNA sequences. These were obtained from para-or topotypes of nine species of Cambalopsidae, including not only five Trachyjulus, but also two Glyphiulus and two Plusioglyphiulus as outgroups. Two members of the family Harpagophoridae from the same order Spirostreptida, as well as two of the family Julidae, order Julida, are also included as more distant outgroups for tree rooting.

Sample collection
Specimens were collected from southern Myanmar and southern Thailand under the Animal Care and Use Protocol Review No. 1723018. The collecting sites were located by GPS by using a Garmin GPSMAP 60 CSx, and all coordinates and elevations were rechecked with Google Earth. Photographs of live animals were taken using a Nikon 700D digital camera with a Nikon AF-S VR 105 mm macro lens. The specimens collected were euthanized by a two-step method following AVMA Guidelines for the Euthanasia of Animals (AVMA 2013). Specimens were then preserved in 95% ethanol for morphological and molecular studies. Ethanol was replaced after 24 hours with fresh 95% ethanol to prevent their defensive chemicals from affecting future DNA extraction. Mostly para-or topotypes of six described species were also used for molecular analyses (Table 1).
The holotype, as well as most of the paratypes are housed in the Museum of Zoology, Chulalongkorn University (CUMZ), Bangkok, Thailand; a few paratypes have also been donated to the collections of the Zoological Museum, State University of Moscow, Russia (ZMUM) and the Natural History Museum of Denmark, University of Copenhagen, Denmark (NHMD), as indicated in the text.

Morphological study
The specimens were examined, measured, and photographed under a Nikon SMZ 745T trinocular stereo microscope equipped with a Canon EOS 5DS R digital SLR camera. Scanning electron micrographs (SEM) were taken with a JEOL, JSM-5410 LV microscope using gold-coated samples, and the material returned to alcohol upon examination. Digital images obtained were processed and edited with Adobe Photoshop CS5. Line drawings were executed based on photographs and specimens examined under a Nikon SMZ 745T trinocular stereo microscope, equipped with a Canon EOS 5DS R digital SLR camera. The terminology used and the carinotaxic formulae in the descriptions follow those in Golovatch et al. (2007aGolovatch et al. ( , 2007bGolovatch et al. ( , 2012a, while body ring counts are after Enghoff et al. (1993) and Golovatch et al. (2007a).

DNA extraction and molecular identification
Total genomic DNA was extracted from the dissected midbody ring tissues using the DNA extraction kit for animal tissue (NucleoSpin Tissue extraction kit, Macherey-Nagel, Germany), following the standard procedure of the manual. Fragments of the mitochondrial cytochrome c oxidase subunit I (COI, 690 bp) gene were amplified using LCO1490 (5'-GGTCAACAAATCATAAAGATATTGG-3'; Folmer et al. 1994) and hco outout (5'-GTAAATATATGRTGDGCTC; Schulmeister et al. 2002) or Nan-  Bogdanowicz et al. 1993); while fragments of the nuclear 28S ribosomal RNA large subunit gene (28S) were amplified using primers 28F2-2 (5'-GCAGAACTGGCGCTGAGGGATGAAC-3') and 28SR2 GAGGCTGTKCACCTTGGAGACCTGCTGCG-3'; Passamaneck et al. 2004). The PCR amplification was performed using a T100™ thermal cycler (BIO-RAD) with a final reaction volume of 20 μL (15 μL of EmeraldAmp GT PCR Master Mix, 1.5 μL of each primer, 10 ng of template DNA and distilled water up to 20 μL total volume). Thermal cycling was performed at 94 °C for 3 min, followed by 35 cycles of 94 °C for 30 s, annealing at 42-56 °C (depending on samples and the primer paired) for 60 s, extension at 72 °C for 90 s, and a final extension at 72 °C for 5 min. Amplification of PCR products were confirmed through 1.5% (w/v) agarose gel electrophoresis before purification by PEG precipitation. Purified PCR products were sequenced in both directions (forward and reverse) using an automated sequencer (ABI prism 3730XL). All nucleotide sequences in this study were deposited in the GenBank Nucleotide sequences database under submission numbers MN893771-MN893781 for COI, and MN897820-MN897826 for 28S. The collecting localities and submission codes of each nominal species are listed in Table 1.
The sequences were edited and aligned using Clustal W, implemented in MEGA7 (Kumar et al. 2016). The aligned sequences were estimated for the best-fit model of nucleotide substitution for each gene separately by KAKUSAN4 (Tanabe 2011). Two phylogenetic methods, maximum likelihood (ML) and Bayesian Inference (BI), were implemented through the on-line CIPRES Science Gateway (Miller et al. 2010). The ML analysis was performed using RAxML v.8.2.10 (Stamatakis 2014) with 1,000 bootstrap replications and GTRGAMMA as the nucleotide substitution model (Silvestro and Michalak 2012). The BI analysis was performed by MrBayes 3.2.6 (Ronquist et al. 2012) using the Markov chain Monte Carlo technique (MCMC). The best-fit evolution models based on the Akaike Information Criterion (AIC: Akaike 1974) were applied: SYM+G for the 1 st COI codon, and HKY85+G for the 2 nd COI codon, the 3 rd COI codon, and the 28S gene. Ten million generations were run with a random starting tree. The resultant trees were sampled every 1,000 th generation and were used to estimate the consensus tree topology; bipartition posterior probability (bpp) and branch lengths after the first 25% of obtained trees were discarded as burn-in. All Effective Sample Size (ESS) values sampled from the MCMC analysis were greater than 1,000 in all parameters. A neighbour-joining tree (NJ) based on K2P-distance was constructed based on the amino acid alignment of peptide sequences corresponding to the mitochondrial COI dataset. Interspecific genetic divergences based on the COI sequence were also evaluated using uncorrected p-distances. The NJ tree and p-distance were implemented in MEGA7 (Kumar et al. 2016 Paratypes. 15 ♂, 20 ♀ (CUMZ), 1 ♂, 1 ♀ (ZMUM), 1 ♂, 1 ♀ (NHMD), same locality, together with holotype.
Name. To emphasize the largest body size of this species compared to all other species known in the genus.
Coloration of live animals red-brown to yellow-brownish (Fig. 1), venter and legs brownish yellow to yellowish, antennae light to pale yellowish, eyes blackish, a thin axial line traceable; coloration in alcohol, after one year of preservation, similar, but body yellowbrownish to light brownish, vertex red-brown to light brown, eyes blackish to brownish.
Legs short, on midbody rings about 2/3 (♂, ♀) as long as body height (Figs 2I,J). Claw at base with a strong accessory spiniform claw almost half as long as main claw (Fig. 2K). Tarsi and tarsal setae very delicately fringed.
Remark. The often striking colour difference between head+collum+ring 2 and the remaining rings observed in SEM micrographs (Fig. 2) is certainly an artifact resulting from unwanted electrical charging.

Phylogenetic analysis
Our concatenated dataset contained 15 individuals, including seven Trachyjulus ingroup and eight outgroup species, and an alignment of approximately 1,501 base pairs (bp). We were unable to obtain sequences of the 28S gene from T. phylloides, G. sattaa, and P. erawan. The final alignment of the COI gene fragment yielded 690 bp (298 variable sites, 270 parsimony informative), while the 28S gene fragment comprised 811 bp (100 variable sites, 45 parsimony informative). The phylogenetic tree estimated by both ML and BI revealed equivalent topologies. As only one position within the outgroup taxa was controversial, solely a ML tree is shown in Figure  5A. The monophyly of the genus Trachyjulus was strongly supported (1 bpp for BI and 96% bootstrap values for ML). Within the Trachyjulus clade, T. singularis was placed in the basal part, followed by T. magnus sp. nov., T. unciger, and a sister clade of T. phylloides and T. bifidus, respectively. All internal nodes were strongly supported (0.99-1 bpp for BI and 97-100% bootstrap values for ML). In addition, nine cambalopsid species were recovered as a monophyletic clade against the analyzed outgroups representing two other families and one order, although only the BI analysis was supported by this grouping and showed a bpp of 0.95. Within the cambalopsid clade, three genera were clustered separately as a monophyletic clade. However, no evolutionary relationship among them was revealed. The NJ tree based on the COI corresponding amino acid sequences also clearly recovered the monophyly of Trachyjulus (73% bootstrap values) (Fig. 5B).
The interspecific divergence of the COI uncorrected p-distance among these nine cambalopsid species was found to be generally high (13.48-24.49%; Table 2). Among the Trachyjulus species concerned, the average distance values ranged from 15.07-23.62%. Trachyjulus singularis showed the highest divergence from the other Trachyjulus species, ranging from 21.16-23.62%. The lowest divergence among Trachyjulus species was 15.07% between T. bifidus and T. phylloides. Five Trachyjulus species showed a long-distance relationship to their closely related genera, Glyphiulus (18.84-23.62%) and Plusioglyphiulus (20.00-24.49%). In addition, the average distances between the members of Glyphiulus and Plusioglyphiulus were also relatively high, ranging from 17.39-21.16%. The interspecific divergence among Glyphiulus and Plusioglyphiulus species amounted to 17.97 and 13.48, respectively.

Discussion
Trachyjulus magnus sp. nov. clearly represents a taxonomically valid species based on both morphological and molecular evidence. In the latest taxonomic review of Trachyjulus, Golovatch et al. (2012a) emphasized and listed the following primary morphological characters deemed useful to distinguishing it from the other cambalopsid genera. The genus Trachyjulus shows a collum which is smooth or nearly smooth at least dorsally, usually not particularly inflated compared to postcollum constrictions; midbody metazonae are strongly carinate, the carinotaxic formulae typically being 11-8/11-8+I/ i+2/2+m/m; male leg 1 is strongly reduced to a broad transverse coxosternum that shows a pair of central, often completely fused coxal processes flanked by rudimentary telopodites; some structures of the gonopods are also unique. It is gonopodal structures, often highly conservative, that usually appear to be especially useful for species delimitations among congeners in the family Cambalopsidae (Golovatch et al. 2007a(Golovatch et al. , 2007b(Golovatch et al. , 2011a(Golovatch et al. , 2011b(Golovatch et al. , 2011c(Golovatch et al. , 2011d(Golovatch et al. , 2012bJiang et al. 2017;Likhitrakarn et al. 2017). Morphologically, the new species looks especially similar to T. unciger, but both are clearly distinguishable (see Diagnosis above). Molecular evidence likewise reveals a sufficiently strong genetic divergence between T. magnus sp. nov. and T. unciger (p-distance = 20.43±1.52) ( Table 2). Compared to other studies, the interspecific distances among Bavarian millipedes range from >5% among members in the same genus and up to 33.18% between the different orders, averaging 14.17% (Spelda et al. 2011). In Thailand, Pimvichai et al. (2016) reported the interspecific divergences of mitochondrial COI as ranging between 2 and 17% in the large-bodied Thyropygus millipedes, family Harpagophoridae. The average interspecific divergences among Trachyjulus species in the present study appear to be even higher than in Thyropygus: 15.07-23.62%. This may be accounted for by the much smaller sizes of Trachyjulus spp., as well as their usually more limited dispersal capacities that make them largely restricted to a particular cave or cave complex (Golovatch et al. 2007a(Golovatch et al. , 2011a. High rates of interspecific genetic differentiation in small-bodied cave-dwelling species have long been reported elsewhere: 8.2-9.2% between two parapatric Callipodida millipedes from the USA, Tetracion tennesseensis Causey, 1959and T. jonesi Hoffman, 1956(Loria et al. 2011. The phylogenetic trees, both rooted and inrooted, and based on the concatenated dataset, provide strong support to the monophyly of the genus Trachyjulus in both ML and BI analyses (bpp = 1.0 for BI and bootstrap value = 100% for ML) (Fig. 5A). In addition, the NJ tree based on the COI corresponding amino acid sequences in which the protein evolves at a slow rate (Drummond et al. 2005) clearly recovered the monophyly of Trachyjulus as well (73% bootstrap values) (Fig. 5B). Therefore, the molecular evidence confirms that all of the Trachyjulus species concerned, including the new species, do belong to the same genus. Because the unrooted phylogram totally failed to alter the topology of the rooted one, only the latter is shown in Figure 5.
Trachyjulus singularis was recovered as the basal clade of the tree. It also showed the highest genetic divergence from the other Trachyjulus species (21.16-23.62%). These results are in accordance with their geographic distributions, as T. singularis occurs only in eastern Thailand, i.e. far away from the congeners in southern Thailand. In addition, T. singularus has retained the ancestral character of a divided promentum of the gnathochilarium, a trait absent from the other members of Trachyjulus, but present in two other related genera, Glyphiulus and Plusioglyphiulus.
In conclusion, we put on record the first results of a molecular phylogenetic study on Trachyjulus, a largely cavernicolous genus, using a combination of the nuclear 28S rRNA and mitochondrial CO1 genes for a total of 1,501 bp. Our results reveal high rates of interspecific divergence among Trachyjulus species and other closely related genera. Given that Thailand and the neighbouring countries are extremely rich in karst and karst caves, there can hardly be any doubt that additional new species of Cambalopsidae generally and Trachyjulus in particular still await discovery. A combination of morphological and molecular studies in Cambalopsidae seems the best to provide further insights into the taxonomy and phylogenetic relationships in this large and widespread group.