A revision of the Thyropygus allevatus group. Part V: Nine new species of the extended opinatus subgroup, based on morphological and DNA sequence data (Diplopoda: Spirostreptida: Harpagophoridae)

. The Thyropygus opinatus subgroup (Diplopoda: Harpagophoridae) of the T. allevatus group in Thailand is revised. Based on a phylogenetic analysis of mtDNA sequence data, it is merged with the T. bifurcus subgroup to form an extended T. opinatus subgroup. Nine new species are described: Thyropygus cimi sp. nov. and T. forceps sp. nov. from Nakhonsrithammarat Province, T. culter sp. nov., T. planispina sp. nov., T. undulatus sp. nov. and T. ursus sp. nov. from Krabi Province, T. mesocristatus sp. nov. from Songkhla Province, T. navychula sp. nov. from Phang-Nga Province and T. sutchariti sp. nov. from Phetchaburi Province.


Introduction
The millipede genus Thyropygus Pocock, 1894, is the most species-rich genus of the family Harpagophoridae (subfamily Harpagophorinae) in Southeast Asia.Currently, it includes 56 named species of which 35 occur in Thailand (Enghoff 2005;Jeekel 2006;Pimvichai et al. 2009a, b;2011a, b).The genus is widely distributed in Southeast Asia and is particularly diverse in Thailand.In a series of previous papers four subgroups of the large Thyropygus allevatus group were revised: the T. opinatus subgroup (Pimvichai et al. 2009a), the T. bifurcus subgroup (Pimvichai et al. 2009b), the T. induratus subgroup (Pimvichai et al. 2011a), and the T. cuisinieri subgroup (Pimvichai et al. 2011b).Yet, several species of the T. allevatus group still have not been assigned to a subgroup.
Although the four currently recognized subgroups of the T. allevatus group appear morphologically distinct, Pimvichai et al. (2014) observed that the T. opinatus and T. bifurcus subgroups did not form separate clades in a mtDNA phylogeny, but instead formed a single mixed clade.This suggested that the two subgroups should not be kept separate, which would also better refl ect their morphological similarity (Pimvichai et al. 2009b).
The mitochondrial cytochrome c oxidase subunit I (COI) is widely used for DNA-based species identifi cation (DNA barcoding) of closely related species (Hebert et al. 2003), while mitochondrial 16S rRNA sequences are often used as an additional marker in diplopod taxonomy and phylogeny (Pimvichai et al. 2014).Hence, in the present study, COI and 16S rRNA sequence data are used to (1) delimit species of the T. opinatus and T. bifurcus subgroups, and (2) assess the monophyly of both these subgroups.As a result we describe nine new species and suggest to abandon the concept of a separate T. bifurcus subgroup.Instead, we propose to merge the T. bifurcus and T. opinatus subgroups into an extended T. opinatus subgroup.

Material and methods
Fresh specimens were hand-collected and preserved partly in 70% ethanol, partly in a freezer at -20 °C for subsequent DNA studies.Specimens were examined from the following collections: CUMZ = Museum of Zoology, Chulalongkorn University, Bangkok, Thailand NHMW = Naturhistorisches Museum, Vienna, Austria ZMUC = Zoological Museum, Natural History Museum of Denmark, University of Copenhagen Partial COI and 16S rRNA gene sequences were used to construct a phylogenetic tree of Thyropygus species.COI was also used to evaluate genetic divergences between species.The DNA sequence data were extracted from Pimvichai et al. (2014) and supplemented with sequences of seven new samples.The procedures for extracting and analyzing COI and 16S rRNA sequences followed Pimvichai et al. (2014).The COI fragment was amplifi ed using the primers LCO-1490 and HCO-2198 (Folmer et al. 1994), while the 16S rRNA DNA fragment was amplifi ed with the universal primers 16Sar and 16Sbr (Kessing et al. 1989).The new sequences are deposited in GenBank under the accession numbers KU306518-KU306531.Collection localities and accession codes for each nominal species are shown in Table 1.
Table 1.Specimens from which partial 16S rRNA and/or COI DNA sequences were obtained.The cuisinieri and induratus subgroups are as defi ned by Pimvichai et al. (2011a, b).Names of provinces in Thailand are shown in capitals.Abbreviated locality names are provided in parentheses after species names.GenBank accession numbers for DNA sequences are indicated for each species.(Table continued on next pages.) Table 1 (concluded).

Alignment and phylogenetic analysis
The resulting forward and reverse sequences were assembled using CodonCode Aligner version 4.0.4(CodonCode Corporation) and checked for errors/ambiguities.The partial nucleotide sequences (COI, 16S rRNA) were analyzed with the Basic Local Alignment Search Tool (BLAST) provided by NCBI to compare them with known reference sequences in GenBank.DNA sequences were aligned using MUSCLE (Edgar 2004).The two original mtDNA alignments consisted of 660 bp (COI) and 514 bp (16S rRNA).Yet, 27 bp, including all the gaps, were excluded from the 16S rRNA because of alignment ambiguities, resulting in a fi nal 16S rRNA alignment of 487 bp.The sequences were checked for ambiguous nucleotide sites, saturation and phylogenetic signal using DAMBE version 5.2.65 (Xia 2013).MEGA version 6 (Tamura et al. 2013) was used to (1) perform Tajima's D test for selective neutrality of mutations (D-test statistic), (2) calculate uncorrected pairwise p-distances among sequences, (3) translate COI protein coding sequences into amino acids, (4) check for stop codons, and ( 5) evaluate transition/transversion rates for COI.Substitution models were inferred independently for each of the gene partitions using jModelTest version 2.1.7(Darriba et al. 2012) applying Akaike (1974) weights as selection criterion.
Phylogenetic trees were constructed using maximum likelihood (ML) and Bayesian inference (BI).The shape parameter of the gamma distribution, based on 16 rate categories, was estimated using likelihood analysis.ML trees were inferred with RAxML version 8.0.0 (Stamatakis 2014) through the CIPRES Science Gateway (Miller et al. 2010) using a GTR+G substitution model and 1000 bootstrap replicates to assess branch support.The concatenated sequence alignment was partitioned by gene, allowing different evolution rates for COI and 16S rRNA.
The BI analysis was run for 5 million generations (heating parameter = 0.08), sampling every 1000 generations.Convergence was confi rmed by verifying that the standard deviations of split frequencies were below 0.01.Then the fi rst 2500 trees were discarded as burn-in, so that the fi nal consensus tree was built using the last 7502 trees.The optimal parameters were assessed by calculating the Bayes factor using Tracer version 1.6.0(Rambaut et al. 2014).Support for nodes was defi ned as posterior probabilities.

Species descriptions
Drawings were made using a stereo microscope.We have mainly focused on adult males, although adult females and a few juveniles were also available.For the gonopods we used the terminology of Pimvichai et al. (2009a: 21, fi gs 1-3) with the following abbreviations (terms not used in Pimvichai et al. 2009a, b in bold): ac = anterior coxal fold: the main part of gonopod in anterior view; confusingly called posterior coxal fold by Demange (1961) and Hoffman (1975) aip = additional spine-like process: between lateral and mesal processes of anterior coxal fold alp = lateral process of anterior coxal fold: the distolateral part of the anterior coxal fold amp = mesal process of anterior coxal fold: an additional projection on the anterior coxal fold, protruding from its mesal margin bp = blepharochaete (pl.-ae): the normal form of apical setae, long, slender, stiffened, and usually pigmented, somewhat reminiscent of the mammalian eyelash (Hoffman 1975) cr = longitudinal crest in gutter of palette: a crest which runs along the middle of the gutter near the tip of the palette fe = femoral spine (also fe 1 and fe 2): a usually long, curved spine on the telopodite, originating slightly distal to the point where the telopodite emerges from the coxa lc = longitudinal crest: a strong longitudinal crest at the mesal margin of amp in posterior view ll = lamellar lobe: a small, slightly folded lobe at the basis of the apical part of the telopodite pa = palette: the distalmost lobe of the apical part, carrying the row of blepharochaetae pc = posterior coxal fold: the main part of gonopod in posterior view, usually shorter than ac and forming shelf for accommodation of telopodite shaft plp = lateral process of posterior coxal fold: the lateral part of the posterior coxal fold, normally digitiform pmp = mesal process of posterior coxal fold: the mesal part of the posterior coxal fold, normally forming a shelf for accommodation of telopodite shaft px = paracoxite: the basal, lateral part of the posterior coxal fold sl = spatulate lobe: a distinct distal, separate lobe at the apical part, spatulate, sometimes with a distal spine-like process sls = slender long spine: an additional slender long spine (much longer than ss) at the base of the apical part of telopodite in posterior view ss = small spine: an additional small spine at the base of the apical part of telopodite in posterior view st = sternum: a small, usually triangular sclerite between the basal parts of the anterior coxal folds ti = tibial spine: a usually long spine on the telopodite, originating distal to the femoral spine, at the basis of the apical part of the telopodite, usually curved in the opposite direction of the femoral spine, the two together forming a circle Apical part: the part of the telopodite distal to the tibial spine Shelf: the distal surface of the posterior coxal fold

DNA sequence data
The COI dataset included 660 bp, while the 16S rRNA dataset included 487 bp.The concatenated dataset therefore comprised 1147 bp.
The aligned 16S rRNA gene fragment (487 bp), had nucleotide frequencies of 0.324, 0.100, 0.210 and 0.366 for A, C, G and T, respectively (29.9% GC content).The uncorrected p-distance between the taxa ranged from 0.00 to 0.15 (Appendix I).
The concatenated data set (1147 bp) had nucleotide frequencies of 0.303, 0.172, 0.182 and 0.343 for A, C, G and T, respectively (35.3% GC content).The uncorrected p-distance between the taxa ranged from 0.01 to 0.18.The estimated value of the shape parameter, as evaluated by MEGA version 6 with the substitution pattern and rates fi tted by the GTR+G model (Kumar & Nei 2000), for the discrete Gamma Distributions (Gu et al. 1995) were 0.1542, 0.1584 and 0.1424 for the COI fragment, 16S rRNA (unambiguously aligned sequences) and the concatenated data set, respectively.
D values for selective neutrality were 1.91, 0.63 and 1.45 for COI, 16S and the combined sequences, respectively.The COI transition/transversion ratio was 2.24.
Clade 1A1, Thyropygus allevatus is very strongly supported (100% for ML bootstrap replicates and a BI posterior probability of 1.00).Abbreviations after species names refer to locality names as shown in Table 1.
Clade 1A2, T. cuisinieri subgroup (here represented by foliaceus and jarukchusri) is very strongly supported (100% for ML bootstrap replicates and a BI posterior probability of 1.00).
Clade 1A3 (coloured highlight), is strongly supported by 89% ML bootstrap replicates and strongly supported by a BI posterior probability of 1.00, and is a monophyletic group consisting of a mixture of species from the T. opinatus and T. bifurcus subgroups (Pimvichai et al. 2009a, b), including the nine new species described in this study.Thus, we treat all species of these subgroups as a single subgroup viz., the extended T. opinatus subgroup (henceforward simply "the T. opinatus subgroup").
Gnathochilarium: mentum smooth, at most with a few small setae distally and a large, horseshoe-shaped ridge opening distally.Stipites densely covered with spine-like setae, except for an irregular oblique band from c. middle of lateral margin to border of lamella lingualis; long setae present on the distolateral part; males (not females) distally with a small sclerotised 'island' with 1-3 spine-like setae in middle of an ovoid, poorly sclerotized, hairless area.Lamellae linguales with three long, apical setae and a number of short, basal, spine-like setae.
Ventral pads on postfemora and tibiae on all legs, except fi rst three pairs.
Gonopods: sternum (st) triangular.Anterior coxal fold (ac) basally slender, becoming broader towards tip, lateral margins diverging; distally with two processes: a lateral process (alp) and a usually smaller mesal process (amp); the shape of these processes is species-specifi c. Posterior coxal (pc) fold much lower than anterior coxal fold, basally with moderately high lateral paracoxites (px), distally variously modifi ed but always with a smooth area over which the telopodite shaft can slide.Telopodite with a single or double, well-developed femoral spine (fe) and a long, slender, curved tibial spine (ti).Many species with a characteristic spatulate lobe (sl) originating under the base of the tibial spine; lobe sometimes distally rounded, spoon-like, sometimes ending in a large, stout spine.Other species instead with a small and slightly folded lateral lamella (ll).Apical palette (pa) simple, forming a broad gutter, sometimes with a longitudinal crest in the concavity; apically with a row of 7-13 brownish blepharochaetae (bp).(Demange, 1986)  Mesal margin of lateral process of anterior coxal fold (alp) with a short curved caudad spine; mesal process of anterior coxal fold (amp) as long as alp, straight ...T. demangei Pimvichai et al., 2009b

Etymology
This species is named after the organization "Centre International de Myriapodologie -CIM" (www.myriapodology.org) in recognition of its immense importance for inspiring and supporting research on myriapods.

DNA barcode
The GenBank accession number of the barcode of the holotype is KU306519 (voucher code CUMZ-D00086).Distribution (Fig. 12) Known only from the type locality.

Remarks
Coexisting with the smaller T. forceps sp.nov.

Diagnosis
A species of the opinatus subgroup.Lateral process of anterior coxal fold (alp) fl attened, slightly curved, its laterodistal margin coarsely dentate.Similar in this respect to T. cristagalli, T. implicatus and T. undulatus sp.nov.Differs from these species by having the mesal process of posterior coxal fold (pmp) very high, pointed-triangular, directed almost straightly distad, and by having the tibial spine (ti) recurved.

Etymology
The name is a Latin noun in apposition, meaning "knife", and refers to the knifelike second femoral spine (fe 2).

DNA barcode
The GenBank accession number of the barcode of one of the paratypes is KC519535 (voucher code CUMZ-D00078).Distribution (Fig. 12) Known only from the type locality.

Diagnosis
A species of the opinatus subgroup.Lateral process of anterior coxal fold (alp) long, slender, regularly curved, tip close to tip of opposite alp, the two together forming a circle.Similar in this respect to T. erectus and T. fl oweri.Differs from the former by having a telopodite lobe (lo).Particularly similar to T. fl oweri, differing from it by having the mesal process of posterior coxal fold (pmp) slender, directed distolaterad, laterally with a digitiform process (plp), and by not having pmp strongly developed along the anterior-posterior axis.

Etymology
The name is a Latin noun in apposition, referring to the forceps-like gonopod coxae.

DNA barcode
The GenBank accession number of the barcode of one of the paratypes is KC519531 (voucher code CUMZ-D00073).Distribution (Fig. 12) Known only from Namwang Srithammasokrach and Tham Pha Deang temple in Nakhonsrithammarat Province.

Remarks
Coexisting with the larger T. cimi sp.nov. at Namwang Srithammasokrach.

Diagnosis
A species of the opinatus subgroup.Lateral process of anterior coxal fold (alp) slender, curving mesad; mesal process of anterior coxal fold (amp) as long as alp, slender, straight, directed distad.Similar in these respects to T. demangei.Differs from this species by having mesal margin of amp forming a strong longitudinal crest (lc) in posterior view, by having only one femoral spine (fe) and by having the spatulate lobe (sl) broad and rounded.

Etymology
The name is a Latin adjective, referring to the longitudinal crest on the mesal process of the anterior coxal fold.

DNA barcode
The GenBank accession number of the barcode of the paratype is KC519534 (voucher code CUMZ-D00077).
Distribution (Fig. 12) Known only from the type locality.slender, regularly curved and by having the mesal process (amp) broadly expanded and forming a strong longitudinal crest (lc) in posterior view.

Etymology
The species is named after the Royal Thai Navy, in recognition of their kind assistance which enabled us to pursue the necessary fi eldwork at the type locality, and also after Chulalongkorn University where pink is the symbolic colour of the university, indirectly referring to the pink legs of the species; the species name is treated as a noun in apposition.

DNA barcode
The GenBank accession number of the barcode of one of the paratype is KU306522 (voucher code CUMZ-D00089).
Distribution (Fig. 12) Known only from the type locality.

Diagnosis
A species of the opinatus subgroup.Differs from all other species in the subgroup by having the lateral process of the anterior coxal fold (alp) regularly curved, terminating in a sharp spine pointing slightly distad, and by having the tibial spine (ti) short and fl attened.

DNA barcode
The GenBank accession number of the barcode of the holotype is KU306521 (voucher code CUMZ-D00088).

Diagnosis
A species of the opinatus subgroup.Mesal process of anterior coxal fold (amp) directed slightly obliquely disto-mesad, tip overlapping tip of opposite amp.Similar in this respect to T. chelatus.Differs from this species by having amp as long as alp and by having the lateral process of the posterior coxal fold (plp) as a massive, broad lobe.

Etymology
The species is named in honour of Chirasak Sutcharit in recognition of his devotion to collecting millipedes.

DNA barcode
The GenBank accession number of the barcode of the holotype is KU306524 (voucher code CUMZ-D00090).
Distribution (Fig. 12) Known only from the type locality.

Diagnosis
A species of the opinatus subgroup.Lateral process (alp) fl attened, slightly curved, its laterodistal margin coarsely dentate.Similar in this respect to T. cristagalli, T. implicatus and T. culter sp.nov.Differs from these species by having the mesal process of the posterior coxal fold (pmp) with two rounded distal lobes, visible in anterior view between alp and amp, and by having the tibial spine (ti) not ending in a sharp point.

Etymology
The name is a Latin adjective, referring to the undulate/coarsely dentate laterodistal margin of the lateral process of the anterior coxal fold.

DNA barcode
The GenBank accession number of the barcode of the holotype is KU306520 (voucher code CUMZ-D00087).

Etymology
The name is a Latin noun in apposition, meaning "bear", and refers to the (somewhat) bearhead-like profi le of the lateral process of the anterior coxal fold (alp).

DNA barcode
The GenBank accession number of the barcode of the holotype is KU306523 (voucher code NMHW-Inv.7855).
Distribution (Fig. 12) Known only from the type locality.

Discussion
In our previous treatment of the T. opinatus subgroup (Pimvichai et al. 2009a) we anticipated that DNA sequence data would help resolve relationships between Thyropygus species.A subsequent study (Pimvichai et al. 2014), which included a mtDNA sequence analysis of 23 Thyropygus species, did indeed shed light over these relationships.The present study expands the taxonomic coverage of this latter analysis by adding 9 new species to the mtDNA dataset.The analysis clearly shows that the opinatus and bifurcus subgroups cannot be maintained as separate groups and hence should be merged into an extended opinatus subgroup.The monophyly of this (extended) opinatus subgroup is wellsupported (Fig. 1).However, the phylogenetic relationships within the opinatus subgroup are often not well-resolved.Thus, still more species and nuclear genes need to be added to improve the phylogenetic resolution.
For the delimitation of species in the opinatus subgroup we used both gonopodal characters and divergences in COI sequences (Table 2).Interspecifi c distances among Bavarian Diplopoda ranged between 0% (in subspecies) and 33.18 % (among different orders).The mean value of the interspecifi c distance for Diplopoda was 14.17% (Spelda et al. 2011).For 21 species of the opinatus subgroup the divergence was between 2 and 17%.The lowest interspecifi c divergence was 2%, between T. cristagalli and T. planispina sp.nov., between T. planispina sp.nov.and T. undulatus sp.nov., and between T. quadricuspis and T. ursus sp.nov.In these cases, speciation may have been very recent but we still need more DNA data and more samples to be able to infer the taxonomic status and relationships among these species properly.However, all these species show unique gonopodal characters and occur in separate localities and hence, for the time being, we treat these "morphs" as separate species.Nevertheless, intraspecifi c divergence of T. induratus was also 2% in the samples of the present study but despite this sequence divergence the overall gonopod confi guration remained the same in the specimens here assigned to T. induratus.
All species of the opinatus subgroup share a conspicuous synapomorphy, viz. the additional projection on the anterior coxal fold.This character separates the opinatus subgroup from other subgroups of the T. allevatus group.
The main gonopodal variations within the opinatus subgroup are: -the number of femoral spines (fe): one or two -presence/absence of a spatulate lobe (sl) -presence/absence of a longitudinal rounded crest (cr) near the tip of the palette -presence/absence of a small spine (ss) or a slender long spine (sls) at the posterior base of the apical part of the telopodite -presence/absence of a lateral (plp) and a mesal (pmp) process of the posterior coxal fold The vast majority of the species we have studied and of the previous records of the T. opinatus subgroup are from the southern part of Thailand (Fig. 12).The remaining material comes from the North and West of Thailand and from some adjacent areas in Myanmar and Malaysia.Most probably further fi eldwork in all these areas will lead to the discovery of further new species of Thyropygus.

Fig. 1 .
Fig. 1.Phylogenetic relationships of Thyropygus species based on maximum likelihood analysis (ML) and Bayesian Inference (BI) of 1147 bp of concatenated gene fragments of COI (660 bp) and 16S rRNA (487 bp).Numbers at nodes indicate branch support based on bootstrapping (ML) / posterior probability (BI).Scale bar = 0.06 substitutions/site.# indicates branches which received < 50% ML bootstrap support, -indicates non-supported branches by posterior probability.Clade memberships and designations are shown as vertical bars; 1A1 = T. allevatus, 1A2 = cuisinieri subgroup, 1A3 = opinatus subgroup and 1A4 = induratus subgroup.The coloured area marks the T. opinatus subgroup.Abbreviations after species names refer to locality names as shown in Table1.

Fig. 12 .
Fig. 12. Known distribution of the species of the T. opinatus subgroup.

Table 2 .
Estimates of evolutionary divergence between COI sequences expressed as p-distance.