Gonopods were photographed with a digital microscope camera (Zeiss Stemi 305). Samples for scanning electron microscopy (SEM: Hitachi TM4000Plus) were air-dried directly from alcohol and sputter-coated for 60 s with gold (Hitachi: MC1000). Scanning electron micrographs were taken at the Central Lab of Mahasarakham University. Drawings were made using a stereomicroscope and photographs. Voucher specimens were deposited in the collections of the Museum of Zoology, Chulalongkorn University, Bangkok, Thailand (CUMZ).

Total genomic DNA was extracted from legs of specimens of Coxobolellus saratani sp. nov. (CUMZ-D00153 and CUMZ-D00153) from Loei Province and Coxobolellus serratoligulatus sp. nov. (CUMZ-D00154-1 and CUMZ-D00154) from Uttaradit Province, Thailand using the PureDireX column based genomic DNA extraction kit (tissue) (Bio-Helix) following the manufacturer’s instructions. PCR amplifications and sequencing of the standard mitochondrial COI and 16S rRNA gene fragments were done as described by Pimvichai et al. (2020). The COI and 16S rRNA gene fragments were amplified with the primers LCO-1490 and HCO-2198 (Folmer et al. 1994) for COI and 16Sar and 16Sbr (Kessing et al. 2004) for 16S rRNA. All new nucleotide sequences have been deposited in GenBank under accession numbers OP580097–OP580100 and OP580512–OP580515 for the partial COI and 16S rRNA fragment sequences respectively. Sample data and voucher codes are provided in Table 1.

The 16S rRNA data included 51 specimens, representing 12 genera and 37 nominal species of ingroup taxa (Table 1). Three species of the order Spirostreptida, viz. Anurostreptus barthelemyae Demange, 1961 (Harpagophoridae), Chonecambala crassicauda Mauriès & Enghoff, 1990 (Pericambalidae) and Thyropygus allevatus (Karsch, 1881) (Harpagophoridae) were used as the outgroup. The same ingroup and outgroup taxa were used for COI, with the addition of: (1) Coxobolellus simplex Pimvichai, Enghoff, Panha & Backeljau, 2020, (2) C. albiceps (Stpl) Pimvichai, Enghoff, Panha & Backeljau, 2020, (3) C. valvatus (TCD) Pimvichai, Enghoff, Panha & Backeljau, 2020, (4) Siliquobolellus amicusdraconis Pimvichai, Enghoff, Panha & Backeljau, 2022, (5) S. constrictus Pimvichai, Enghoff, Panha & Backeljau, 2022, (6) S. prasankokae Pimvichai, Enghoff, Panha & Backeljau, 2022, (7) Paraspirobolus lucifugus (Gervais, 1836), (8) Narceus annularis Rafinesque, 1820, (9) Trachelomegalus sp., (10) Macrurobolus macrurus (Pocock, 1893), (11) Atopochetus anaticeps Pimvichai, Enghoff, Panha & Backeljau, 2018, and (12) Benoitolus birgitae (Hoffman, 1981). Hence, the COI data set included 63 specimens, representing 18 genera and 47 nominal species of ingroup taxa (Table 1).

CodonCode Aligner (ver. 4.0.4, CodonCode Corporation) was used to assemble the forward and reverse sequences and to check for errors and ambiguities. All sequences were checked with the Basic Local Alignment Search Tool (BLAST) provided by NCBI and compared with reference sequences in GenBank. They were aligned using MUSCLE (ver. 3.6, see http://www.drive5.com/ muscle; Edgar 2004). The COI alignment consisted of 660 bp and that of 16S rRNA consisted of 458 bp (gaps were excluded by complete-deletion). The sequences were checked for ambiguous nucleotide sites, saturation and phylogenetic signal using DAMBE (ver. 5.2.65; see http://www.dambe.bio.uottawa.ca/DAMBE/dambe.aspx; Xia 2018). MEGA (ver. X, see http://www.megasoftware.net; Kumar et al. 2018) was used to (1) check for stop codons, (2) translate COI protein-coding sequences into amino acids, and (3) calculate uncorrected pairwise p-distances among sequences.

Phylogenetic trees were constructed using maximum likelihood (ML) and Bayesian inference (BI) applied to either the COI and 16S rRNA data separately or combined. The shape parameter of the gamma distribution, based on 16 rate categories, was estimated using maximum likelihood analysis. ML trees were inferred with RAxML (ver. 8.2.12, see http://www.phylo.org/index.php/tools/raxmlhpc2_tgb.html; 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. BI trees were constructed with MrBayes (ver. 3.2.7a, see http://www.phylo.org/index.php/tools/mrbayes_xsede.html; Huelsenbeck and Ronquist 2001). Substitution models were inferred using jModeltest (ver. 2.1.10, see https://www.github.com/ddarriba/jmodeltest2/releases; Darriba et al. 2012) applying Akaike Information Criterion weights as selection criterion. This yielded as best models GTR+I+G (–lnL = 17117.3442, gamma shape = 0.6820) for the combined dataset, TrN+ I+G (–lnL = 12352.6533, gamma shape = 0.4820) for COI, and TIM3+ I+G (–lnL = 6495.3310, gamma shape = 0.8870) for 16S rRNA.

BI trees were run for 2 million generations for the combined dataset, 3 million generations for the separate 16S rRNA, and 4 million generations for the separate COI datasets (heating parameter: 0.01 for the combined dataset and 16S rRNA, and 0.02 for COI), sampling every 1000 generations. Convergences were confirmed by verifying that the standard deviations of split frequencies were below 0.01. Then the first 1000 trees were discarded as burn-in, so that the final consensus tree was built from the last 3002 trees for the combined dataset, 4502 trees for the 16S rRNA, and 6002 trees for the COI datasets. Branch support was assessed by posterior probabilities.

For ML trees we consider branches with bootstrap values (BV) of ≥ 70% to be well supported (Hillis and Bull 1993) and < 70% as poorly supported. For BI trees, we consider branches with posterior probabilities (PP) of ≥ 0.95 to be well supported (San Mauro and Agorreta 2010) and below as poorly supported.

The uncorrected p-distance between the COI sequences (660 bp) ranged from 0.00 to 0.26 (Suppl. material 1). The mean interspecific sequence divergence within Coxobolellus was 0.12 (range: 0.06–0.15). The mean intergeneric sequence divergence between Coxobolellus and Pseudospirobolellus was 0.21 (range: 0.20–0.23), between Coxobolellus and Siliquobolellus it was 0.17 (range: 0.14–0.20), and between Coxobolellus and Benoitolus birgitae it was 0.21 (range: 0.20–0.23).

The uncorrected p-distance between the 16S rRNA (458 bp) sequences ranged from 0.00 to 0.30 (Suppl. material 2). The mean interspecific sequence divergence within Coxobolellus was 0.11 (range: 0.05–0.15). The mean intergeneric sequence divergence between Coxobolellus and Pseudospirobolellus was 0.21 (range: 0.18–0.24). Currently, there are no 16S rRNA sequences for Siliquobolellus and Benoitolus.

The uncorrected p-distance between the combined COI + 16S rRNA (1118 bp) sequences ranged from 0.00 to 0.27 (Suppl. material 3). The mean interspecific sequence divergence within Coxobolellus was 0.11 (range: 0.05–0.14). The mean intergeneric sequence divergence between Coxobolellus and Pseudospirobolellus was 0.21 (range: 0.19–0.23).

The ML and BI trees based on the separate and combined datasets (COI, 16S rRNA and COI + 16S rRNA) were largely congruent with respect to the well-supported branches (by visual inspection of the branching pattern). The combined COI + 16S rRNA tree is used for further discussion (Fig. 1). The separate COI and 16S rRNA trees are presented in Suppl. materials 4, 5, respectively.

Figure 1. 

Phylogenetic relationships of Coxobolellus species based on maximum likelihood analysis (ML) and Bayesian inference (BI) of 1118 bp in the combined COI + 16S rRNA alignment. Numbers at nodes indicate branch support based on bootstrapping (ML) / posterior probabilities (BI). Scale bar = 0.07 substitutions/site. # indicates branches with <50% ML bootstrap support and <0.95 Bayesian posterior probability. - indicates non-supported branches. The colored areas mark the subclades of Coxobolellus and are labelled as in Pimvichai et al. (2020).

The genus Coxobolellus (Clade 1) is consistently well supported in the three trees and can be subdivided into five generally well supported subclades (1A–E; cf. Pimvichai et. al. 2020), one of which is formed by the two new species (subclade 1E). Yet, this latter subclade is only well supported by ML of the 16S rRNA (and combined) data, but not by BI and the COI data. Moreover, the position of two other species, viz. C. fuscus (COI and 16S rRNA) and C. simplex (COI only), remains ambiguous, even if both species seem to be somehow associated with subclades 1C and 1D.

Subclade 1A comprises C. nodosus and C. valvatus, two species that are distributed in western and northern Thailand, respectively. They differ in their posterior gonopod telopodite: in C. nodosus the telopodital part (pt) is as long as the coxal part (pcx), while in C. valvatus the telopodital part (pt) is much shorter than the coxal part (pcx). Their COI sequence divergence is 0.06–0.07.

Subclade 1B comprises the sympatric species C. albiceps and C. transversalis. They mainly differ in the following two aspects: (1) in C. albiceps the tip of the anterior gonopod coxa is apically obliquely truncated, while in C. transversalis it is transversely truncated; (2) the telopodital part (pt) of the posterior gonopod telopodite is fairly long in C. transversalis, while in C. albiceps it is short. Their COI sequence divergence is 0.08.

Subclade 1C comprises C. compactogonus and C. tenebris, two species that are distributed in northeastern and central to western Thailand, respectively. They mainly differ in the following two aspects: (1) the anterior gonopod telopodite (at) ends in a rounded lobe in C. compactogonus, while in C. tenebris there is a tiny triangular laterad denticle near the tip of the anterior gonpod telopodite (at); (2) the telopodital part (pt) of the posterior gonopod telopodite is with three processes in C. compactogonus, while in C. tenebris this part ends in a rounded, smooth lobe. Their COI sequence divergence is 0.10 (range: 0.09–0.10).

Subclade 1D comprises C. serratus and C. tigris, two species that are distributed in southern Thailand. They mainly differ in the following two aspects: (1) in C. serratus the posterior gonopod telopodite with a long coxal part (pcx), while this part is short in C. tigris; and (2) the telopodital part (pt) laterally is with a serrate margin in C. serratus, while in C. tigris only the apical part with a serrate margin. Their COI sequence divergence is 0.12.

Subclade 1E (only supported by ML of 16S rRNA) comprises the two new species: C. saratani sp. nov. and C. serratoligulatus sp. nov., which are distributed in northeastern and northern Thailand, respectively. They differ in both the anterior gonopod and the posterior gonopod telopodite, and will be treated in detail further below. Their COI sequence divergence is 0.08 (range: 0.07–0.08).

The combined COI + 16S rRNA (Fig. 1) and the separate 16S rRNA (Suppl. material 5) trees showed Pseudospirobolellus as a well-supported sister group of Coxobolellus. However, this sister group relation was no longer supported when the genus Siliquobolellus was included, i.e. in the separate COI tree (Suppl. material 4).