Genetic and morphological analyses indicate that the Australian endemic scorpion Urodacus yaschenkoi (Scorpiones: Urodacidae) is a species complex

Biological material

Samples of Urodacus yaschenkoi were obtained from field and museum collections (Table 1). Live specimens were collected from eight locations (approximately 500 m²) in the semi-arid and arid regions of Central Australia in December 2010 and October 2011 (Table 1 and Fig. 1). Individuals were collected at night from pitfall traps set in front of their burrows, and those outside their burrows were detected using ultraviolet (UV) lamps that reveal soluble fluorescent components (β-carboniles) in the scorpion exoskeleton (Stachel, Stockwell & Van Vranken, 1999). Captured scorpions were kept alive and transported to the laboratory for morphological identification according to Koch (1977). Key diagnostic feature that distinguishes U. yaschenkoi from other Urodacus species is a very small terminal prolateral tarsus unguis. All specimens were handled according to good animal practices defined by the Government of Australia, and all institutions and museums involved approved the animal handling work. Scorpions were anaesthetized by cooling in a refrigerator (4 °C) for 5 min before removing ∼1 mm² of leg muscle tissue, which was stored in 90% ethanol at 4 °C or −20 °C for subsequent DNA extraction. Additional samples were obtained from collections at the South Australian Museum (SAM) and Western Australian Museum (WAM) containing specimens collected between 2000 and 2010 (Table 1).

DNA extraction, amplification and sequencing

Total DNA was extracted from the stored muscle tissue using the DNeasy Blood and Tissue Kit (Qiagen, Venlo, Netherlands) following the manufacturer’s instructions. Two mitochondrial loci (cytochrome oxidase subunit I, COXI; large ribosomal subunit, 16S) and a single nuclear locus (28S) were amplified by PCR with a reaction volume of 20 µl containing 0.5 ng of template DNA, 10 µl of Go Taq Master Mix (Promega, Madison, Wisconsin, USA), 0.5 µl of 10 nM primers and 7 µl of RNase-free water (Qiagen). The primer sequences and PCR amplicon sizes are summarized in Table 2.

Primers previously designed for the insect COXI gene (Simon et al., 1994; Tanaka et al., 2001) were used to amplify a 630-base pair (bp) fragment from the 3′ end of the locus. The amplification conditions comprised an initial denaturing step at 95 °C for 5 min followed by 35 cycles of denaturing at 94 °C for 30 s, annealing at 52 °C for 40 s, and extension at 72 °C for 45 s, and a final extension phase at 72 °C for 5 min. For the mitochondrial 16S gene, the scorpion-specific primer pairs modified by Gantenbein et al. (2005) were used to amplify a 425-bp region at the 3′ end of the locus. The amplification conditions comprised an initial denaturing step at 94 °C for 4 min followed by 30 cycles of denaturing at 94 °C for 30 s, annealing at 47.5 °C for 30 s, and extension at 72 °C for 30 s, and a final extension phase at 72 °C for 7 min. The COXI and 16S gene fragments were also amplified from three specimens keyed out as Urodacus manicatus (Um2714, Um1814) and U. novaehollandiae (Un2112, Table 1). Sequences from these taxa were used as outgroups in downstream phylogenetic reconstruction. Primer pairs R1S and R1AS, and R2S and R2AS, designed by Arabi et al. (2012), were used to amplify 1,158-bp and 1,246-bp fragments of the 28S locus, respectively. Each set of primers amplifies a different region of the gene, which overlaps by 327 bp, and their sequences were concatenated to form a larger product of 2,076 bp. The amplification conditions for both sets of primers comprised an initial denaturing step at 94 °C for 4 min, followed by 30 cycles of denaturing at 94 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s, and a final extension phase at 72 °C for 7 min.

Museum specimens that were not stored under ideal conditions for preservation failed to yield COXI amplicons suitable for direct sequencing. To address this issue, additional PCR primers were designed to amplify smaller fragments for COXI locus (Table 2), resulting in amplicons of 150 bp that were used for subsequent analysis. For the SAM specimens, the amplification of the 28S nuclear gene failed entirely and these samples were excluded from further analysis of the nuclear gene variation. All amplicons were sequenced in both directions using the PCR amplification primers, and carried out on an Applied Biosystems 3130 genetic analyzer by Macrogen Inc. (Seoul, South Korea).

Sequences were aligned and edited in Geneious Pro v6.1 (Biomatters Ltd.) using the MUSCLE alignment option with default parameters. All chromatograms were checked for the presence of multiple peaks (which indicate heterozygosity), and authenticity of the COX1 coding gene was validated by checking for indels and premature stop codons. After this editing process, the alignment of the mitochondrial gene fragments yielded 616-bp and 396-bp products for the COX1 and 16S genes respectively, and the final 28S alignment was 2,076 bp in length. The final dataset contained 68 sequences for each of the mitochondrial genes and 27 sequences for the 28S locus (Table 1, (GenBank accession # KP176717 –KP176786)). Shared haplotypes were identified and the uncorrected pairwise genetic distances (%) were calculated using Geneious Pro v6.1 (Biomatters Ltd.). This simple distance measure was implemented to achieve reliable estimates of both intraspecific and interspecific genetic variation.

Phylogenetic analysis

Phylogenetic reconstructions and divergence dates among lineages were calculated using unique haplotypes and Bayesian Inference (BI) methods implemented in BEAST v2.1.3 (Bouckaert et al., 2014). We used jModeltest v0.1.1 (Posada, 2008) to select the best-fit model of evolution, based on Akaike Information Criteria (AIC) (Akaike & Company, 1981) for each of the mitochondrial and nuclear genes (GTR + G in each case). Mitochondrial loci were combined for analysis due to their similar modes of evolution (GTR + R), as indicated by the incongruence-length difference (ILD) tests (Farris et al., 1995) implemented in PAUP_4.0b10 (Swofford, 2002). The nuclear gene (28S) was analyzed independently due to inconsistencies in taxon sampling (Table 1).

Operators were auto-optimized, and five independent Markov Chain Monte Carlo (MCMC) runs were performed using a Yule (speciation) tree-prior, each running for 5 ×10⁶ generations, sampling every 10,000 states. Log files were examined with Tracer v1.5 (Rambaut et al., 2014) to ensure that runs were sampling from the same posterior distribution, to determine appropriate burn-in, and to ensure that effective sample sizes (ESSs) of parameters of interest were greater than 1,000. Tree files of independent runs were then combined using LogCombiner v2.1.3 (Drummond et al., 2012), discarding the first 20% and re-sampling at a lower frequency of 15,000. The maximum clade credibility (MCC) tree was recovered from a sample of 10,000 posterior trees, and branch support was annotated using TreeAnnotator v2.1.3 (Drummond et al., 2012). Each analysis started with a random starting tree and seed with no root specified. Sequence data from species of the same genus (U. manicatus and U. novaehollandiae) were used to estimate the root of the mitochondrial gene tree.

Additional phylogenetic constructions were also performed using a truncated COXI alignment to test the influence of missing data on the final tree topology. Because numerous museum collections yielded short COXI gene products, we trimmed the alignment to 150-bp to exclude regions of the alignment with high levels of missing data. This exercise demonstrated that the inclusion/exclusion of missing data had little influence on the phylogenetic reconstructions. Consequently, all results presented from this point reflect those from the non-truncated COX1 alignment.

Species delineation based on molecular data

We implemented two DNA taxonomy approaches to evaluate the presence of cryptic species. First, the general mixed Yule coalescent (GMYC) approach (Fujisawa & Barraclough, 2013; Pons et al., 2006) was applied to an ultrametric tree (produced using BEAST) in R Development Core Team (2008) with the Splits package (http://splits.r-forge.r-project.org). The GMYC model is a process-based approach that detects the threshold in a gene tree at which within-species processes (i.e., coalescence) shift to between-species processes (i.e., speciation and extinction). Second, we combined the Poisson Tree Processes model for species delimitation (PTP) and a Bayesian implementation of PTP (bPTP) to infer putative species boundaries on a given phylogenetic input tree (Zhang et al., 2013). The PTP/bPTP model, unlike the GMYC model, requires a bifurcated phylogenetic tree rather than an ultrametric tree. PTP/dPTP models speciation or branching events in terms of the number of substitutions. The following parameters were used: MCMC, 500,000 generations; thinning, 100; burn-in, 0.1; seed, 123, and assessed convergence in each case to ensure the reliability of the results.

Delineation based on the analyses of morphological measurements

Proportions of several characters that approximate body shape were assessed in 39 female adult specimens that were keyed out as U. yaschenkoi (according to Koch, 1977) and were collected at 26 locations (Table 1, Fig. 1). Gender was determined by examining the genital opercula of adult scorpions, with males having a small finger-like projection known as the genital papilla. Because our collection contained only three males, the analyses were done only with females.

The following traits were measured under a microscope using an ocular ruler with 1-mm precision: carapace length (CL), metasoma segment V length (MVL), telson length (SL), pedipalp length (PL), chela length (ChL), pecten length (PecL) and pecten width (PecW). Ratios of traits (e.g., CL/MVL, SL/PL etc.) gave in total 21 variables scored in each individual (Supplemental File 4). These variables were treated as predictors in the Linear Discriminant Analysis (LDA) implemented in the R package “MASS” (Venables & Ripley, 2002). LDA was used to test whether the linear combination of 21 variables (ratios of morphological measurements) can predict individual’s membership to a mitochondrial lineage (putative species). Strong predictive power of morphological variation on the observed molecular divergence would provide additional support for a species complex in U. yaschenkoi.

Divergence time estimation

The mitochondrial gene tree was time calibrated with divergence times of nodes inferred from 95% highest posterior density (HPD) intervals. Scorpion-specific mutation rates of 0.007 substitutions/site/million years for COXI and 0.005 substitutions/site/million years for 16S (Gantenbein et al., 2005; Gantenbein & Largiadèr, 2003) were used to calibrate the tree. These estimates are derived from buthid scorpions and have been used to estimate divergence times among various scorpion lineages including non-buthid taxa (Bryson et al., 2013a; Bryson et al., 2013b; Graham, Oláh-Hemmings & Fet, 2012). Substitution rates were set in BEAUti v1.7.3 (Drummond et al., 2012) using relaxed clock log normal priors. Tracer was then used to obtain parameter estimates for time to most recent common ancestor (tMRCAs) for nodes within the gene tree.

Phylogenetic analysis

Mitochondrial markers

Bayesian inference analysis of the mitochondrial dataset identified several genetically divergent lineages (three major lineages represented as black, red and green clades in Fig. 2), with strong statistical support for their respective monophyly (posterior probability > 0.95). Sublineages within the black clade are broadly distributed across Victoria, South Australia and Western Australia, whereas the red and green clades are restricted to Western Australia (Fig. 1). From this point forward we will refer to the black, red and green clades as the south-central (SC), western (W) and central-western (CW) lineages, respectively.

Mean uncorrected pairwise genetic distances between the three major lineages (SC, CW and W) ranged from 6.4 to 6.9% (overall mean ± standard deviation = 6.6% ± 0.9%). The mean sub-lineage distances ranged from 2.2% ± 0.4% and 0.8% ± 0.2%, respectively (not calculated for the W lineage due to only a single recorded haplotype). Mean uncorrected distances between the three major lineages and the outgroups ranged from 9.3 to 10.3% (mean ± standard deviation = 9.4% ±  1.4%).

Nuclear marker

Despite low level of variation in the 28S dataset, Bayesian analysis produced a nuclear gene topology that was largely concordant with the mitochondrial gene tree. Three genetically divergent clades were identified, corresponding to those from the mitochondrial dataset (SC, CW and W, Fig. 3). In each case, strong statistical support for the monophyly of each clade was found (posterior probability >  0.95). The unresolved interrelationships among lineages within each clade in the nuclear gene tree prevented any reliable inferences of phylogeographic patterns.

Molecular-based species delineation

Among the 31 unique mitochondrial haplotypes described above, the GMYC model identified nine entities and the PTP/bPTP approach identified seven, each representing putative species (Table 3). The assignment of haplotypes to putative species groups is shown in Fig. 2, where conspecifics share a common number. Species assignments were highly consistent when comparing each of the methods, but we presented the PTP/bPTP results as they are more accurate when the evolutionary distances between lineages are small (Zhang et al., 2013). In summary, SC, W and CW clades were recognized as putative species groups, as were the sub-lineages within the SC ancestral grouping (SC-1 to 5, Fig. 2).

Discriminant power of morphological variation

None of the U. yaschenkoi specimens that were characterized at 21 morphological ratio variables were assigned to the W mitochondrial clade, hence the LDA was done on 39 females assigned to the SC and the CW clades. Individuals were categorized into four groups (putative species) based on the results of the PTP/bPTP molecular species delineation analysis: 18 females from SC-1, 12 from SC-3, three from the SC-4, and six from the CW clade (Fig. 2). Because our dataset contained four groups, we could find a maximum of three discriminant functions that separate these groups.

The first discriminant function (LD1) achieved 93.7% of the separation, reflecting the morphological distinction of the CW clade from the SC clade (Fig. 4). Further separation of the three putative groups within the SC clade was weak (LD2-3, Fig. 4). We then grouped samples into two putative species (CW and SC clade) and tested the accuracy of prediction using 100 jackknife resampling steps. The grouping into two molecular clades based on morphological variation was 100% accurate (33/33) for the SC clade and 83.3% accurate (5/6) for the CW clade. Therefore, our results indicate strong predictive power of body proportion variation on the observed molecular divergence, and suggest the existence of at least two distinct taxa within U. yaschenkoi.

The most discriminating uncorrelated proportions were of the telson and chela length (SL/ChL) and pedipalp and pecten length (PL/PecL). Overall, members of the CW clade tend to have disproportionately shortened chela and enlarged pecten when compared to the members of the SC clade.

Divergence dating

Our time calibrated mitochondrial phylogeny suggested that the split between the major U. yaschenkoi clades (SC, CW and W lineages) occurred during the late Miocene/early Pliocene (4–7 MYA) (Fig. 2). Lineage diversifications within SC appear to have occurred during the Pliocene and early Pleistocene (1.8–4 MYA), while finer-scale phylogeographic patterning within the sub-lineages arose during the late Pleistocene (<1 Mya). Divergence time estimates should be interpreted with some caution, as the nucleotide substitution rate was derived from a different scorpion family (Buthidae) and there are large errors margins around 95% HPD estimates.

Biogeographic patterns

The SC lineage showed substantial geographic structure. The most divergent sub-lineage (SC-5) was found in Western Australia in sympatry with the CW lineage (Fig. 1). SC-1 was found west of the Central Ranges, through to the Eyre Peninsula in South Australia, while SC-3 had a distribution extending from the Central to Mt Lofty Ranges in South Australia, and across to north-western Victoria. SC-4 had a narrow north-south distribution in the central inland and coastal regions of South Australia (Fig. 1).

Revising the U. yaschenkoi taxonomy—future directions

Our results provide solid baseline data on the historical and spatial extent of diversification in U. yaschenkoi and offer some guidelines for future integrative taxonomic approaches in delimiting species within this taxon. We found an agreement among disciplines (morphology, nuclear and mitochondrial genetic information) during a primary exploration, which strengthens the argument for a taxonomic revision (Pante, Schoelinck & Puillandre, 2014; Schlick-Steiner et al., 2009). Congruent morphological and molecular phylogenetic signals are particularly compelling for a scorpion taxon, given that this is not the case in many scorpion lineages (Sharma et al., 2015).

The level of mitochondrial sequence divergence observed between U. yaschenkoi lineages satisfy the requirements for species delineation based on the principles of the phylogenetic species concept (De Queiroz, 2007; Wheeler, 1999), The three major lineages (SC, CW, W) can be considered the putative species. Because genetic ‘yardstick’ approaches provide crude taxonomic measures and nucleotide substitution rates often vary considerably between taxonomic groups, some caution is needed when considering findings of these analyses alone. Additional DNA-based species delineation approaches (GMYC and bPTP) indicated extensive cryptic speciation in U. yaschenkoi (Fig. 2). The GMYC method has been criticized for over-splitting species with a pronounced genetic structure (Satler, Carstens & Hedin, 2013), yet several recent studies have shown that it is highly robust (Fujisawa & Barraclough, 2013; Talavera et al., 2013) The obvious next step is to characterize the nuclear genome-wide variation in U. yaschenkoi sampled extensively within the “type” locality (28°35′S, 138°33′E), as well as western and northern parts of the distribution. We certainly advise against a pool-sequencing phylogenomic approach (e.g., samples from the same location are pooled to achieve cost-efficiency), given that the putative species have been found in sympatry.

The proportions of various morphological characters are routinely used in species descriptions or identification keys, particularly for arthropods where morphologically similar species often differ significantly in body proportions but not in qualitative characters (Baur & Leuenberger, 2011). Arguably, the results of multivariate analyses summarizing the overall body shape differences between groups are not easily interpreted. Yet, our initial results suggest that further analyses of e.g., chela shape might reveal more easily quantifiable diagnostic characters for U. yaschenkoi. Several parameters of chala shape were found to be correlated with the amount of strain stress they can withstand. Specifically, slender chela morphologies may be less suitable for high-force functions such as burrowing and defence (Van der Meijden, Kleinteich & Coelho, 2012). Given that U. yaschekoi putative species (SC and CW) show marked shape differences involving chela, further exploration of burrowing behavior or pray preference might provide additional characters to describe the U. yaschenkoi species complex.

Finally, it is important to note that we cannot exclude the possibility that some of the cryptic lineages have already been described as species, and we are not able to compare our genetic data against other Urodacus sequences as none published at the time of our study. Also, our sampling did not cover the exact “type” locality (28°35′S, 138°33′E). The samples closest to this area belong to the SC clade and likely represent the “type” lineage These data gaps would need to be addressed in further studies aiming to revise the taxonomy of the Australian desert scorpion U. yaschenkoi.