Taxonomy is a key discipline for describing and understanding biodiversity and is critical for documenting and inventorying undescribed species before they reach extinction due to climate change or habitat loss (Wheeler 2018, 2020). Even though taxonomy has been recently recognized as extremely important (Bond et al. 2022), it has been regarded by some as less relevant and with little intellectual content, resulting in what has been called “taxonomy in crisis” (Agnarsson & Kunter 2007). It is true that many taxonomic groups, such as mygalomorph spiders (e.g., Rix et al. 2017; Harvey et al. 2018; Opatova et al. 2020) often remain poorly resolved due to the occurance of multiple cryptic species that offer little morphological variation that can be used by taxonomists in species description. In this study, we show how multivariate statistical analysis of trait shape can assist in species delimitation in taxonomic studies of mygalomorph spiders.

The genus Proshermacha Simon, 1908 (Mygalomorphae: Anamidae) was recently resurrected (Harvey et al. 2018) to include nine species (World Spider Catalog 2021) from southern Australia (Main 1982; Raven 2000). Members of the genus are relatively gracile mygalomorph spiders that construct silk-lined burrows with an open entrance (Harvey et al. 2018). Molecular data indicate the existence of numerous undescribed species within its distributional range (Harvey et al. 2018). One locality of interest for the genus is the Stirling Range, a mountain formation estimated to have originated ∼1.2 billion years ago (Rasmussen et al. 2002). The Range is home to ancient endemic invertebrate species including land snails, onychophorans, assassin spiders, and trapdoors spiders (Rix et al. 2015), all of which are considered relictual groups that can be traced back to the Gondwanan super-continent (Cooper et al. 2011). It has been suggested that biodiversity within the Range has been driven by multiple mechanisms of diversification, such as the remnants of relictual fauna, vicariant isolation, or in situ speciation (Rix et al. 2015). Several studies have explored the biogeography and evolution of Stirling Range fauna. For instance, five species of the millipede genus Atelomastix Attems, 1911 (Diplopoda: Spirostreptida: Iulomorphidae) (Edwards & Harvey 2010), three species of the genus Bertmainius Harvey, Main, Rix & Cooper, 2015 (Araneae: Migidae) (Harvey et al. 2015), and four species of the trapdoor spider genus Cataxia Rainbow, 1914 (Araneae: Idiopidae) (Rix et al. 2017) have been identified and are known to be endemic to the Range. In each of these genera, individual species are found on different peaks within the range, supporting the idea that the Stirling Range acted as a natural refuge for relict invertebrates and, given its isolated sky-island nature, the peaks represent an important driver of allopatric speciation for species with limited-dispersal (Main 1993). Although little is known regarding the habitat requirements, biogeography, and speciation patterns of the genus Proshermacha, like all mygalomorph spiders they are likely to be dispersal-limited (Buzatto et al. 2021), and we might therefore expect to find a similar pattern of vicariant isolation in the Stirling Range.

Extracting both qualitative and quantitative data from morphological characters of mygalomorph spiders can often be difficult. Morphological traits might suggest geographically isolated populations that represent a single, undifferentiated species, which when analyzed through molecular methods result in distinct molecular species (Bond & Stockman 2008). Molecular data has without doubt contributed to an increasing appreciation of the extent of biodiversity (Franzini et al. 2013; Castalanelli et al. 2014; Teixeira Jr. et al. 2016; Hupalo et al. 2020) by recognizing molecular divergence when there are seemingly few morphological differences between taxa (Bond et al. 2001), and by solving delimitation of species complexes (Bond & Stockman 2008) or differentiating between cryptic species (Leavitt et al. 2015). Molecular data have helped develop modern-day systematic classifications of mygalomorph spiders (e.g., Bond et al. 2012; Opatova et al. 2020; Harvey et al. 2018, 2020) that provide a better understanding of the evolutionary relationships amongst groups and lineages.

Taxonomy has improved greatly since its initial species descriptions based on comparative morphology, with new and complementary methods being develped for delimiting species boundaries (Dayrat 2005; Padial et al. 2010). Although molecular sequencing provides more reliable, reproducible, universal, and scalable data (Sharma et al. 2021), both natural history and morphological data are crucial, and should be used in combination to support species hypotheses derived from molecular data (Giribet 2015; Muster & Michalik 2020). Traditional and modern approaches should be combined into an integrative taxonomy (Bond et al. 2022). The role of morphology becomes even more evident for those species groups that have not been revised and are known only from museum material which can not be analyzed with molecular methods due to inadequate preservation (Derkarabetian et al. 2019). In the field of arachnology, most molecular studies work in tandem with qualitative morphological descriptions (with some exceptions; see Bond 2012), but few incorporate multivariate statistical methods during data analysis to support morphological differences between molecularly distinct species (Wilson et al. 2021).

Geometric morphometrics (GM) is a technique that relies on landmarks placed over specific anatomical structures from which quantitative shape variation can be obtained (Bookstein 1991; Rohlf & Marcus 1993). Soon after its development (Bookstein 1991), GM was adopted by numerous biologists and has improved morphological studies by allowing multivariate analysis to be conducted to capture shape variation (Adams et al. 2004). GM has been used widely in taxonomic groups such as vertebrates, reinforcing species delimimations based on mutilocus data in milksnakes (Ruane 2015), to maximize shape coverage in 2D models of varanid lizards (Openshaw et al. 2017), and even to explore optimal sample sizes and sampling error when utilizing GM using horse teeth as a model (Cardini 2014). Within arachnids GM has been used to compare phylogenetic and morphological divergences, reveal sexual dimorphism in shape allometry, quantify intraspecific variation, and to delimit species within a species-complex (Crews & Hedin 2006; Fernádez-Montraveta & Marugán-Lobón 2017; Torres et al. 2018; Wilson et al. 2021). Like many fields, GM has grown rapidly with new techniques becoming available such as the use of semi-landmarks that allow capturing curvature data (Adams et al. 2013), the development of specialized programs for landmark acquisition and manipulation (Rohlf 2013, 2017, 2021), and specific statistical packages to analyze landmark data (Adams et al. 2016). The use of GM provides taxonomists a quantitative tool with which to test for morphological differences among species (Wilson et al. 2021). Although there is the added cost of greater effort and identification time, quantitative shape analysis may be the only means available when discrete morphological features are unavailable. We had access to numerous alcohol-preserved specimens belonging to the genus Proshermacha that were collected during the early 1990s. These samples provided an ideal study group for applying quantitative geometric morphometric analysis to test for distinct morphotypes and their distributions through the Stirling Range.

METHODS

Study site.—Samples used in this study came from five different localities in the Stirling Range National Park in southwest Western Australia (Fig. 1). Mondurup and Talyuberlup Peaks, in the western and southern parts of the range, respectively, are dominated by thicket plant formations in the higher areas, and mallee and woodlands at lower elevations. Mount Magog and Toolbrunup, in the northern and eastern parts of the range, reach higher elevations, and both are dominated by a combination of thicket and mallee (Keighery 1993). Two Mile Lake is located on the lower elevations at the eastern edge of the park and is dominated by open landscapes with remnants of bushes and short mallee (Keighery 1993).

Figure 1.

Map of morphotype distributions and sampling localities. ☉ represents Morph1 and ☆ represents Morph2. Notice sympatry of morphotypes at Mondurup, Mt. Magog (separated altitudinally), and Toolbrunup.

Material examined.—Samples were collected in pitfall traps with ethyl-glycol (70-30 mix) as preservation liquid and stored in 70% ethanol in the Western Australian Museum (WAM). Given the lack of female samples across all localities, only adult males were used for this study. A total of 39 specimens from an undescribed species were analyzed, coming from the localities as follows: Mondurup n = 5; Mt. Magog Low n = 4; Mt. Magog High n = 4; Talyuberlup n = 5; Toolbrunup n = 15; Two Mile Lake n = 6. Across all specimens, the right palpal bulb was separated from the cymbium (Figs. 2A–C), and the left leg I was detached from the body for capturing images of the tibia (Figs. 2E,F) and metatarsus (Figs. 2H, I) to allow for accurate landmark placement. Structures were fully submerged in 70% ethanol on a petri dish and held in place by colorless hair styling gel when capturing images, orienting the bulb to capture the ventral side and the retrolateral side for leg I.

Figure 2.

Visually distinct types of bulbs, ventral view (A–C); tibia, prolateral side of leg (E–F); metatarsus, prolateral side of leg (H–I); and detail of landmark placement (D, G, and J). Arrow in B points to the flanged process at the tip of the embolus. Scale bar represents 1mm for all images.

Image acquisition was conducted with the Leica Application Suite v4.6 using the multi-focus setting, visualized through a Leica DFC 500 digital camera attached to a Leica MZ16A microscope. Image scaling and landmark digitalization were performed using tpsDig2 v2.3.1 (Rohlf 2017) and following a counterclockwise direction for each anatomical structure. Landmarks consisted of Type I fixed and Type II sliding semi-landmarks and followed the typology of Bookstein (1991) and Zelditch et al. (2004) selection criteria and are described in Table 1 (see Fig. 2D, G, and J for landmark placement).

Table 1.

Landmark placements and their contribution to shape variation for each structure. Type I are fixed landmarks, Type II are sliding semi-landmarks (Bookstein 1991, Zelditch et al. 2004).

Morphometrics and statistical analysis.—Landmark coordinates were aligned and superimposed by Generalized Procrustes Analysis (GPA) with the software tpsRelw (Rohlf 2013) which provided data including centroid size (CS), partial warps (PW), and relative warps (RW). Measurements of cephalothorax length were obtained using the linear measurement tool in the Leica Application Suite v4.6 and were used as the standard measure of body size (Size) for each sample. Repeatability analyses were conducted for each of the three traits (bulb, tibia, and metatarsus) to determine if the selected landmarks provided repeatable estimates of trait size and shape. Landmarks were placed on images from 10 randomly selected individuals (2 per site), on 3 different occasions. Relative warps and centroid sizes were extracted and the repeatability analyses were conducted as Gaussian data type with 1000 bootstrap iterations using the rptR package (Stoffel et al. 2017).

Multivariate analysis of variance (MANOVA) was conducted using those RWs that represented 95% of the cumulative shape variance of each structure. Locality was used as the explanatory variable, and both Size and CS were used as covariates to control for body size and trait size respectively during analysis. Pairwise comparisons between populations were made using only RW1 for each structure as it accounted for most of the shape variation (. 60%). All values were set to 0.05 significance and were conducted in RStudio (RStudio 2019) using R v.4.0.5 (2021).

RESULTS

Shape analysis.—Repeatability analysis showed that landmark placement provided repeatable estimates for bulb centroid size (P ,0.001, R = 0.98 [CI = 0.941, 0.993]) and shape RW1 (P ,0.001, R = 0.96 [CI = 0.876, 0.986]), tibia centroid size (P ,0.001, R = 0.98 [CI = 0.955, 0.996]) and shape RW1 (P ,0.001, R = 0.988 [CI = 0.957, 0.995]), and metatarsus centroid size (P ,0.001, R = 0.99 [CI = 0.978, 0.998]) and shape RW1 (P ,0.001, R = 0.99 [CI = 0.979, 0.998]). For the three structures, RW1 explained most of the trait shape variance (Table 2). Visualisation of RW1 and RW2 revealed the existence of two clusters of specimens for bulb (Fig. 3A) and tibia (Fig. 3B), each of which contained specimens from multiple localities, but no clear clusters could be observed for metatarsus (Fig. 3C).

Table 2.

Percentage of variance in trait shape explained by each Relative Warp (RW). Only those RW corresponding to .95% of trait shape variance per structure are included.

Figure 3.

Scatterplots showing positive and negative values for RW1 and RW2. Notice the clear formation of two distinct clusters for (A) bulb and (B) tibia, but not so clear differences for (C) metatarsus.

Statistical analysis and data visualization.—MANOVA revealed significant differences in mean shape for all structures across localities: bulb (P ,0.001, F_30,160 = 3.48), tibia (P ,0.001, F_30,160 = 3.75), and metatarsus (P ,0.001, F_35,155 = 3.02). Data from RW1 showed differences between localities for bulb (Fig. 4A), tibia (Fig. 4B), and size (Fig. 4D), but no clear differences could be observed with metatarsus RW1 (Fig. 4C). Shape differences corresponding to extremes in RW1 were visualized using thin-plate spline deformation grids (Fig. 5). For the bulb, positive RW1 values represent a long, elongated structure without a flanged tip (Fig. 5A), while negative RW1 values represent a shorter, flanged structure (Fig. 5C). For tibia, positive RW1 values represent a narrower structure (Fig. 5D), and negative RW1 values represent a wider structure (Fig. 5F). For metatarsus, positive RW1 values represent a slim structure with a short metatarsal depression (Fig. 5G), and negative RW1 shows a broad structure with a long metatarsal depression (Fig. 5I), though variation in shape along this axis was not statistically significant.

Figure 4.

Boxplots for quantitative values of RW1 per locality for (A) bulb, (B) tibia, (C) metatarsus, and (D) size (measured as cephalothorax length).

Figure 5.

Thin-plate spline deformation grids showing extreme positive RW values for (A) bulb, (D) tibia, and (G) metatarsus; consensus values for (B) bulb, (E) tibia, and (H) metatarsus; and extreme negative RW values for (C) bulb, (F) tibia, and (I) metatarsus. Consensus deformation grids (B, E, H) show landmark placement at each structure (see Table 1 for descriptions).

Locality comparisons.—Samples from Mt. Magog Low, Mondurup, and Two Mile Lake show positive RW1 values for bulb, negative RW1 values for tibia, and positive RW1 values for metatarsus. Those from Mt. Magog High, Talyuberlup, and Toolbrunup show the opposite, with negative RW1 values for bulb, positive RW1 values for tibia, and negative RW1 values for metatarsus (Fig. 5). Bulb data (Table 3A) suggests there is no difference between samples from Mt. Magog High, Talyuberlup, and Toolbrunup, and shows a strong difference between Mt. Magog Low & Mt. Magog High, Talyuberlup, and Toolbrunup. Two Mile Lake samples differ from Mt. Magog High, Mt. Magog Low, Talyuberlup, and Toolbrunup, while Mondurup did not differ significantly from any other locality.

Table 3.

Pairwise comparisons between localities using RW1 for each structure. Bold values represent significant differences.

Similar results were found with tibia data (Table 3B), with the exception that no difference was found between Mt. Magog Low and Two Mile Lake. Metatarsus showed almost no differences across all localities, with differences only detected between Two Mile Lake samples and those from Talyuberlup and Toolbrunup (Table 3C).

DISCUSSION

Morphotype recognition and distribution.—Based on the two statistically significant traits, bulb and tibia, two distinct morphotypes could be recognized amongst the five localities. Morph1 is characterized by having a long, elongated bulb (Fig. 2A), a wide tibia (Fig. 2F), and was found on individuals from Mt. Magog Low, Mondurup, and Two Mile Lake (Fig. 1). Morph2 has a short, flanged bulb (Fig. 2B), a narrow tibia (Fig. 2E), and was found on individuals from Mt. Magog High, Talyuberlup, and Toolbrunup (Fig. 1).

Biogeographical patterns could not be accurately estimated given the limited sample size across localities. Considering the overlapping distributions of both morphotypes, and the fact that there is no apparent altitudinal separation between them, the morphotypes are expected to show some level of sympatry. This contradicts what has been found for other dispersal-limited invertebrate species within the Stirling Range (Edwards & Harvey 2010; Harvey et al. 2015; Rix et al. 2017). It is possible that even though both morphotypes of Proshermacha can be found on the same mountain, they inhabit isolated areas with little to no overlap between them and are separated by unknown ecological constraints. Further sampling efforts should be conducted to address this topic and provide a better understanding of the morphotype biogeography within the Stirling Range.

Trait effectiveness for morphotype recognition.—Out of the three structures, bulb and tibia showed more evident differences between localities, while metatarsus provided no statistical separation and does not appear to differentiate localities clearly (Fig. 4). Copulatory structures in spiders have been widely used in spider systematics to differentiate between morphologically similar species (Huber 2004). Interestingly, not only do primary copulatory organs evolve rapidly but also secondary copulatory structures such as the tibial spurs are considered useful when differentiating between closely related species (Simmons 2014). Indeed, genital morphology is one of the best traits when it comes to understanding species-level divergence (Eberhard 2010), and the rapid evolution of the copulatory organs is the most likely reason why our results suggest stronger evidence for recognizing two distinct morphotypes when using bulb and tibia (tibial spur). Both structures are known to aid directly in mating by affecting insemination and grasping female fangs/legs/coxae during copula, respectively (Foelix 2011; Pérez-Miles & Perafán 2017). In contrast, little is known about how the metatarsus might influence mating in spiders.

An initial visual examination of samples might suggest there to be three different bulb types among the available samples (Figs. 2A–C). However, our data analysis recognized only two morphotypes (Figs. 2A,B). Through a qualitative approach, the individual in Fig. 2C may appear visually distinct, but under a quantitative analysis of shape, it appears to represent a size variation of Morph1. Further study including molecular data will provide information to determine whether there is a substantial genetic divergence between morphotypes to be considered different valid species. Additional studies including female samples need to be made, and given female genitalia can diverge as rapidly as males (Simmons & Fitzpatrick 2019), it is expected that morphological examination of female samples should reveal a similar pattern of divergence across localities.

Geometric morphometrics as a taxonomical tool.—With further methodological advances, different lines of evidence can be integrated onto species delimitation (Agnarsson & Kunter 2007), combining quantitative data along with descriptive characters. Here we provide an example where GM appears as a useful tool for arachnid taxonomy. In mygalomorphs, contradicting results have been found when contrasting morphological and molecular data. Some studies have distinguished between different morphological and molecular lineages within cryptic species (Wilson et al. 2021), while others have found no morphological differences despite clear molecular divergence (Bond & Stockman 2008). Bond & Stockman (2008) found no morphological differences when using male mating claspers (leg I) in a group of trapdoor spiders, and as extensive as their quantitative measures were, our results suggest that the assessment of other reproductive structures might yet reveal differences between species.

It has been suggested that morphological data have become accessory to molecular phylogenetics when documenting speciation (Bond & Stockman 2008; Satler et al. 2013; Sharma et al. 2021). However, morphological data should remain part of an integrative approach for species descriptions, providing one of many methods contributing to species delimitations and part of the transformation in taxonomy to a collaborative, integrative information science (Bond et al. 2022). An approach like GM can provide statistical support in the recognition of morphological divergence among samples and could become a fundamental statistical approach in the arachnologists' toolkit.

Limitations of Geometric Morphometrics.—Landmark placement and orientation can become complicated tasks when dealing with intricate morphological structures, and it increases in complexity with the addition of semi-landmarks for analyzing outlines or curvature, such as in the palpal bulb. Consistent orientation of the structures is crucial for the landmarks to effectively capture any shape variation, and this can sometimes become a complicated mission when dealing with small structures such as spider genitalia.

When studying closely related populations, a sample size of several dozen individuals has been suggested (Cardini et al. 2015) and limitations on sample size can limit the power of shape analysis. Confidence intervals and variance estimates, as measures of uncertainty, are negatively correlated with sample size, where smaller samples exhibit larger confidence intervals and variance estimates. Not only statistical results, but also the visualization of shape variation might also be affected by small sample size, resulting in estimates of mean shapes not being truly representative of true biological variation (Cardini & Elton 2007). There is a possibility that the “third” bulb type identified from the initial visual examination (Fig. 2C) might represent a distinct morphotype but was not recognized in our statistical analysis because the program failed to fully eliminate size variation from a bimodal distribution of male size due to limited sample size, even after a Generalized Procreustes Analysis was conducted. However, the capacity of our data to differentiate between sympatric individuals suggests that this was not the case and that the specimens matching Fig. 2C morphology are a variation of Morph1. Nonetheless, further sampling should be conducted on the different Stirling Range populations to corroborate these findings.

Evaluating a single plane of view (ventral) from 2D images alone might be a limiting factor in the identification of morphological divergence between populations, and alternate planes of view could provide additional information. Other studies have compared the effectiveness of 3D over 2D images in shape analysis, suggesting that 2D images may underestimate shape variation, resulting in erroneous conclusions (Cardini 2014). This is likely to be especially problematic for complex structures with greater variation (Buser et al. 2018). Nevertheless, this is not always the case, and for some structures a traditional 2D approach can be effective (McWhinnie & Parsons 2019).

In conclusion, even though GM can provide quantitative evidence for improving taxonomic studies it should not be implemented alone, but rather used as a template to which other morphological datasets such as MRI (Ziegler et al. 2011) or lCT scans (Semple et al. 2018) can be added. Our findings demonstrate that 2D GM can identify morphological divergence between sympatric populations, and as others have shown (Crews & Hedin 2006; Franzini et al. 2013; Wilson et al. 2021) they can be paired with molecular methods to resolve species hypotheses with quantitative morphological support.