Sample design in biodiversity studies matters: a fine-scale study of Lawrence’s velvet worm, Peripatopsis lawrencei (Onychophora: Peripatopsidae), reveals hidden diversity

Sample collection

A total of 69 new Peripatopsis lawrencei sensu lato (s.l.). specimens was collected from the known distribution range of the species in the Western Cape province, South Africa (Fig. 1). These specimens were combined with data collected from two previous studies (Daniels et al. 2009; McDonald and Daniels 2012), yielding a total of 119 specimens (Table 1). Peripatopsis lawrencei s.s. comprised 116 of the samples, whereas 3 samples were suspected to be the unknown sympatric lineage at Oubos (Riviersonderend Mountains). A handheld GPS was used to record the locality information and a permit (CN44-87-22079) was obtained from CapeNature for velvet worm sampling. Specimens were collected in Afromontane forest patches under decaying logs, leaf litter and stones (McDonald et al. 2012). These specimens were preserved directly in absolute ethanol and stored in a refrigerator at 4°C once returned to the laboratory.

Fig. 1.

Map showing all sampling localities for Peripatopsis lawrencei s.l., Western Cape Province, South Africa. The red markers indicate sampled localities with numbers corresponding to Table 1.

DNA sequencing

A tissue biopsy was taken from each specimen (2–3 mm) and a Macherey–Nagel DNA extraction kit was used, following the manufacturer’s protocol. DNA extractions were kept at 4°C until required for use in polymerase chain reaction (PCR). Two loci were targeted for amplification, namely mitochondrial (mtDNA) cytochrome c oxidase subunit I (COI) and nuclear 18S rRNA subunit (18S). These two loci have been extensively used in velvet worm systematics studies of Peripatopsis, allowing us to combine our data with DNA sequences generated during earlier studies (Daniels et al. 2009, 2013, 2016; McDonald and Daniels 2012; Barnes et al. 2020; Grobler et al. 2023). The COI and 18S rRNA primer pairs (specifically 5F and 7R for the latter locus) used in the amplification of the two fragments were obtained from Folmer et al. (1994) and Giribet et al. (1996) respectively. All newly collected specimens were sequenced for COI. Only a single representative of each locality was sequenced for the 18S locus, except for the Oubos locality at which the two lineages (A and B) occurred in sympatry. A total of 14 new 18S rRNA sequences was generated during our study and combined with three 18S rRNA sequences generated by McDonald and Daniels (2012).

PCRs were conducted for a 25-μL reaction and the samples held at 15°C afterwards. The reactions both contained 14.9 μL of deionised water, 2.5 μL of 10× Mg²⁺ free buffer, 0.5 μL of a 10-mM dNTP solution and 0.5 μL of the primer sets at 10 mM, 0.1 units of JMR SuperTherm Taq polymerase and 1 μL of extracted DNA. The PCR was conducted on a geneAmp PCR System Thermocycler. PCR conditions for the COI locus were as follows: 94°C for 4 min, 94°C for 30 s, 42°C for 35 s and 72°C for 40 s for 36 cycles, with an extension of 72°C for 10 min. For the 18S locus, PCR conditions were as follows: 94°C for 4 min, 94°C for 30 s, 48°C for 40 s and 72°C for 40 s for 34 cycles, with an extension of 72°C for 10 min. Following successful PCR amplification, PCR products were electrophoresed in a 1% agarose gel for 180 min, at 100 V and gel-purified using a Bio flux gel purification kit, following the manufacturer's instruction. An ABI 3730 machine was used to sequence the products at the DNA sequence facility of the University of Stellenbosch.

Phylogenetic analyses

Data collected from the sequences were utilised with the DNA sequences for P. lawrencei s.l. from Daniels et al. (2009) and McDonald and Daniels (2012). Forward and reverse strands were used to compute a consensus sequence and check for base pair ambiguity, whereas sequence alignment for both the COI and 18S rRNA data was computed using ClustalX (ver. 2.0, see https://www.genome.jp/tools-bin/clustalw;Thompson et al. 1997). For 18S rRNA, highly variable indel regions that could not be aligned with accuracy were excluded from the phylogenetic analyses using Gblocks (ver. 0.91b, see http://phylogeny.lirmm.fr/phylo_cgi/one_task.cgi?task_type=gblocks; Talavera and Castresana 2007). Using jModelTest2 (ver. 2.2.10, see https://github.com/ddarriba/jmodeltest2; Posada 2008), the Akaike information criterion (AIC) (Akaike 1973) was utilised to determine the best-fit DNA substitution model on XSEDE (ver. 2.0, see https://www.xsede.org/) in CIPRES (see http://www.phylo.org/; Miller et al. 2010). Bayesian inference (BI) and maximum likelihood (ML) approaches were used to construct phylogenies. Analyses were conducted on CIPRES Science Gateway (Miller et al. 2010), with the Bayesian analysis taking place in MrBayes (ver. 3.2.6, see https://github.com/NBISweden/MrBayes/; Ronquist et al. 2012) and the phylogeny for ML in RAxML (ver. 7.2.8 alpha, see https://github.com/stamatak/standard-RAxML; Stamatakis et al. 2008). For the part of the study focusing on P. lawrencei s.s., the ML analyses were conducted on the IQ-TREE web server (ver. 2.2.2.6, see http://iqtree.cibiv.univie.ac.at/; Trifinopoulos et al. 2016) and BI analyses through CIPRES, to create the tree topology for the COI data set. For the phylogenetic placement of the new lineage, a combined tree topology for COI and 18S rRNA was reconstructed using the aforementioned procedures, and DNA sequences for both loci for all described Peripatopsis species were included to determine the phylogenetic placement of the Riviersonderend specimens (Daniels et al. 2009, 2013; Ruhberg and Daniels 2013; Barnes et al. 2020; Grobler et al. 2023). Uncorrected COI and 18S rRNA p-distances were calculated in PAUP* (ver. 4.10, see http://phylosolutions.com/paup-test/; Swofford 2002) and compared to interspecific sequence values recorded among velvet worms.

Branch support values were calculated through an ultrafast Bootstrap analysis with 10 000 pseudo replicates and bootstrap (BS) values >75% were taken as support for the nodes using IQ-TREE. For Bayesian inference, burn-in was set to 20% and a posterior probability (PP) >0.95 was considered statistically well supported, estimated through percentage of time spent on node recovery (Barnes and Daniels 2022). With a Bayesian framework, a divergence time estimation was conducted for the COI and 18S rRNA loci, using a probability model to describe molecular sequences of divergent lineages. The ages of the clades were subsequently estimated using the Markov Chain Monte Carlo that was run for 50 million generations, sampling the chain every 1000th generation. In total, 20% of generations were discarded as burn-in. Convergence was assessed visually by examining the performance of the chain, as in Zhang et al. (2013). A relaxed molecular clock was implemented, running the analysis through BEAST (ver. 2.7.6, see https://www.beast2.org/; Drummond and Rambaut 2007) and a multiple coalescent model was used (Heled and Drummond 2010). Mutation rates of 1.5–2.3% per million years were used for the COI locus and the mutation rate for the 18S rRNA locus was estimated through the use of a non-informative (1/x) prior. The aforementioned use of mutation rates has been applied in several Peripatopsis divergence time estimations (Barnes et al. 2020; Myburgh and Daniels 2022; Grobler et al. 2023). Finally, for further support we conducted a Shimodaira–Hasegawa test (Shimodaira and Hasegawa 1999) to evaluate the likelihood of a constrained tree topology, where clades A and B are combined into a single clade, rather than being represented as two separate entities.

Phylogeographic analysis of P. lawrencei s.s.

A haplotype network was constructed through the use of TCS (ver. 1.21, see https://bioresearch.byu.edu/tcs/), with a 95% confidence interval (Clement et al. 2000) and population genetic structure was computed in ARLEQUIN (ver. 3.5.2.2, see http://cmpg.unibe.ch/software/arlequin35; Excoffier and Lischer 2010) using the COI data for P. lawrencei s.s. The sample localities were subjected to a hierarchical Analysis of Molecular Variance (AMOVA), and F_ST values were calculated, along with standard diversity indices including nucleotide diversity, number of haplotypes and haplotype diversity at each locality.

Morphology and scanning electron microscopy (SEM)

A Canon EOS camera was used to photograph the specimens, and gross morphological traits were examined using a Leica MZ 7.5 stereomicroscope, noting dorsal and ventral integument colour, number of leg pairs and the presence of a claw on the posterior leg pair. A digital calliper was used to measure total length of the ethanol-preserved specimen from the anterior-most point to the posterior side, and additional comparisons between P. lawrencei s.s. and the new lineage were undertaken using a Leica S9i digital stereo microscope at the South African Museum of Cape Town. The rarity of the novel lineage (clade B) poses limitations on the methods of examination. Only two individuals of the new lineage were available for examination and due to the destructive nature of the process, traditional scanning electron microscopy (SEM) was not possible. However, the South African Museum of Cape Town has an environmental SEM, Hitachi TM4000Plus that allows for the analysis of specimens without requiring critical point drying or sputter coating. This allowed us to preserve our limited specimens for later use as holotype and paratype. Through the use of this SEM technique, dorsal and ventral scale ranks for representatives of each of the two main clades were examined. Specimens were accessioned into the entomology collection of the South African Museum of Cape Town (SAM-ENW-C013473–SAM-ENW-C013487). The new species name was deposited in ZooBank (urn:lsid:zoobank.org:act:71276915-79C6-4759-9244-20BEF6303F58).

Phylogenetic analysis of COI for P. lawrencei s.l.

We amplified and sequenced a 644 base-pair (bp) fragment of the COI locus from 69 P. lawrencei s.l. specimens, two of which were the suspected novel lineage. Novel sequences were deposited in GenBank (Accession numbers OR423143–OR423211). We combined COI data generated during our study with 50 P. lawrencei s.l. sequences, of which 1 is known to be the novel species, from two previous studies to yield a total of 119 sequences, 116 of P. lawrencei s.s. and 3 of the novel lineage (clade B). Using corrected AIC, we selected TIM1 + G as our DNA substitution model (–lnL = 2224.64). The base frequencies for the COI locus were as follows: A = 25.00%, C = 12.56%, G = 16.29% and T = 46.14%, whereas the rate matrix included R(a) [A–C] = 1.00, R(b) [A–G] = 19.17, R(c) [A–T] = 2.72, R(d) [C–G] = 2.72, R(e) [C–T] = 8.76 and R(f) [G–T] = 1.00. Gamma (G) shape distribution was 0.15. Both Bayesian inference and maximum likelihood analyses yielded near-identical tree topologies, hence only the Bayesian inference topology is shown (Fig. 2). The tree topology shows two highly divergent statistically well-supported monophyletic clades (A & B) (>75% bootstrap support, BS/>0.95 posterior probability, PP). Within P. lawrencei s.s. (clade A) (84% BS/<0.95 PP), the most basal lineage occurs at Rondevlei NR and Jonkershoek NR (77% BS/0.95 PP). A clade comprising Villiersdorp, Jonaskop, Greyton, Riviersonderend and Oubos A (haplocluster A2; Fig. 2, 3) is sister group to a clade found exclusively at Caledon (92% BS/<0.95 PP), likely reflecting a recent colonisation history between the Riviersonderend Mountains and Klein Swartberg at Caledon (Fig. 2). A coastal clade exists at Vogelgat and Fernkloof NRs (haplocluster A1; Fig. 2, 3) that is sister group to the groups found further south-east at Grootbos, Linden Hill, Sandies Glen and Napier (haplocluster A3; Fig. 2, 3). These groups are found to be sister group to a clade comprising Grabouw, Dappat se Gat, Franschhoek, High Noon and Kogelberg Biosphere NR (haplocluster A3; Fig. 2, 3). Clade B, the novel lineage, was exclusive to the Riviersonderend Mountains and sympatric with P. lawrencei s.s. at Oubos B.

Fig. 2.

Bayesian phylogenetic tree topology derived from the COI locus, showing the evolutionary relationship across the P. lawrencei s.s. population and separation from the sympatric novel lineage. Statistical support for the nodes is shown as posterior probability (>0.95 PP) above the nodes and bootstrap values (>75% BS) below the node. Posterior probability and bootstrap values <0.95 PP and <75% BS respectively are indicated with a hash (‘#’). Clades A and B have been indicated, showing Peripatopsis lawrencei s.s. and the novel lineage respectively. Clade A (P. lawrencei s.s.) is divided into three haploclusters (A1–A3) and corresponds to the haplotype network in Fig. 3.

Fig. 3.

The haplotype network based on the COI sequence data showing the genetic distribution of Peripatopsis lawrencei s.s. across 19 localities in the Western Cape, South Africa. A total of 57 haplotypes was identified, with circle size corresponding to the frequency of the haplotypes. The codes A1–A3 represent the clades (haploclusters) present on Fig. 2. Missing or unsampled haplotypes are indicated with a solid black circle. Haplotype numbers correspond to those in Appendix A1.

The uncorrected p-distance for the COI locus between P. lawrencei s.s. (clade A) and clade B ranged from 9.13 to 10.62%. The uncorrected p-distance between P. capensis and P. overbergiensis was 6.58% (our two outgroups), and between P. overbergiensis and P. lawrencei s.s., this ranged from 5.93 to 7.49%. The uncorrected p-distance ranged from 10.53 to 10.73% between clade B and P. overbergiensis. The marked p-distance between the two clades (A and B) of P. lawrencei s.l. further corroborates the genetic distinction of clade B.

Phylogeographic analysis of P. lawrencei s.s. (clade A)

The TCS analyses of 116 COI sequences revealed 57 haplotypes, comprising 3 large unconnected haploclusters that were genetically and geographically discrete (Fig. 3). Haplocluster A1 was exclusive to a narrow coastal strip at Hermanus (Fernkloof and Vogelgat Nature Reserves), whereas haplocluster A2 was restricted to the Riviersonderend Mountains (Oubos A, Greyton, Villiersdorp, Jonaskop and Riviersonderend). Finally, haplocluster A3 was geographically more widespread from Rondevlei to Napier. Within the aforementioned three haploclusters, we observed evidence for maternal dispersal exclusively between Fernkloof NR and Vogelgat, separated by 18 km because these two localities shared one haplotype (Appendix A1). Overall no evidence of shared haplotypes existed among the remaining sample localities, indicating a lack of maternal dispersal based on the COI sequence data (Appendix A1). The lack of maternal dispersal is validated by the AMOVA results. The AMOVA demonstrated that 83.06% of the genetic variation was present among populations (Va = 6.32, d.f. = 18, s.s. = 693.91, P < 0.01), whereas 16.94% of the genetic variation was present within populations (Vb = 1.29, d.f. = 97, s.s. = 125.05, P < 0.01). Furthermore, pairwise and statistically significant F_ST values showed marked to moderate levels of genetic differentiation among conspecific populations (Table 2, Fig. 4). Non-significant results in Sandies Glen, Riviersonderend, Jonaskop, Villiersdorp and Rondevlei can likely be attributed to small sample sizes (n = 1) at these localities.

Fig. 4.

Heatmap diagram with the colour showing the strength of the F_ST value for each respective pairwise F_ST combination. Darker colours show higher levels of F_ST between different localities and are therefore more distinct. The stronger the signal, the more easily the different populations can be distinguished from one another by looking at the individual's DNA. The lighter the colour, the less distinct. Blocks marked by a superimposed X indicate non-significance.

The number of samples, haplotypes, polymorphic sites and amounts of haplotype (h) and nucleotide diversity (π) are reported in Table 3 for all 19 localities. The number of haplotypes ranged from one to seven, with most occurring between one and three. The highest haplotype and nucleotide diversities were found in Greyton (h = 1.0000 ± 0.2722; π = 0.009081 ± 0.007451) despite the comparatively low sample size (n = 3) and the second highest was found in Linden Hill (h = 0.9000 ± 0.1610; π = 0.008696 ± 0.005872). The lowest haplotype diversity (with n_h > 1) was found in Grabouw (h = 0.4643 ± 0.2000) and the lowest nucleotide diversity (with n_h > 1) in Franschhoek (π = 0.000776 ± 0.000963). Rondevlei NR, Riviersonderend, Sandies Glen, Villiersdorp and Jonaskop all had only one haplotype. The highest diversity was present in Greyton and Linden Hill.

Combined DNA sequence (COI and 18S rRNA) phylogeny and divergence time estimation

We amplified a 519-bp fragment of the 18S rRNA locus for 14 P. lawrencei s.l. specimens, one of which was the suspected novel lineage (clade B). The novel 18S rRNA data were deposited in GenBank (Accession numbers: OR421492–OR421505). The uncorrected p-distance between the two sympatric Oubos lineages (A and B) was 5.58%. To determine the phylogenetic placement of the novel lineage (clade B), we combined all 18S rRNA sequences with COI sequence data from published sequence data for Peripatopsis.

The DNA substitution model for 18S rRNA, selected using AIC, was TIM1 + I + G (–lnL = 1242.71). The base frequencies for the 18S rRNA locus were as follows: A = 19.37%, C = 27.27%, G = 30.28% and T = 23.09%, whereas the rate matrix included R(a) [A–C] = 1.00, R(b) [A–G] = 0.94, R(c) [A–T] = 2.08, R(d) [C–G] = 2.08, R(e) [C–T] = 4.02 and R(f) [G–T] = 1.00. The gamma (G) shape distribution was 1.00. For COI of Peripatopsis, the selected DNA substitution model was GTR + I + G (–lnL = 6422.14) and base frequencies were A = 34.69%, C = 3.71%, G = 9.53% and T = 52.07%, with the rate matrix being R(a) [A–C] = 8.47, R(b) [A–G] = 52.93, R(c) [A–T] = 2.01, R(d) [C–G] = 24.23, R(e) [C–T] = 75.49 and R(f) [G–T] = 1.00. The gamma (G) shape distribution was 0.35. The combined COI and 18S rRNA data yielded a 1163-bp fragment, and the topologies produced from maximum likelihood and Bayesian inference were nearly identical, therefore only the maximum likelihood tree topology is shown (Fig. 5).

Fig. 5.

Divergence time estimation based on COI and 18S rRNA sequence data for Peripatopsis. The 95% HPD interval bars are shown at each node, with a timescale bar representing the relevant epochs. Clade A shows the monophyletic P. lawrencei s.s. whereas clade B indicates the ancestral novel lineage, denoted as P. aereus sp. nov. Blue bars indicate HPD ranges. Posterior probability values >0.95 PP are shown above each node and bootstrap values >75% BS are shown below the node. Posterior probability and bootstrap values <0.95 PP and <75% BS respectively are indicated with a hash (‘#’).

Peripatopsis was retrieved as monophyletic (>75% BS/>0.95 PP) (Fig. 5). The most basal lineage was P. alba, found in the Wynberg Caves on Table Mountain. The latter species was sister group to a clade of P. balfouri exclusively from the Kogelberg and Cape Peninsula. Peripatopsis balfouri was in turn sister group to a clade where P. bolandi Daniels, McDonald & Picker, 2013 was sister to P. purpureus Daniels, McDonald & Picker, 2013. The next clade comprised P. ferrox Barnes, Reiss & Daniels, 2020 sister to P. edenensis Barnes, Reiss & Daniels, 2020 and these two species were sister group to a clade comprising P. tulbaghensis Barnes, Reiss & Daniels, 2020, sister to P. clavigera Purcell, 1899 and the latter species was sister to P. cederbergiensis Daniels, McDonald & Picker, 2013. The novel lineage, P. aereus sp. nov. (clade B) was basal to the remaining Peripatopsis species, affirming the evolutionary distinction from the sympatric P. lawrencei s.s. specimens at Oubos A. In the next clade, we observed a monophyletic P. lawrencei s.s. (clade A) as sister group to P. capensis. The latter species is in turn sister group to P. overbergiensis. The latter clade, containing the three aforementioned species P. lawrencei s.s., P. capensis and P. overbergiensis was sister group to a larger clade of species distributed from the southern Cape into the Eastern Cape and KwaZulu–Natal. Within the last clade, P. sedgwicki Purcell, 1899 was basal to two other clades. Within these, the first clade comprised P. janni Ruhberg & Daniels, 2013 as sister group to three species P. birgeri Ruhberg & Daniels, 2013, sister group to P. polychroma with these two species being sister group to P. hamerae Ruhberg & Daniels, 2013. The second sister clade comprised P. storchi Ruhberg & Daniels, 2013 as sister group to P. moseleyi. Results of the Shimodaira-Hasegawa test indicated that the difference between the unconstrained topologies and constrained topology enforcing the monophyly of the Oubos clades (A sister to B) is highly significant (Δ –lnL = 185.11, P < 0.01), further validating the distinction of the two clades.

Divergence time estimates (Fig. 5) based on the combined COI and 18S data set suggested that Peripatopsis originated 28.83 Ma (95% highest posterior density, HPD 21.27–37.49 Ma) during the Oligocene–Eocene. Peripatopsis alba originated 13.62 Ma (95% HPD: 10.37–17.42 Ma), whereas P. balfouri separated shortly afterwards, 13.38 Ma (95% HPD: 9.53–15.73 Ma). Sister clades P. purpureus and P. bolandi diverged 11.21 Ma (95% HPD: 8.62–14.44 Ma), whereas P. ferox and P. edenensis diverged 10.23 Ma (95% HPD: 6.62–14.35 Ma). Peripatopsis tulbaghensis originated 11.07 Ma (95% HPD: 7.90–14.70 Ma) and P. clavigera diverged from P. cederbergiensis 8.72 Ma (95% HPD: 5.94–11.88 Ma). Furthermore, the novel lineage (clade B) diverged from the remaining Peripatopsis lineages 13.68 Ma (95% HPD: 10.46–17.73 Ma) during the Miocene. Within P. lawrencei s.s., sister clades P. overbergiensis and P. capensis diverged 6.02 Ma (95% HPD: 4.34–8.11 Ma), whereas P. lawrencei s.s. originated 5.23 Ma (95% HPD: 3.75–7.02 Ma) during the late Miocene to the early Pliocene. In a similar time period, within the more derived members of Peripatopsis, P. sedgwicki diverged 8.69 Ma (95% HPD: 6.59–11.11 Ma), whereas P. moseleyi diverged 5.17 Ma (95% HPD: 3.79–6.82 Ma). During the Pliocene, we observed the origin of P. janni 5.26 Ma (95% HPD: 3.97–6.74 Ma), P. storchi 3.77 Ma (95% HPD: 2.59–5.18 Ma), P. hamerae 4.60 Ma (95% HPD: 3.35–6.96 Ma) and the split between P. polychroma and P. birgeri 3.41 Ma (95% HPD: 2.32–4.81 Ma).

Morphology of P. lawrencei s.l.

Diagnostic gross morphological differences were observed between clades A and B, including the number of leg pairs and ventral colour (Table 4). In both clades, dorsal colour ranged from navy blue to slate black, to black with tints of orange and various shades of orange. The ventral colour for P. lawrencei s.s. was predominantly cream white and sometimes pale orange or yellow (Fig. 6a, c). For the new lineage (clade B) the ventral colour was a bronze-like golden brown (Fig. 6b, d). The presence of two claws was noted on the back legs of all individuals across both clades. The number of leg pairs differed between the two clades, as P. lawrencei s.s. had 17 leg pairs and clade B had 18, with the last leg pair being reduced in size. The SEM micrographs revealed marked differences in the number of scale rank counts on the primary papillae between the two lineages, with clade B having more than double the number of scale rank counts for both the dorsal and ventral papillae. P. lawrencei s.s. had seven scale rank counts on the dorsal papilla and four on the ventral papilla (Fig. 7a, c), whereas clade B had fifteen scale ranks for the dorsal papilla and nine on the ventral papilla (Fig. 7b, d).

Fig. 6.

Photographs of recently collected velvet worm specimens suspended in ethanol, taken with a Leica light microscope. (a) Ventral view of P. lawrencei s.s. (b) Ventral view of the novel lineage, P. aereus sp. nov. (clade B). Enhanced ventral views showing the ventral organs in (c) P. lawrencei s.s. and (d) the novel lineage, P. aereus sp. nov. Scale bars: 2 mm.

Fig. 7.

Scanning electron micrographs of dermal papillae on dorsal and ventral surfaces of P. lawrencei s.s. and P. aereus sp. nov. Black dots represent a single scale rank. (a) Dorsal papilla for P. lawrencei s.s. with seven scale ranks. (b) Dorsal papilla for the novel lineage (clade B, P. aereus sp. nov.) with 15 scale ranks. (c) Ventral papilla for P. lawrencei s.s. with four scale ranks. (d) Ventral papilla for the novel lineage (clade B, P. aereus sp. nov.) with nine scale ranks. Scale bars: 100 μm.

Novel species delineation

Using the ‘unified species concept’ advocated by de Queiroz (2007), we find P. lawrencei s.s. (clade A) and the sympatric lineage at the Riviersonderend Mountains (clade B) to be separately evolving metapopulations, with multiple lines of evidence substantiating the classification of clade B as a divergent species. Following conventional taxonomic criteria, this can be identified as a unique species based on discernible morphological differences (Fig. 6a–d and 7a–d, Table 4). Using both gross morphology and scanning electron microscopy the two clades within P. lawrencei s.l. could be delineated. In clade A (P. lawrencei s.s.), the ventral colour was highly variable (Table 4) however in clade B, based on two samples, the ventral surface colour was bronze. In velvet worm species, the dorsal colour pattern can notably be highly variable. For example, in Peripatopsis polychroma the dorsal colour surface ranged from olive green, through slate black to rusty brown (Grobler et al. 2023). Similarly, in P. overbergiensis, the dorsal surface colour can range from rusty brown to bright red (Daniels pers. obs). Although leg pair count is a variable indicator of species boundaries (Daniels et al. 2009, 2013, 2016; Barnes et al. 2020), the number of leg pairs differed between the species in our study, as 17 were present in P. lawrencei s.s. and 18 in clade B. McDonald et al. (2012) reported that in P. capensis and P. overbergiensis the number of leg pairs ranged from 17 to 18 respectively, with limited overlap in this character. By contrast, Grobler et al. (2023) reported that between sister species P. birgeri and P. polychroma the number of leg pairs always overlapped and ranged from 20 to 23. Leg pair numbers, together with the colour of the dorsal surface should be interpreted cautiously and preferably within a phylogenetic context to limit artificial splitting or lumping of species. The arrangement of dermal papillae has more often been used as a diagnostic feature in earlier velvet worm studies (Oliveira et al. 2011; McDonald and Daniels 2012; Daniels et al. 2013, 2016; Sato et al. 2018; Grobler et al. 2023). The novel lineage (Oubos B, clade B, P. aereus sp. nov.) was shown to have a reduced number of scale ranks on the papillae for both the dorsal and ventral surfaces.

Incorporating DNA sequence data into our argument for the recognition of a new species, we observed several lines of evidence that aid in the delineation. Marked uncorrected sequence p-distance values, using the COI locus between the two clades falls into a similar range as reported among velvet worm species. For example, the uncorrected p-distances between P. capensis, P. overbergiensis and P. lawrencei s.s. ranged between 5.93 and 7.49%, however, in our study the uncorrected p-distance COI distance value between P. lawrencei s.s. and the novel lineage (clade B) was between 9.13 and 10.62%. Similarly, between clade B and P. overbergiensis, the uncorrected p-distance was between 10.57 and 10.73%. The values found between clade B and P. lawrencrei s.s. are higher than what has been reported for sister species of velvet worms in prior studies (McDonald and Daniels 2012), suggesting a more distant relationship; a result corroborated by our phylogeny for Peripatopsis (Fig. 5). Barnes et al. (2020) used similar p-distance values ranging from 8.46 to 10.64% to delineate velvet worm species. Sequence divergence values as low as 3% have been used to delineate novel velvet worm species (Rockman et al. 2001). In addition, the 5.58% uncorrected p-distance difference between the two sympatric phylogenetically distinct Oubos lineages (A and B) observed for the highly conserved 18S rRNA locus further corroborates the genetic distinction of the two velvet worm species. Barnes et al. (2020) also used uncorrected p-distance differences in the 18S rRNA locus to validate the recognition of species in the P. clavigera species complex. For the 18S locus the latter authors noted p-distance values ranging from 0.68 to 3.40% to aid the recognition of novel species. The 18S rRNA p-distance value we report is higher than those from earlier velvet worm studies, providing corroborative nuclear differences for the recognition of a novel lineage. To further strengthen our argument for distinct taxonomic status, a multi-locus tree topology (incorporating both mt and nuDNA sequence data) of the South African Peripatopsis suggests a Miocene divergence between the two Oubos clades (A and B) (Fig. 5). Given the multifaceted evidence for the distinction of this novel lineage in clade B, we are confident in the delineation and recognition of the new species.

The divergence time estimation indicates an Oligocene origin of Peripatopsis, 28.83 Ma (95% HPD: 21.27–37.49 Ma), followed by rapid diversification in the mid-Miocene to the Pleistocene. The Miocene divergence of the stem lineage, 13.68 Ma (95% HPD: 10.46–17.73 Ma) suggests a more ancient origin of this velvet worm and a wider distribution than previously considered (Daniels et al. 2009). This can be attributed to the different analytical approaches employed by the two studies. The Miocene epoch was characterised by humid subtropical forests in the Cape region of southern Africa (Steinthorsdottir et al. 2021) but was succeeded by major climatic ameliorations such as the development of the proto Benguela current along the west coast. This resulted in increased aridification, possibly contributing to cladogenesis (Siesser 1980). Cladogenesis was also observed at a comparable time period, 13.63 Ma (95% HPD: 9.53–15.73 Ma), between the sister taxa P. bolandi and P. balfouri (Daniels et al. 2013). Notably, these have a similar distribution, encompassing the Western Cape region and would have experienced similar climatic pressures. The dry conditions of the late Miocene persisted throughout the Plio–Pleistocene (Mucina and Rutherford 2006; Sepulchre et al. 2006). Contemporary populations of Peripatopsis generally reside at high altitudes, suggesting biome contraction in lowlands. This could have driven many lineages upward to forest refugia, further separating populations. As velvet worms are extremely mesic-adapted animals, the long-distance gene flow is highly restricted. In such cases, local adaptation or genetic drift is expected to play an important role in species divergence and genetic differentiation.

Phylogeography of P. lawrencei s.s.

The phylogeographic analysis of P. lawrencei s.s. using the COI dataset yielded results congruent with the tree topologies observed thus far (Fig. 2 and 5). Velvet worms are known to be highly fragmented and generally restricted to high elevation forest patches (Daniels et al. 2009, 2016; Barnes and Daniels 2019; Barnes et al. 2020; Grobler et al. 2023). The observed distribution of P. lawrencei s.s. is consistent with such habitat restriction. Haplocluster A3 had a wide distribution, A2 was confined to the Riviersonderend Mountains and A1 showed exclusivity to a narrow coastal region at Hermanus (Fig. 2 and 3). The geographically isolated haploclusters underscore the presence of strong barriers to gene flow and limited maternal dispersal among populations, as observed in prior velvet worm studies (Daniels and Ruhberg 2010; Barnes et al. 2020; Barnes and Daniels 2019, 2022). There is a clear interplay of geographical factors, such as physical barriers or habitat fragmentation and genetic differentiation across these populations. This distribution may have been influenced by localised extinctions in neighbouring or connecting populations, resulting in the observed genetic drift and allelic fixation. The F_ST values further support this, showing considerable genetic differentiation among the conspecific populations, indicating high genetic structuring and low rates of dispersal (Table 2). The lack of connectedness promotes genetic divergence, as supported by the AMOVA results, highlighting that a substantial portion of genetic variation occurred among populations. The higher diversity values in certain localities such as Greyton and Linden Hill could reflect the presence of refugia, allowing for accumulation of genetic diversity (Table 3). These findings suggest the need for conservation units within Peripatopsis lawrencei s.s., as the Afromontane habitats are highly reduced and highly threatened. The scarcity of samples for the novel lineage further supports the need for habitat maintenance.

Fine-scale sampling and future directions

Fine-scale sampling plays a crucial role in uncovering hidden alpha-taxonomic diversity, especially in poorly studied saproxylic taxa. These organisms are often overlooked or understudied due to habitat complexity or crypticity. Saproxylic zones have been observed to harbour a wealth of taxonomic diversity, as these comprise numerous microhabitats allowing for highly specific niche divergence among taxa that may live in proximity (Daniels et al. 2009; Oliveira et al. 2011; Sato et al. 2018). Furthermore, weak taxonomic classifications can arise through poor sampling, leading to divergent sympatric species being lumped together into a single clade (Barnes and Daniels 2019). Many saproxylic taxa are rare and potentially endangered due to the nature of the habitats that are prone to habitat degradation and fragmentation. Fine-scale sampling allows researchers to detect these rare species and assess the conservation status precisely. A study by Barnes and Daniels (2019) demonstrated the difference between large spatial scale- and fine-scale sampling in detecting novel lineages. Fine-scale sampling, as undertaken here, enables refinement of taxonomic classifications where diversity has previously been underestimated. In our case, fine-scale sampling revealed an ancient novel lineage (clade B) living in sympatry with a well-documented species, P. lawrencei s.s. in the absence of these detailed studies, such species might be unnoticed, impeding conservation efforts. A fine-scale approach is highly recommended when investigating species of limited dispersal capacity and those that inhabit saproxylic zones. Future studies could intensify field surveys and assess the conservation status of this novel lineage. Population size estimates and ongoing monitoring of this lineage and the populations of P. lawrencei s.s. can aid in evaluating the vulnerability and informing targeted conservation strategies, including the identification of potential conservation units. Genome-wide analyses could be conducted to pinpoint specific genes responsible for morphological divergence, while also providing insights into broader evolutionary and developmental processes in velvet worms within these rapidly changing landscapes.

Family PERIPATOPSIDAE Bouvier, 1907

Genus Peripatopsis Pocock, 1894

(Fig. 2, 5, 6a, b, 7a, b, Table 4.)

Peripatopsis lawrencei s.s. McDonald, Ruhberg & Daniels, 2012.

Peripatopsis aereus Daniels & Nieto Lawrence, sp. nov.

(Fig. 2, 5, 6c, d, 7c, d, Table 4.)

ZooBank: urn:lsid:zoobank.org:act:71276915-79C6-4759-9244-20BEF6303F58