Camel spider trait evolution demonstrates repeated patterns of convergence (Arachnida: Solifugae: Eremobatidae)

Abstract

Morphology has long been used to classify and identify living organisms. However, taxonomic descriptions are often limited to qualitative descriptions of size and shape, making identification difficult due to the subjective language used to describe complex shapes. Additionally, for some taxa, there are few reliable qualitative characters available for delimitation that have yet to be tested objectively in a phylogenetic context. Solifugae is one such example. The order, Solifugae, is recognized from the other arachnid orders by the possession of large, powerful jaws or chelicerae. Male cheliceral morphology is the leading diagnostic character system in solifuge systematics and is the basis for much of solifuge current taxonomy. Female chelicerae, on the other hand, are reportedly deeply conserved and much of the species identification is based on female operculum morphology. To elucidate patterns of chelicerae and opercula trait evolution within the solifuge family, Eremobatidae, we used a 2-dimenstional morphological analysis using an Elliptical Fourier approach for closed outlines, in addition to an analysis of traditionally used measures in a phylogenetic context. Using ancestral state reconstruction and ultra-conserved elements, we assessed the taxonomic utility of female cheliceral and opercular morphology, and we evaluated which male morphological characters reflect shared, derived ancestry. Investigation into ubiquitously used character sets, in addition to newly proposed characters herein, illustrates the complex evolution of traits with high levels of convergence. Our results provide taxonomic insight into future, higher level taxonomic revisions of Eremobatidae.

Solifuge, chelicerae, Elliptical Fourier analysis, Gaussian mixture modeling, morphometrics

Introduction

The primary goals of phylogenetics and taxonomy are to describe, document, and understand the descendance of both extant and extinct groups of life. These 2 fields are complementary: descent with modification is a position of central importance for modern taxonomy. Utilizing a collection of evidence from multiple independent sources is integral to achieving the goals of phylogenetics and taxonomy. One such piece of information that has been long used to classify and identify orders of life is morphology. Together with phylogenetics, analyzing morphological traits can not only help to classify life, it may also open researchers to evolutionary-based questions regarding the mechanisms behind functional morphology (Wood et al. 2016, Wood 2020), the relationship between morphology and feeding strategies (Waloszek et al. 2007, Johnson et al. 2018), or the emergence of convergent evolution with adaptive radiations (Losos and Ricklefs 2009, Gillespie et al. 2018, Cerca et al. 2023). However, analyzing the evolution of morphological traits is incredibly difficult for some organisms as their body plan has been deeply conserved for millions of years (Cognato and Grimaldi 2009, Bryson et al. 2018) and the number of informative characteristics for delimitation is sparse. Another taxonomic challenge is that species descriptions are often restricted to qualitative characterizations of morphological size and shape. Shapes of morphological structures are often too complex to be captured with descriptive language (Bonhomme et al. 2014). For example, subtle descriptive differences found in taxonomic keys, such as fondal notch broadly C-shaped vs. J-shaped vs. irregularly U-shaped (Muma 1962) could make comparisons of 2 closely related species challenging for a confident species identification. Moreover, limiting shape descriptions to a single descriptor excludes the consideration for variation that may exist intraspecifically. Thus, to robustly examine variation within taxa, morphological variation must be considered both qualitatively and quantitatively. Assessing the extent of morphological variability of prominent characters within a species can facilitate future species designations and, in the presence of an evolutionary framework, can provide insight into the evolutionary history of commonly utilized characters for taxonomic revision.

To investigate the morphological evolution of shape through objective classification, we used solifuges as a model system. Solifuges, commonly known as camel spiders, are an understudied, diverse group of mostly nocturnal animals. Despite their recognized diversity, solifuge taxonomy is incredibly challenging due to the limited diagnostic characters available for confident species identification and taxonomic delimitation. Additionally, much of the available resources for taxonomic identification are nearly between 60 and 90 yr old for some solifuge families and are based on a small set of highly variable, and widely controversial characters (Muma 1951, Lawrence 1955, Maury 1982, Gromov 2000). The solifuge family, Eremobatidae Kraepelin 1899, is endemic to western North America and as per current taxonomy based on morphology, this family consists of 183 species, 2 subfamilies, 8 genera, and 18 species groups (Harvey 2003) making them one of the most diverse solifuge families. However, this diversity is reflective of a small subset of characteristics. In 1951, Muma revised Eremobatidae primarily based on male cheliceral morphology and associated setae. As of now, the gross male cheliceral morphology continues to be the basis of eremobatid taxonomy, due to the presence of conspicuously variable, diagnostic characters. Given the number of characters available for delimitation and the reported diversity, eremobatids are an ideal model system for evaluating taxonomic boundaries based on objective classification.

Solifuges are recognized from other arachnid groups for their large, powerful chelicerae. Chelicerae functionally serve both males and females for feeding (Turner 1916, Muma 1966b, 1967, Cloudsley-Thompson 1977), digging and burrowing (Cloudsley-Thompson 1977), or antagonistic interactions (Muma 1967). Despite some overlapping utility in both sexes, chelicerae are sexually dimorphic. Male chelicerae vary in dentition, shape, seta types, and shape of the mesoventral groove across multiple taxonomic levels (Fig. 1; Muma 1951, 1970, Cushing et al. 2018). Variability in male chelicerae can be attributed to their role in reproduction. Camel spiders, like many other arachnids, lack a direct intromittent organ and instead rely on a secondary morphological structure, specifically their chelicerae, for sperm transfer and for positioning the female during copulation (Muma 1966a, Rowsell and Cushing 2020, Peretti et al. 2021). Species-specific variation among males is hypothesized to be attributed to strong sexual selective pressure, like other sexually dimorphic characters such as body size (Fox et al. 1995), coloration (Punzalan et al. 2010), or ornamentation (Schilthuizen 2003) found in other arthropods that help to promote reproductive success. Female solifuge cheliceral morphology, on the other hand, is reportedly conserved as deep as the family level (Bird et al. 2015) and has historically been suspected to be conserved within Eremobatidae. However, female eremobatid cheliceral morphology (Fig. 2) has yet to be thoroughly characterized and analyzed in an evolutionary context. Female genital plate morphology (otherwise known as operculum morphology), on the other hand, has been demonstrated to be useful in species identification, due to diagnostic variation between species (Muma 1951). In 2015, Cushing et al. published the first molecular phylogenetic hypothesis for Eremobatidae. Notably, the findings of Cushing et al. (2015) demonstrated that Eremobatidae needs extensive taxonomic revision since the recognized delimiting morphological characters recover many of the genera either poly- or paraphyletic. Thus, a comprehensive phylogenetic framework for the taxa that comprise the group is necessary for studying and evaluating character systems that were used to establish the current taxonomic relationships.

Fig. 1.

Male chelicera variation from the family Eremobatidae. A) Hemerotrecha parva (DMNS ZA.42131), B) Hemerotrecha macra (DMNS ZA.35652), C) Eremochelis branchi (DMNS ZA.38315), D) Eremochelis flexacus (DMNS ZA.16134), E) Chanbria regalis (DMNS ZA.25444), F) Eremorhax magnus (DMNS ZA.42055. Photograph credit: Cole Logan), G) Hemerotrecha hanfordana (DMNS ZA.37358), H) Horribates spinigerus (DMNS ZA.40086), I) Eremobates leechi (DMNS ZA.37070. Photograph credit: Quincy Hansen). All scale bars refer to 1 mm, except F, which refers to 2.5 mm.

Fig. 2.

Female cheliceral morphology from representatives of the family Eremobatidae. A) Hemerotrecha delicatula (DMNS ZA.41812), B) Eremochelis striodorsalis (DMNS ZA.18968), C) Horribates spinigerus (AMNH Holotype), D) Chanbria rectus (DMNS ZA.25456), E) Eremochelis andreasana (DMNS ZA.40822), F) Hemerotrecha serrata (DMNS ZA.25449), G) Hemerotrecha sp. (CIDA107821), H) Hemerotrecha prenticei (DMNS ZA.18169), I) Eremochelis bilobatus (DMNS ZA.17685). All scale bars refer to 1 mm.

As an exploratory study on the evolution of commonly utilized character systems in eremobatids, the aims of the present study were 2-fold: to establish a comprehensive evolutionary hypothesis and to quantify morphological variability in eremobatid taxa to ultimately assess the phylogenetic integrity of such morphological traits to aid in future taxonomic revisions. First, we sought to establish a new phylogenetic hypothesis from phylogenomic data with a focus on 2 of the most diverse eremobatid genera: Eremochelis and Hemerotrecha. As many of the deeper-level relationships were weakly supported in Cushing et al. (2015), we sought to bring more resolution to the evolutionary history of eremobatids by including taxa that were previously underrepresented. The morphological focus of this contribution was an attempt to understand the level of morphological variation among the focal taxa and to objectively characterize morphological traits. Specifically, we sought to independently analyze morphological traits to discover consistent morphological patterns. As an extensive examination of the selected morphological characters, we demonstrate the reliability of prevalent morphological traits for future taxonomic utility, highlight important patterns of trait evolution, and thus provide a foundation for a monographic revision of several genera within Eremobatidae.

Materials and Methods

Taxon Sampling

Taxon sampling for phylogenetic estimation included 196 eremobatid individuals representing 6 out of 8 eremobatid genera: Hemerotrecha, Eremochelis, Chanbria, Horribates, Eremobates, and Eremorhax. We focused rigorous sampling efforts on the 2 of the most speciose, yet historically neglected genera: Eremochelis and Hemerotrecha (Table 1). For outgroups, we selected four individuals within 3 outgroups families: Solpugidae, Daesidae, and Ammotrechidae.

Table 1.

Taxonomic breakdown of eremobatid taxa used in this study partitioned by genus, species group, the fraction of taxonomic coverage, and unidentified/putative new species. The fraction of taxonomic coverage indicates the number of species out of the total number represented by that group that we used as exemplars in our study. Not all genera have species groups described. Detailed taxonomic breakdowns for all terminals can be found in Supplementary Information 1

Genus	Species group	Taxonomic coverage	Unidentified species/putative new species
ChanbriaMuma 1951	N/A	3/4 Species
EremorhaxRoewer 1934	N/A	1 Species
HorribatesMuma 1962	N/A	1/3 Species
Putative new genus	N/A	1 Species	1 Putative new genus
Eremochelis Roewer 1934		25/37 Species
	andreasana	1/2 Species
	bilobatus	12/16 Species	2 Putative new species
	branchi	6/14 Species	3 Putative new species
	imperialis	4/5 Species	2 Putative new species
	striodorsalis	1/1 Species
Hemerotrecha Banks 1903		23/34 Species
	banksi	7/9 Species	1 Putative new species
	branchi	7/9 Species	1 Putative new species
	denticulata	4/6 Species	3 Putative new species
	serrata	1/1 Species
	simplex	2/7 Species	3 Putative new species
	texana	1/1 Species
Eremobates Banks 1900		5/77 Species
	palpisetulosus	3/41 Species
	pallipes	0/17 Species	1 Unidentified species
	aztecus	1/1 Species
	No species group placement	1 Species

We examined specimens from the following institutions: Denver Museum of Nature and Science (DMNS), California Academy of Sciences (CAS), American Museum of Nautral History (AMNH), University of Arizona (UA), College of Idaho Orma Smith Natural History Collection (CIDA), Colección Arachnológia del Centro de Investigación Biologicas Del Noroeste, S.C. (CARBIO), Essig Museum of the University of California Berkeley (ESS), Colección Naciónal Arácnidos Instituto de Biología de la Universidad Naciónal Autónoma de México (CNAN), the San Diego Natural History Museum (SDNHM), and the Florida State Collection of Arthropods (FSCA). For morphometrics analyses, when available, we examined multiple male, and female representatives from the 6 focal eremobatid genera to help capture intraspecific variation. Images of retrolateral cheliceral views and ventral operculum views were taken using a Leica M125C Stereo Microscope and an Olympus SZX12 Microscope. Individual images were rendered using Leica Application Suite X (LAS X) software and Helicon Focus to produce 1 composite image for each morphological view. Traditional measures (Fig. 3A) were either measured directly in the LAS X software, or if images taken using the Olympus SZX12, structures were subsequently measured in Fiji (Schindelin et al. 2012). To characterize qualitative male traits, we visualized and imaged prolateral views of male chelicerae using a Hitachi TM4000Plus II Tabletop Scanning Electron Microscope.

Fig. 3.

Summary of the quantitative and qualitative traits considered in this study for eremobatid males and females. A) Right, retrolateral view of male Eremochelis insignatus (DMNS ZA.33749) chelicera. Measurements considered in this study are illustrated. B) Scanning electron microscope (SEM) image of papillated texturing inside flagellar groove of Eremochelis branchi (DMNS ZA.37135). C) SEM image of prolateral view of male fixed finger and associated flagellar groove belonging to Eremochelis branchi (DMNS ZA.37135). D) Right, retrolateral view of female chelicera of Hemerotrecha delicatula (DMNS ZA.41812) depicting the fixed finger tooth pattern and tip (apex) of fixed finger to fixed finger distal tooth (FF-to-FD) measure. Coloration and abbreviations follow Bird et al. (2015). Terminology abbreviations are FD = fixed finger, distal teeth, FSD = fixed finger, subdistal teeth, FM = fixed finger, medial teeth, FSM = fixed finger, submedial teeth, FP = fixed finger, proximal teeth. E) Operculum length, width, and diagonal measures on female Eremochelis plicatus (DMNS ZA.41887).

The taxonomic breakdown for each species considered for traditional morphometrics, EFA, and UCE analysis can be found Supplementary Information 1, 2, and 3 respectively.

DNA Extraction, Library Preparation, and Bioinformatics

Genomic DNA was extracted from fresh leg or cheliceral tissues using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA) following the manufacturer’s protocol. For poorly preserved specimens and specimens stored at room temperature, we extracted genomic DNA following the protocol for degraded tissues outlined in Tin et al. (2014). We used a Qubit Flourometer 4.0 (Life Technologies, Inc.) to determine DNA concentration of each extracted sample. To determine whether extracts required fragmentation for downstream library preparation, a necessary step for next-generation sequencing, we ran select extracts on a gel. We ran select extracts to avoid exhausting DNA volumes necessary for downstream library preparation, especially for those low DNA input samples. DNA isolates with high molecular weight were fragmented using the KAPA HyperPlus Kit at half the reaction of the manufacturer’s protocol.

Genomic libraries, target DNA that is processed to be compatible with next-sequencing technologies, were prepared at the DMNS. Due to COVID-19 restrictions, some genomic libraries were prepared and sequenced with Daicel Arbor Biosciences (Ann Arbor, MI, USA). For libraries prepared at DMNS, we used the KAPA Hyper Prep Kit (Kapa Biosystems) following Faircloth et al. (2015) but using a modified protocol from Derkarabetian et al. (2019) for use with the iTru dual-indexing adapter system (Glenn et al. 2019). Libraries were enriched via PCR in a reaction summarized in Derkarabetian et al. 2019 using a unique Illumina dual-indexed (i5 and i7; Glenn et al. 2019) motif for each sample. To ensure amplification success, amplified libraries were quantified using a Qubit 4 fluorometer. Enriched libraries were strategically pooled based on the enriched library yield. Target enrichment of libraries, by way of hybrid capture, was achieved via the MYbaits Arachnids 1.1K version1 kit (Daicel Arbor Biosciences) following the myBaits Hybridization Capture for Targeted NGS protocol provided by Arbor Biosciences (https://arborbiosci.com/wp-content/uploads/2018/04/myBaits-Manual-v4.pdf), with adjustments following Derkarabetian et al. (2019) for museum samples. Posthybridization, all pooled libraries were quantified via Qubit to verify that the hybridization step was fulfilled. Hybrid enriched libraries were multiplexed and sequenced on a partial lane using the Illumina NovaSeq 6000 at Novogene (University of California Davis, CA).

Raw, demultiplexed reads were processed primarily using the PHYLUCE pipeline (Faircloth 2016), with additional software to complement data processing. First, adapter sequences were trimmed from raw sequencing data via Illumiprocessor, (Faircloth 2013). Following adapter trimming, to identify other undesired sequences such as remnant Illumina primers or artificial Poly G tailing, a familiar issue for Illumina HiSeq Series machines that rely on a two-color chemistry system (Chen et al. 2018), we generated sequence summary reports for each sample using FastQC (Andrews 2010). Summary reports include (but are not limited to) per base sequence quality, per sequence quality scores, per base N content, and an indication of overrepresented sequences. From the overrepresented sequence hits, we generated a single FASTA file containing specific contaminant sequences that matched Illumina primer and sequences containing more than 10 guanine repeats. We removed contaminant sequences from adapter-cleaned sequences using fastp (Chen et al. 2018) and the used single FASTA file as input.

Clean reads were assembled de novo, using metaSPAdes assembler (Nurk et al. 2017) and default parameters. UCE loci were identified from resultant contigs using the arachnid UCE sequence probes and the match_contigs_to_probes function within PHYLUCE. Recovered UCE loci were aligned using MAFFT (Katoh and Standley 2013) and were trimmed internally using TrimAI (Capella-Gutiérrez et al. 2009) with a gap-threshold value of 0.2 as recommended by Portik and Wiens (2021). To minimize the effect of missing data, taxa that rendered less than 10 UCE loci were not considered for downstream analyses. Additionally, older museum samples that produced long branches in preliminary phylogenetic analyses were not considered in the final analyses. The final data matrix for this study was generated using alignments with ≥75% locus occupancy attained using the function phyluce_align_get_only_loci_with_min_taxa. And finally, resultant alignments were concatenated using the phyluce_align_concantenate_alignments function.

Phylogenetic Estimation

Maximum likelihood (ML) analyses utilizing the concatenated 75% locus occupancy matrix was performed in IQ-TREE 2 (Minh et al. 2020) with ultrafast bootstrap approximation and best-fit models of evolution were identified using ModelFinder (Kalyaanamoorthy et al. 2017). For coalescent-based phylogenetic estimation, individual gene trees using the internally trimmed alignments were also estimated in IQ-TREE 2 with ModelFinder. Inferred gene trees were then used as input in ASTRAL-III v 5.7.8 (Zhang et al. 2018), and nodal support was summarized by local posterior probabilities approximated by gene tree quartet frequencies (Sayyari and Mirarab 2016).

Elliptical Fourier Analysis for Closed Outlines

For 2D geometric morphometric analysis of male and female chelicerae and opercula, we implemented an Elliptical Fourier analysis (EFA) for closed outlines (Kuhl and Giardina 1982) using the Momocs package in R (Bonhomme et al. 2014, R Core Team 2022). The EFA mathematically characterizes complex 2D shapes from sine and cosine functions. In this method, contours in an x and y coordinate plane are described by a series of related ellipses, called harmonics, that can decompose an object outline through additive periodic functions (Kuhl and Giardina 1982, Bonhomme et al. 2014, Caple et al. 2017). Using EFA, geometric shapes are quantified and can be analyzed using multivariate approaches (Bonhomme et al. 2014). This method is favorable over other traditional methods, such as the analysis of linear measurements and ratios because any outline can be described using an appropriate number of harmonics; moreover, equally spaced points are not necessary for the analysis (Rohlf and Archie 1984, Crampton 1995, Renaud and Michaux 2003). Additionally, in contrast to non-Fourier methods, EFA eliminates the need for establishing homologous landmarks across disparate taxa.

All structures were converted to a black silhouette on a white background, saved as a.jpg file, and used as input in Momocs. Binary images, otherwise known as silhouette images, of each generated composite image were rendered using GIMP 2.10 (The GIMP Development Team 2019) and Adobe Photoshop Software (Adobe Inc. 2021). Due to the asymmetry of female opercula and discontinuity of the plates, only the left plate was considered for downstream shape analysis. Following import, image outlines were smoothed using 200 smoothing iterations to eliminate pixel artifacts using the coo_smooth function and interpolated using 1,000 points using the coo_interpolate function in the Momocs package. The taxonomic breakdown and associated unique identifiers of the outlines considered in this study can be found in Supporting Information 2..

The size of a structure, in conjunction with overall shape, may be an important evolutionary consideration. Therefore, for each image, we rescaled each outline based on a pixel/mm ratio for each respective image using the coo_scale function in Momocs. All chelicera outlines were of the retrolateral views, and all outlines in this study were scaled to reflect actual sizes. Following rescaling, all centroids were aligned and centered prior to estimating the number of harmonics needed to conduct an EFA. The optimal number of harmonics, estimated by capturing 99.9% of variation in the dataset, was inferred using the calibrate_harmonicpower_efourier function. Once the ideal number of harmonics was determined, we conducted an EFA using the efourier function and the resulting Fourier coefficients from the optimal number of harmonics was used downstream in principal components analyses (PCA). The first 3 principal components (PCs) were retained for downstream model-based clustering (see Morphospace clustering section). Shape variation in morphospace was summarized using the first 2 PCs and visualized using the plot function in R as these 2 PCs summarized >50% of the variation (see Results). We also performed an allometric test on the male dataset to determine whether there was a relationship between size and centroid shape as in Kallal et al. (2022). We calculated scaled centroid size using the coo_centsize found in the Momocs package and used the first 3PCs as dependent variables for quantifying the shape of the uniformly scaled outlines. We then performed a linear test using log-transformed centroid size as the predictor using the lm.rrpp function in the RRPP package (Collyer and Adams 2018). We decided to use our male dataset only since this consisted of the greatest sampling.

Measures for Traditional Morphometrics

We performed traditional morphometrics on common cheliceral measurements as found in eremobatid literature, those being cheliceral length (CL), height (CH), fixed finger height (FFH), and aspect ratios (Muma 1951, 1963, 1971, Brookhart and Muma 1981, Brookhart and Cushing 2002, 2004, 2008, Ballesteros and Francke 2008, Cushing et al. 2018). The purpose of this aspect of this study was to evaluate the reliability of such commonly used characters for delimitation of taxonomic boundaries. Additionally, we sought to investigate these commonly reported measurements to inform the extent of morphological variation that may be present within and between taxa and to use as continuous character for objective classification via model-based clustering analyses.

Aspect ratios, introduced to solifuge taxonomy with the intention of describing relative shape, generally concern length-to-width or length-to-height comparisons. For example, the “jaw index,” a length/height (CL/CH) ratio proposed by Cloudsley-Thompson (1961) and synonymous with the cheliceral length/width (CH/CW) ratio introduced by Brookhart and Muma (1981), is used to characterize cheliceral dimensions (Bird et al. 2015). This ratio is meant to capture relative shape, such that large CL/CH ratios suggest a narrow cheliceral morphology, while a more robust morphology is represented by a smaller ratio. Other ratios have been proposed for taxonomic utility, yet only 1 other index is inclusive of all eremobatid species. The fondal notch (FN), defined here as the distinct indentation extending beyond the fondal axis and located between the fixed finger (FF) and fondal teeth, is absent in some species of Eremobatidae (Supplementary Fig. S3). The cheliceral height/fixed finger height (CH/FFH) ratio was introduced to eremobatid taxonomy to provide measurements for those species for which the FN was lacking (Brookhart and Cushing 2002); thus, we adopted the usage of this inclusive aspect ratio for all representatives used in this study. To investigate the utility of the female genital operculum length-to-genital operculum width ratio (GOL/GOW), a ratio first introduced by Brookhart and Cushing (2002), we recorded GOL, GOW, and genital operculum diagonal (GOD)—determined from measuring approximately 45° from the length starting point to the base outline of a single plate (Fig. 3E). We also considered novel measurements first reported here and visualized in Fig. 3A and C, those including tip of FF-to-FF distal tooth (FD) for female chelicera, tip of FF to the nearest point on the FF (used to quantify curvature in males), tip of the moveable finger (MF) to MF proximal tooth (MD), and tip of MF-to-MF medial tooth (MM).

From the measures recorded, we performed 3 independent PCAs for male cheliceral measures, female cheliceral measures, and female operculum. Final morphometric matrices for a male chelicera included 6 continuous characters: cheliceral length I (CLI), CH, FFH, tip of MF-to-MD, and tip of MF-to-MP, tip of FF to nearest point of FF (Fig. 3A). Matrices for a female chelicera included length I, CH, FFH, tip of MF-to-MD, tip of MF-to-MM, and tip of FF-to-FD (Fig. 3A and D). Female operculum measures included GOL, GOW, and GOD (Fig. 3E).

For both male and female chelicerae, few measures were not taken due to specimen damage, thus final data matrices contained some missing measurements for the taxa considered in this study. Most PCAs, however, can only be applied to complete data sets, therefore, to resolve the issue of missing data, we estimated missing values using the missMDA package in R (Josse and Husson 2016). This package permits the use of PC methods for an incomplete data set by predicting missing values based on parameter estimation, considering similarities between individual measures and the relationship between the categorical variables of interest (Josse and Husson 2016). Missing values are estimated using multiple imputation procedure, then iterative PCA algorithm or regularized iterative PCA algorithm (Josse and Husson 2016). We found this method suitable as some individuals had a near complete dataset, and therefore, we sought to keep these individuals to maintain adequate taxon sampling for some species representatives. To estimate missing values, we first implemented the estim_ncpPCA function under a K-fold cross-validation to infer the number dimensions to be used for predicting missing values. Next, we estimated missing values using the regularized iterative PCA algorithm in the imputePCA function, as this method is best to avoid overfitting when many values are missing (Josse and Husson 2016). After predicting missing values in our data sets, we performed a PCA using the prcomp function in R. Morphological variation for male chelicerae, female chelicerae, and female opercula, in morphospace was summarized by genus and visualized using the first 2 PCs. Similarly, the first 3PCs were used for downstream model-based clustering analyses as they explained most of the variation in each data set (see Results).

In this study, we were motivated to evaluate some of the prominent qualitative male characteristics used to establish eremobatid generic boundaries. The purpose of investigating such qualitative characteristics was to assess their utility in supporting new taxonomic boundaries based on common descent. In total, we scored 7 qualitative morphological characters for male representatives included in this study and those characters and character descriptions are summarized in Table 2. We were also interested in investigating the taxonomic utility of female cheliceral morphology as this has yet to be examined; thus, we scored 2 female qualitative character states (Table 2). First, we recorded female tooth formula as described in Bird et al. (2015) from the FF subdistal teeth (FST) up until the FF proximal tooth (FP; Fig. 3F). Next, we recorded female tooth count from the FD to the FP. Scored character states used for analysis can be found in Supporting Information 1.

Table 2.

Summary of the qualitative characteristics used in this study and the descriptions of each character. Characters used for each sex are indicated

Character	Character description of qualitative characters
Flagellar groove (FG) ♂	The flagellar groove is a concave structure located on the prolateral side of the fixed finger. See Fig. 3C
Fondal notch (FN) ♂	The fondal notch is located ventrally with respect to the fixed finger when viewing the chelicera laterally. It is an indentation extending significantly past the fondal axis. The fondal axis is a reference line that includes the linear base of the fondal teeth. Any indentation beyond this axis is considered a notch. See Supplementary Fig. 3.
Relative flagellar complex length (FCL) ♂	The flagellar complex length here only considers the length of the basal area of the complex with respect to the fixed finger, not the length of the setae themselves. See Supplementary Fig. 1.
Papillated texture ♂	Papillated texture are small, rounded protuberances located on the fixed finger, often within the flagellar groove. See Fig. 3B.
Papillated texture position ♂	Papillated texture, defined above, is variable in position with respect to the fixed finger. Texture fields are generally limited to proximal or distal positions, or in some cases, are distributed throughout the flagellar groove.
Medial tooth, moveable finger (MM) ♂	The medial tooth, moveable finger is a large anterior tooth located on the fixed finger, yet smaller than the proximal tooth. See Supplementary Fig. 2.
Central carina (CC) ♂	A keel-shaped structure, or ridge, located within the flagellar groove. See Fig. 4 and Supplementary Fig. S18 for ASR.
Tooth count ♀	Tooth count quantifies the number of teeth starting from the fixed finger, distal tooth to the fixed finger, proximal tooth. See Fig. 3D for tooth terminology.
Tooth formula ♀	The tooth formula is in accordance with the terminology established in Bird et al. (2015). See Fig. 3D and figure legend for tooth terminology. Any preceding numbers before the terminology abbreviation indicate the number of teeth for that tooth characterization.

Table 3.

Characteristics and GMM clustering analysis of the continuous variables quantified in this study. The best-fit geometric models, number of clusters, and number of individuals per cluster before and after uncertainty thresholds are summarized. The best-fit models for each data input are as follows: VEE = ellipsoidal, equal shape and orientation, VEI = diagonal, equal shape, VVE = ellipsoidal, equal orientation

Data set variable	Best-fit model	Number of clusters	Number of individuals per cluster membership	Number of individuals per cluster membership (0.25 uncertainty cut off)	Number of individuals per cluster membership (0.05 uncertainty cut off)
Male chelicera measures	VEE	7	Cluster 1: 26 Cluster 2: 15 Cluster 3: 56 Cluster 4: 51 Cluster 5: 14 Cluster 6: 24 Cluster 7: 74	Cluster 1: 24 Cluster 2: 14 Cluster 3: 48 Cluster 4: 41 Cluster 5: 11 Cluster 6: 20 Cluster 7: 63	Cluster 1: 24 Cluster 2: 11 Cluster 3: 35 Cluster 4: 23 Cluster 5: 3 Cluster 6: 12 Cluster 7: 24
Female chelicera measures	VEE	3	Cluster 1: 62 Cluster 2: 4 Cluster 3: 29	Cluster 1: 60 Cluster 2: 4 Cluster 3: 25	Cluster 1: 53 Cluster 2: 4 Cluster 3: 14
Female opercula measures	VEI	3	Cluster 1: 11 Cluster 2: 57 Cluster 3: 33	Cluster 1: 11 Cluster 2: 55 Cluster 3: 31	Cluster 1: 10 Cluster 2: 51 Cluster 3: 22
Male chelicera outlines	VEI	4	Cluster 1: 39 Cluster 2: 11 Cluster 3: 102 Cluster 4: 85	Cluster 1: 34 Cluster 2: 10 Cluster 3: 73 Cluster 4: 73	Cluster 1: 0 Cluster 2: 7 Cluster 3: 25 Cluster 4: 49
Female chelicera outlines	VEI	3	Cluster 1: 25 Cluster 2: 25 Cluster 3: 36	Cluster 1: 24 Cluster 2: 23 Cluster 3: 32	Cluster 1: 20 Cluster 2: 11 Cluster 3: 22
Female opercula outlines	VVE	2	Cluster 1: 86 Cluster 2: 13	Cluster 1: 86 Cluster 2: 13	Cluster 1: 82 Cluster 2: 10

Morphospace Clustering

Model-based clustering using morphometric data was performed in R using the mclust package (Fraley and Raftery 2002) to determine the inherent structure of the morphological data in morphospace with the goal of using any detectable structure to inform taxonomic boundaries within Eremobatidae. To determine the optimal number of clusters from the input data, we executed the Mclust function, which implements both the Expectation-Maximization (EM) algorithm and the Bayesian Information Criterion (BIC) to determine the best-fit clustering model. Under this approach, large BIC values provide support for a specific model, in addition to the optimal number of clusters (Fraley and Raftery 1999). We inferred cluster membership of each individual outline using an uncertainty threshold of 0.05 and 0.25 since uncertainty values range from 0 (highly certain) to 1 (highly uncertain) as in Kallal et al. (2022). Mean shapes of resulting clusters were determined by the function MSHAPES found in the Momocs package.

Ancestral State Reconstruction

One of the common assumptions of estimating ancestral states is that phenotypic change is proportional to the amount of time that has elapsed along a branch (Cusimano and Renner 2014). Contrary to traditional approaches, some studies on character evolution have suggested that traits and the rate of molecular change can be correlated (Davies and Savolainen 2006, Smith and Donoghue 2008), therefore ancestral states inferred from a chronogram, where branch lengths reflect changes through time, can lead to inaccurate ancestral states. Recent studies, however, have explored the effect of branch length on estimating ancestral states. Of the handful of studies that have investigated the sensitivity of branch length choice, it is suggested that researchers estimating ancestral states from continuous traits should opt for the branch lengths that return the highest phylogenetic signal as it increases the accuracy of ancestral state estimation (Litsios and Salamin 2012). For ancestral states modeled for discrete states, on the other hand, Wilson et al. (2022) suggest that the most accurate ASRs are estimated from branch lengths that renders the best model fit as determined by the lowest Akaike information criterion score corrected for small sample sizes (AICc; Burnham and Anderson 2002). Therefore, to determine whether to use a chronogram versus a phylogram for inferring ancestral states, we estimated the phylogenetic signal for continuous characters and compared AICc model statistics (see below) for the best-fit models for discrete traits using our pruned chronograms and phylograms. Results of this comparative evaluation are summarized in Table 4 and Supplementary Table S1 for continuous traits, and Table 5 for discrete states. The exact P-values for Table 4 are uploaded on Figshare. Using the concatenated topology estimated from IQ-tree, which was congruent with the topology rendered from ASTRAL for all major clades (Supplementary Fig. S5), we pruned tips to represent single species hypotheses using the function drop.tip function in the ape package (Paradis and Schliep 2019). Because an updated time-calibrated phylogeny is not available for the increased taxon sampling presented here, we scaled branch length by evolutionary rate using the chronopl package in the ape package (Paradis and Schliep 2019) with an intermediate smoothing rate value (lambda) of 0.1 as it is suggested to work best for topologies with pruned branches (Sanderson 2002).

Table 4.

Summary of eremobatid male and female continuous characters quantified in this study and associated phylogenetic signal statistics calculated using our estimated chronogram. Blomberg’s K value with asterisk refer to a recovered P-value < 0.001 from associated test. Pagel’s λ, log-likelihood, likelihood ratio rest are also reported. Pagel’s λ values with asterisks also indicates an associated P-value < 0.001. Specific P-values are uploaded onto Figshare

Continuous character	Blomberg’s K	Pagel’s λ	Pagel’s λ log-likelihood	Pagel’s λ likelihood ratio test
EFA Chelicera PC1 ♂	1.791*	1.232*	−19.906	51.314
Chelicera length (CL) ♀	1.563*	1.085*	−59.241	27.727
Chelicera height (CH) ♀	1.415*	1.069*	−30.209	22.866
Chelicera length (CL) ♂	1.401*	0.979*	−77.124	33.535
Genital operculum width (GOW) ♀	1.325*	1.286*	−8.569	19.386
EFA Chelicera PC1 ♀	1.308*	1.011*	−35.966	21.206
Chelicera fixed finger height (FFH) ♀	1.235*	1.014*	2.279	17.292
Chelicera height (CH) ♂	1.186*	0.847*	−39.404	25.432
Genital operculum diagonal (GOD) ♀	1.156*	1.267*	−17.776	13.367
Chelicera fixed finger height (FFH) ♂	1.139*	0.901*	32.042	22.233
Genital operculum length (GOL) ♀	1.123*	1.262*	−13.494	12.213
Tip of Fixed Finger to Fixed Finger, Distal Tooth (FF-to-FD)	1.119*	0.897*	7.763	13.765
Chelicera length to chelicera height (CL-to-CH) ♀	1.037*	1.156*	6.763	11.086
Chelicera length to FFH (CL-to-FFH) ♂	1.003*	1.247*	−10.337	20.355
Chelicera length to chelicera height (CL-to-CH) ♂	0.976*	1.255*	3.905	22.776
Genital operculum length-to-width (GOL-to-GOW) ♀	0.935	0.932	−64.475	5.148
EFA Opercula PC1♀	0.922	0.882	1.189	4.125
Chelicera length to FFH (CL-to-FFH) ♀	0.85	0.927	18.38	4.126

Table 5.

Summary of eremobatid male and female discrete characters quantified in this study. Comparative analyses for ancestral state reconstructions are summarized by log-likelihood. Asterisks indicate the best-fit model as determined by the highest LnL score and lowest AICc score from each of the 3 models and topology considered. δ statistic and D-statistic (binary characters only) for measuring phylogenetic signal for discrete characters and associated P-value

Discrete character by sex	Coding	Model: Log-likelihood (LnL), AICc, number of estimated parameters (Chronogram)	Model: Log-likelihood (LnL), AICc, number of estimated parameters (Phylogram)	δ statistic (P-value under a random model) for best model	D-statistic (probability of Brownian phylogenetic structure)
EFA Opercula cluster membership ♀	2 states	*ER: −13.261, 28.647, 1 SYM: −13.261, 28.647, 1 ARD: −12.519, 29.425, 2	ER: −13.261, 28.647, 1 SYM: −13.261, 28.647, 1 ARD: −12.519, 29.425, 2	22.086 (0.9)	−4.152 (0.938)
Flagellar groove (FG) ♂	2 states	*ARD: −37.521, 79.267, 2 ER: −38.816, 79.707, 1 SYM: −38.816, 79.707, 1	ARD: −37.521, 79.267, 2 ER: −38.812, 79.699, 1 SYM: −38.812, 79.699, 1	0 (0.005)	−0.611 (0.734)
Fondal notch (FN) ♂	2 states	*ER: −37.763, 77.601, 1 *SYM: −37.763, 77.601, 1 ARD: −36.996, 78.219, 2	ER: −37.643, 77.36, 1 SYM: −37.643, 77.36, 1 ARD: −37.048, 78.322, 2	2.216 (0.075)	−0.819 (0.796)
Central carina in FG ♂	2 states	*ER: −14.924, 31.922, 1 *SYM: −14.924, 31.922, 1 ARD: −14.41, 33.046, 2	ER: −14.51, 31.094, 1 SYM: −14.51, 31.094, 1 ARD: −14.41, 33.046, 2	66.928 (0.905)	−2.66 (0.861)
Moveable finger, medial tooth (MM) ♂	3 states	*ER: −61.257, 124.588, 1 SYM: −60.15, 126.763, 3 ARD: −58.919, 131.552, 6	ER: −61.522, 125.119, 1 SYM: −61.027, 128.515, 3 ARD: −58.907, 131.528, 6	2.692 (0.02)	NA
FG papillated texture ♂	3 states	*ARD: −32.851, 79.489, 6 ER: −39.155, 80.387, 1 SYM: −38.373, 83.227, 3	ARD: −33.438, 80.662, 6 ER: −40.505, 83.088, 1 SYM: −39.364, 85.208, 3	10.28 (0.015)	NA
EFA Chelicera cluster membership ♀	3 states	ER: −39.291, 80.691, 1 SYM: −38.245, 83.175, 3 ARD: −37.529, 89.683, 6	*ER: −38.842, 79.793, 1 SYM: −37.948, 82.581, 3 ARD: −37.029, 88.682, 6	2.114 (0.04)	NA
Chelicera measures cluster membership ♀	3 states	*SYM: −24.594, 55.915, 3 ER: −27.145, 56.405, 1 ARD: −22.757, 60.313, 6	SYM: −25.038, 56.803, 3 ER: −27.524, 57.163, 1 ARD: −23.256, 61.311, 6	12.898 (0.02)	NA
Operculum measures cluster membership ♀	3 states	ER: −28.159, 58.446, 1 SYM: −26.05, 58.927, 3 ARD: −22.857, 60.944, 6	*ER: −26.19, 54.51, 1 SYM: −24.3, 55.427, 3 ARD: −23.013, 61.257, 6	13.486 (0.11)	NA
Flagellum complex length (FCL) ♂	4 states	*ER: −46.624, 95.324, 1 SYM: −42.475, 98.7, 6 ARD: −40.416, 112.261, 12	ER: −47.388, 96.852, 1 SYM: −42.754, 99.258, 6 ARD: −40.627, 112.682, 12	15.355 (0)	NA
FG papillated texture position ♂	4 states	*ER: −40.069, 82.215, 1 SYM: −35.827, 85.442, 6 ARD: −34.373, 100.355, 12	ER: −41.552, 85.18, 1 SYM: −37.386, 88.559, 6 ARD: −34.234, 100.078, 12	22.549 (0.07)	NA
EFA Chelicera Cluster Membership ♂	4 states	SYM: −49.824, 113.515, 6 ER: −56.932, 115.945, 1 ARD: −46.078, 124.157, 12	*SYM: −49.823, 113.513, 6 ER: −57.045, 116.171, 1 ARD: −46.037, 124.073, 12	11.761 (0)	NA
Tooth count ♀	6 states	*ER: −57.929, 117.953, 1 SYM: −48.556, 144.256, 15 ARD: −43.653, 290.382, 30	ER: −56.949, 115.992, 1 SYM: −48.758, 144.659, 15 ARD: −43.899, 290.874, 30	11.27 (0.075)	NA
Chelicera measures cluster membership ♂	7 states	ER: −85.879, 173.84, 1 SYM: −73.636, 221.135, 21 ARD: −68.418, 672.336, 42	*ER: −86.34, 174.761, 1 SYM: −75.183, 224.227, 21 ARD: −70.527, 676.554, 42	2.87 (0.015)	NA
Tooth formula ♀	9 states	ARD: −59.594, −99.294, 72 ER: -83.704, 169.503, 1 SYM: −67.957, 588.485, 36	*ARD: −58.588, −101.307, 72 ER: −83.743, 169.582, 1 SYM: −68.017, 588.605, 36	5.03 (0.02)	NA

Due to the intraspecific variation existing among quantitative and qualitative characters considered, we calculated the mean measurement and identified the modal state for each terminal on the pruned phylogeny. For several terminals, however, we lacked information due to some species descriptions being limited to a single sex, thus the corresponding sex may be unknown. Prior to performing character mapping for discrete states, maximum likelihood model evaluations were performed to determine the best-fit model for the characters of interest and associated topology. We fit ASR models using equal rate (ER), symmetric rate (SYM), and all rates are different (ARD) models using the fitDiscrete function found in the ape package. Best-fit models were selected by identifying the model that produced the lowest AICc score.

Maximum likelihood character mapping of discrete states using the best-fit models was executed using the ace function within the phytools package (Revell 2012). For character mapping of continuous states estimated from maximum likelihood, we implemented the anc.ML method within the contMap function also found in phytools as this method can support the usage of missing data.

Phylogenetic Signal

For continuous characters specifically, several measures of phylogenetic signal have been proposed that either employ a spatial autocorrelation approach (Moran 1950, Abouheif 1999) or a Brownian motion model of evolution approach (Pagel 1999b, Blomberg et al. 2003, Fritz and Purvis 2010). Of the available approaches, we estimated Pagel’s λ and Blomberg’s K using the phylosig function in phytools with 999 simulations. For both statistics, values can reflect weak signal (<1), be congruent with Brownian motion (=1), or strong signal (>1). Phylogenetic signal statistics for discrete characters, on the other hand, are limited to 2 options: D-statistic for binary character states (Fritz and Purvis 2010) and the δ statistic which can extend beyond binary states (Borges et al. 2019). Aside from the difference in the number of discrete states necessary for statistical estimation, the latter statistic is estimated by adopting the concept of Shannon entropy (Shannon 1948), which quantifies the amount of variability associated with a random variable. Values of D can indicate an overdispersion of traits (>1), random evolution (=1), Brownian evolution (=0), or highly conserved trait evolution (<0) (Fritz and Purvis 2010). For the δ statistic, on the other hand, the higher value of this statistic represents greater phylogenetic signal (Borges et al. 2019). Thus, for discrete traits, we estimated D for binary traits using the phylo.d function in the caper package (Orme et al. 2013) and for estimation of the δ statistic we used the ape package and the function delta (Borges et al. 2019) for all qualitative traits considered in this study. Each phylogenetic signal statistic is associated with a test that maps the observed phylogenetic signal statistic to a simulated distribution representing no phylogenetic signal (random distribution of traits along a phylogeny) to determine if the observed trait is statistically different from expected under a random model. We calculated this probability by generating a probability distribution of 200 random δ statistics and used the default settings for generating a random D-statistic distribution. All considered statistics for phylogenetic signal and tests were estimated using the best topology (chronogram vs. phylogram) as determined by lowest recovered AICc score (Table 5) for discrete states. For continuous traits, on the other hand, we calculated phylogenetic signal using both a chronogram and phylogram.

Results

Phylogenetic Relationships of Eremobatidae

Original raw sequence reads are submitted to the Sequence Read Archive (SRA: BioProject ID PRJNA982881). The final 75% occupancy alignment matrix consisted of 225 UCE loci with 166 terminal eremobatid taxa and 4 outgroup taxa for a total of 105,239 sites, 45,691 of those sites being parsimony informative. For coalescent-based phylogenetic reconstruction in ASTRAL, we utilized 916 alignments representing the same number of eremobatid individuals and 4 non-eremobatid taxa as the concatenated analysis. Taxon coverage is summarized in Table 1, and tree files were uploaded onto Figshare. Number of UCE loci recovered for each ingroup taxon is included in Supplementary Information 3. Consistent with the findings of the multilocus analysis of Eremobatidae by Cushing et al. (2015) and the phylogenomic analysis of Solifugae (Kulkarni et al. 2023), both concatenated and coalescent-based phylogenetic analyses recovered Eremobatidae with full support (Fig. 4; Supplementary Fig. S5). All speciose genera, specifically Eremobates, Eremochelis, and Hemerotrecha, were recovered as paraphyletic, which is concordant with the phylogenetic hypothesis of Cushing et al. (2015). Horribates and Chanbria, contrary to the other ingroup genera, were rendered as monophyletic. Most of the major clades in both topologies recovered clades with moderate (ML > 70/0.70) to full support (100/1), except for 7 clades. Five of those 7 clades lacked support from the concatenated analysis, however, the same clades recovered strong support (>0.90) for the coalescent-based analysis. Only 2 clades lacked support entirely. The first was the clade possessing a putative new eremobatid genus and Eremochelis andreasana (Muma 1951). In all preliminary and final analyses for this study, this relationship was consistently recovered with low resolution. The second unresolved clade was a recently divergent clade containing Eremochelis morrisi (Muma 1951), Eremochelis striodorsalisMuma 1962, Eremochelis kastoni Rowland 1974, and members of the Hemerotrecha banksi species group.

Fig. 4.

Maximum likelihood topology of Eremobatidae inferred using IQ-TREE 2 and the 75% occupancy concatenated UCE matrix with taxa that rendered >10 UCE loci with outgroups removed. Nodal support for major clades refer to maximum likelihood bootstrap values (ML) on the left and local posterior probabilities estimated based on gene tree quartet frequencies in ASTRAL (MSC) on the right. Tip colors illustrate generic designations as per current eremobatid taxonomy. Male (left column) and female (right column) structures of eremobatid species are shown on the right of the topology. Images are ordered phylogenetically. A) Eremochelis andreasana (♂, DMNS ZA.41919; ♀, DMNS ZA.40822). B) Hemerotrecha marathoni (♂, CNAN loan; ♀, CNAN loan). C) Hemerotrecha xena (♂, DMNS ZA.16469; ♀, DMNS ZA.42101). D) Eremochelis bilobatus (♂, DMNS ZA.17629; ♀, DMNS ZA.17359). E) Horribates spinigerus (♂, DMNS ZA.40086; ♀, AMNH_Holotype). F) Chanbria regalis (♂, DMNS ZA.25443; ♀, DMNS ZA.25442). G) Hemerotrecha serrata (♂, DMNS ZA.16059; ♀, DMNS ZA.22169). H) Eremochelis flexacus (♂, DMNS ZA.16134; ♀, DMNS ZA.23593). I) Hemerotrecha parva (♂, DMNS ZA.42131; ♀, DMNS ZA.42131). J) Eremochelis larreae (♂, CASENT9033550; ♀, DMNS ZA.33721). K) Eremochelis insignatus (♂, DMNS ZA.28260; ♀, DMNS ZA.26390). L) Eremochelis striodorsalis (♂, DMNS ZA.31773; ♀, DMNS ZA.23593). M) Hemerotrecha californica (♂, DMNS ZA.18288; ♀, DMNS ZA.32217). All scale bars represent 1 mm.

Four independent clades comprising multiple Eremochelis taxa were recovered with full support. The first early diverging Eremochelis clade included Eremochelis albaventralis Brookhart & Cushing 2005, Eremochelis cochiseae Muma 1989, and Eremochelis bilobatusMuma 1951. Second, a later diverging and larger clade consisted of Eremochelis kerni Muma 1989, E. giboi Muma 1989, E. flexacusMuma 1963, E. acrilobatusMuma 1962, E. fuscellus, E. nr bechteli Muma 1988, E. branchiMuma 1951, and E. medialisMuma 1951. Next, E. nudusMuma 1963, E. bidepressusMuma 1951, E. oregonensisBrookhart & Cushing, 2002, E. larreaeMuma 1962, and E. undulus Muma 1989 formed a monophyletic grouping. And lastly, the latest diverging and smallest Eremochelis clade included a putative undescribed species, E. sp., and E. insignatusRoewer 1934.

The earliest diverging clade contained species of the Hemerotrecha branchi species group. A more inclusive Hemerotrecha branchi species group clade would have been recovered, yet 1 anomalous taxon, Eremochelis truncus Muma 1987, was nested in between the other Hemerotrecha branchi species group taxa. We hypothesize that the placement of this species could be erroneous due to the lack of UCE data as this specimen only recovered 25 UCE loci.

The other strongly supported clade included members of the Hemerotrecha denticulata species group, however, consistent with Cushing et al. (2015), a male identified in this study as a member of the Hemerotrecha denticulata species group with 506 UCE loci, was nested within Eremobates. Next, Hemerotrecha elpasoensisMuma 1962, Hemerotrecha fruitanaMuma 1951, putative new species Hemerotrecha sp., and other morphologically similar specimens formed a strongly supported clade. Lastly, most of the Hemerotrecha banksi species group representatives formed a monophyletic group, apart from 2 specimens. Hemerotrecha prenticei (DMNS ZA.16782) and Hemerotrecha californica (DMNS ZA.19103) were nested outside of the main Hemerotrecha banksi species group clade. Like E. truncus, we suspect that the placement of H. prenticei is due to missing data as this specimen recovered 54 UCE loci, a value considerably lower than the other Hemerotrecha prenticei specimen, which yielded 252 UCE loci. As for H. californica, we believe this is a misidentification error, as the operculum of H. californica and E. kastoni are similar.

There were several notable species included in this study that was described from single sexes only: specifically, Eremochelis acrilobatus (♀), E. flexacus (♂). Eremochelis acrilobatus and E. flexacus were recovered as a single monophyletic clade with full support for both phylogenetic approaches. Given this result, these 2 species are likely conspecifics; thus, the 2 species were considered as a single species in downstream ASR analyses.

Morphometrics and Morphospaces

For traditional morphometrics, we accumulated nearly 2,700 individual measures for both male and female eremobatid taxa representing the 6 focal genera (listed in Table 1; Supplementary Information 1). A total of 118 females were examined in this study. We measured 7 cheliceral variables and 3 opercular variables resulting in 95 representatives for female chelicerae and 102 representatives for female operculum. Due to the extensive variability present in male cheliceral morphology, we sought to focus our sampling efforts on male representatives. Our comprehensive sampling of male chelicerae totaled 1,734 measurements partitioned among 7 continuous variables for 264 individuals. PC1 explained 95.94%, 97.56%, and 80.80% of the variation in the male cheliceral, female cheliceral, and female opercular measures, respectively. Morphospace plots for all associated measures for each respective sex illustrate widespread overlap between Eremochelis and Hemerotrecha taxa (Supplementary Fig. S6). As for PC2, 2.00% of the variation was explained for male chelicera, 1.52% for female chelicera, and 15.82% for female opercula. Eremorhax, on the other hand, was the only genus that was consistently isolated in morphospace across male cheliceral, female cheliceral, and female opercular measurements. The general pattern regarding the distribution of samples in morphospace generally supported trends in allometric structure. Characteristically large species, like species in the genus Eremorhax, were distributed on the negative end of PC1, whereas comparatively smaller eremobatids, like the putative new genus was distributed to the right of PC1 (Supplementary Fig. S6) for male chelicerae, female chelicerae, and operculum measures.

We generated 237 male chelicera outlines, 86 outlines of female chelicera, and 99 outlines of a single operculum plate for a total of 436 outlines for EFA (Supplementary Information 2). For the male data set, 16 harmonics explained 99.9% of the variation, resulting in 64 PC axes. The female chelicera data set recovered 19 harmonics, and 13 harmonics for the opercula data set summarizing 99.9% of the variation, resulting in 76 and 52 PC axes, respectively. PC1 explained 78.83% of the variation in male chelicera, 89.13% for female chelicera, and 57.74% for female opercula (Fig. 5). PC2, on the other hand, explained 11.06% of the variation in male chelicera, 7.29% for female chelicera, and 29.31% in female opercula (Fig. 5). The examination of shape for the focal structures corroborated the overall trend observed from the traditional morphometrics data such that large species, such as those in the genus Eremorhax, were oriented on the opposite end of smaller representatives. Moreover, we observed a broad overlap for all considered shapes between Eremochelis and Hemerotrecha.

Fig. 5.

Morphospace plots based on PCA analysis of EF coefficients and GMM cluster assignments recovered from the first 3 PCs as input. Individuals with uncertain cluster membership are displayed in gray. Top row pertains to male chelicera, central row depicts female chelicera, and bottom row refers to female opercula. All PCA plots summarize PC1 and PC2 and illustrate: A) male cheliceral shape by genus, B) male cheliceral shape by cluster membership with mean shapes of each cluster to the right (cluster 1 was not supported after implementing the 0.95 uncertainty threshold), C) female cheliceral shape by genus, D) female cheliceral shape by cluster membership with mean shapes of each cluster to the right, E) female opercula shape by genus, F) female opercula shape by cluster membership with reconstructed operculum from the single opercula median shape. All shapes are not to scale.

The linear model test for allometry regarding the relationship between size and shape resulted in an R² value of 0.111, a Z value (effect size) of 5.395 and a P-value of 0.0009. This result suggests that there is a moderately positive relationship between size and shape.

Model-Based Cluster Designations

Results of the GMM on male cheliceral measurements supported 7 independent clusters in morphospace, which was the largest number of clusters among all the other data sets analyzed (Supplementary Fig. S6; Table 3). Under the 0.25 uncertainty cut-off, otherwise uncertainty values greater than 0.75, a total of 39 taxa could not be confidently assigned to a cluster in the male data set. The stricter threshold of 0.05 (uncertainty > 0.95) resulted in the uncertain placement of 128 male individuals. Cluster 1 and Cluster 2 yielded the most stable clusters as these 2 clusters did not exhibit a dramatic loss of members across the 2 applied thresholds (Table 5). Both GMM analyses of female cheliceral and opercular data sets suggested 3 inherent clusters in morphospace (Supplementary Fig. S6C and F) and exhibited less of an uncertainty loss compared to the male chelicera measurements data set. Under the 0.25 cut-off for female chelicerae, there were 6 individual representatives that could not confidently be placed into a respective cluster, and 4 for the female opercula. With the stricter uncertainty threshold, 24 individuals for female chelicerae and 18 female opercula, could not be confidently assigned a cluster designation. For the female chelicerae, Cluster 2 was the only stable cluster across the enforced thresholds. Clusters 1 and 2 were the most consistent clusters for female opercula with a loss of less than 10 individuals.

The results of the GMM clustering analysis of shapes using the first 3 PCs ranged from 2–4 clusters (Fig. 5; Table 3). First, like the results for traditional measures, male cheliceral shape data recovered the most clusters with support for 4 clusters as optimal. For male cheliceral shape, Cluster 2 was the most stable across the applied thresholds (Table 3). Of the 237 male chelicerae outlines generated, 156 individuals lacked support for their respective cluster, which was 66% of the input data. Moreover, no individuals were supported for Cluster 1 membership after we applied the stringent threshold (Table 3). The GMM algorithm for the female cheliceral shape data set supported 3 clusters as the best-fit model (Fig. 5D), with 7 individuals unable to be confidently placed within a cluster under the 0.25 uncertainty cutoff and 33 for the 0.05 threshold. Only a single cluster, cluster 2, was virtually unchanged across thresholds. And lastly, the female opercula shape data supported 2 clusters given the best-fit model (Fig. 5F), with both clusters remaining relatively consistent across the different thresholds (Table 3).

Phylogenetic Signal and Character Evolution of Quantitative Variables

For males, we quantified 3 chelicera single measures, 2 chelicera ratios, and the PC1 loadings resulting from our EFA (Table 4). For females, on the other hand, we used 6 single measurements, 3 ratios, and 2 PC1 loadings (Table 4). The spectrum of variables for both males and females ranged from displaying strong signal (1.791; PC1 male chelicera) to weak signal (0.849; L/FFH female), with few variables producing indices that suggest Brownian motion evolution under Blomberg’s K. Statistics for Pagel’s λ ranged from operculum width indicating the strongest signal of the traits examined (1.286) to PC1 loadings for operculum (0.881) with the weakest signal (Table 4). Our associated tests for the 2 quantified phylogenetic signal statistics suggest that all continuous variables, in the context of our chronogram, did not agree with a random model of trait evolution (Table 4).

Phylogenetic Signal and Character Evolution of Qualitative Variables

We quantified 7 discrete character states for male chelicerae and 2 qualitative variables for females. To explore the utility of cluster assignments in a phylogenetic context, we treated each cluster assignment resulting from the analysis of the EFA and traditional measures data sets for male chelicera, female chelicera, and opercula as discrete states. For cluster membership identity, when variable, we maintained the modal cluster identity for each terminal. Therefore, in total, we attained 9 discrete variables for males and 6 for females. Most model comparisons supported the equal rates (ER) model as the best-fit model (Table 5). The FN and the central carina variables were equally supported by the ER and symmetrical rate model (SYM); however, we implemented the ER for these variables for consistency.

For males, the δ-statistic ranged from 2.216 for FN to 66.928 for the central carina male characteristic (Table 5). However, despite the large phylogenetic signal for the central carina character, this value was not different than a random distribution of this statistic (P = 0.905), therefore, this character exhibits no phylogenetic signal under this model. The other 8 male discrete characters quantified in this study are suggested to have some degree of phylogenetic structure given that the probability of observing the estimated phylogenetic statistics in a random distribution of statistics was quite unlikely (Table 5). The FG papillation position, or the position of the rough texture found in the FG, recovered the highest δ-statistic, suggesting that this was the discrete character with the highest phylogenetic signal.

The phylogenetic signal statistics for females ranged from 2.114 for EFA female chelicera cluster membership to 22.086 for female cheliceral measures cluster membership (Table 5). However, 2 out of the 6 characters quantified resulted in phylogenetic statistics with high probability of occurrence under the random (no phylogenetic signal) model. The 4 characters that resulted in low probability of being in accordance with the random model were (i) female EFA chelicera cluster membership, (ii) female chelicera measures cluster membership, (iii) female tooth count, and (iv) female tooth formula. As a result, female cheliceral measures cluster membership was the variable that reflected the highest phylogenetic signal.

Among the discrete characters, we coded 4 binary character variables and subsequently estimated 4 D-statistics. Of the 4 binary character states, all states demonstrated phylogenetically conserved evolution, meaning that traits we examined were more similar among closely related taxa, as determined by the recovered negative D values (Table 5; Fritz and Purvis 2010). Contrary to the δ-statistic, the D-statistic and associated test suggested that opercula cluster membership is highly conserved and is parallel with Brownian phylogenetic structure (Table 5).

Ancestral State Reconstructions of Quantitative Characteristics

Ancestral state reconstructions (ASR) of mean PC1 loadings for male chelicerae illustrate 4 independent evolutions for negative PC1 values (Fig. 6). The specific clades that are supported by a negative PC1 ancestral state were the putative new genus, representatives of the Hemerotrecha branchi group, representatives of the Hemerorecha simplex group, and a most recently derived clade that includes Eremochelis nudus to Hemereotrecha kaboomiBrookhart & Cushing 2008. The clade that includes sister taxa, Eremochelis cochiseae and Eremochelis albventralis were taxa that had slightly larger character states than the observed negative PC1 value clades, with sister taxon, Eremochelis bilobatus showing a character state toward the median of the PC1 range (Fig. 6). The clade that included Eremochelis imperialis to Hemereotrecha nr neotena all reflected moderate to large PC1 values. Additionally, the clade that included members of the Hemerotrecha denticulata group also resulted in a similar pattern. Mapping of PC1 loadings of male chelicera and female opercula illustrated similar evolutionary trends, however, many of the terminals in our pruned phylogeny were missing data for female opercula. For example, species like Eremochelis albaventralis, Hemerotrecha texanaMuma 1951, Eremochelis fuscellus, Eremochelis nudus, among others are species that were described from males only. Therefore, we suspect that the results of this ASR were affected by the amount of missing data.

Fig. 6.

Ancestral state reconstruction of EFA coefficients. A) PC1 from the PCA of male chelicera EF coefficients and B) PC1 from the PCA of female opercula EF coefficients.

Mapping of male and female cheliceral length suggests that the eremobatid ancestor chelicera was of a moderately smaller length (Supplementary Fig. S7 and S8). Shorter chelicera length evolved 4 independent times among the putative new genus, the Hemerotrecha branchi species group, the Hemerotrecha simplex species group, and the earliest diverging clade in our phylogeny (Supplementary Fig. S7 and S8). From the small cheliceral length ancestral condition, our results modeled an increasing cheliceral length for both males and females in the clade that included Eremorhax joshui, of which was the taxon with the greatest cheliceral length.

ASR of Male Qualitative Characteristics

The presence of the FG, the groove located on the prolateral side of the FF in male chelicera, was the suggested ancestral state condition for all ingroup eremobatids with full probability (Fig. 7A). Through time, the probability of FG absence increased for many of our observed major clades. The absent state was coded for 5 clades comprising more than 2 species, 1 clade with 2 independent species, and 2 singleton lineages; thus, the results suggest that the FG was lost 8 independent times. For the present state, on the other hand, 5 multispecies clades and 4 single species lineages possess the present state (Fig. 7A). Of the 56 terminals in our working phylogeny, 35 individuals were observed for an absent state, therefore making it the most prevalent state among the taxa considered in our study (62.5%).

Fig. 7.

Ancestral state estimation and mapping of binary states of male A) male flagellar groove and B) fondal notch. Absent = no flagellar groove/no fondal notch, present = flagellar groove/fondal notch.

The male FN, found only in male chelicerae and defined as the indentation extending posteriorly past the fondal axis (Supplementary Fig. S3), independently evolved 9 times among the ingroup taxa (Fig. 7B). More than half of the terminal taxa (35 of 56; 62.5%) were scored as possessing an absent FN character state. Within our phylogeny, 7 clades with more than 2 species and 2 lineages representing a single species were coded for an absent FN character state. The probability of an absent FN character state was the highest for nearly all nodes, except for 4 clades and 3 singleton lineages. Regarding clades, first, a clade containing species Eremochelis kerni, Eremochelis giboi, Eremochelis flexacus, Eremochelis branchi, etc. was a clade that supported an ancestral node with strong probability for a present FN character state (Fig. 7B). Second, the clade comprising species belonging to Eremobates and the single considered species for Eremorhax, was inferred to having the FN character state as the most probable ancestral state. And lastly, the smaller clades consisting of Eremochelis insignatus and Eremochelis sp. and the other clade represented by Eremochelis striodorsalis and Eremochelis kastoni were supported by a FN character state. Overall, from the current character states for the terminal taxa in this study, the FN was gained 9 times.

The medial tooth, moveable finger (MM) was coded for 3 states: present (P), reduced (R), and absent (A; Supplementary Fig. S2). According to the δ-statistic and associated test, this character resulted in weak phylogenetic signal with a value of 2.692 (Table 5). The most probable ancestral state for this character was the P state, with the A state being the second most likely, and the R state as the most improbable state. Of the 56 taxa considered for our ASR analysis of the MM, 26 taxa (46%) were coded as P, 16 were coded for A (29%), and 14 were recorded as R (25%; Supplementary Fig. S9). The most inclusive scored character was the P character state due to this state being a consistent shared, derived condition among 4 clades throughout the phylogeny. Specifically, members of the Hemerotrecha branchi species group, sister species Eremochelis albaventralis and E. cochiseae, E. branchi and E. medialis, E. insignatus and E. nr insignatus, and taxa belonging to the Hemerotrecha banksi species group (e.g., Hemerotrecha banksi, H. prenticei, H. californica, H. vetteri, etc.). Both the R and A conditions only represented 2 clades respectively. The A-coded state was a shared, derived state for sister species Eremochelis kerni and E. giboi, and for E. nudus, E. bidepressus, E. larreae, and E. undulus. Lastly, the R character state was a synapomorphy for Chanbria and associates of the Eremobates palpisetulosus species group (Eremobates leechi Muma & Brookhart 1988 and E. gracilidensMuma 1951) and Eremorhax joshui Brookhart & Muma 1987.

Within Eremobatidae, the flagellar groove (FG), specifically the presence/absence of the flagellar groove, shape, and position, was a foundational character set used for establishing and delimiting generic boundaries. However, after our initial survey of the FG under SEM we observed texturing inside the FG (Fig. 3B). This texturing, we describe here as papillated, had not yet been previously reported for Eremobatidae. With one exception, the papillated texturing was observed to be present primarily inside the FG and varied in position within the FG. Therefore, we explored the possibility that this texturing may be useful for taxonomic revision and tested the phylogenetic utility of the papillated texturing and position. FG papillated texture position was coded as 4 states: 1 = smooth FG/texturing absent, and 2 = papillated texturing, and 3 = no FG/smooth FF. Meanwhile, the FG papillated texture discrete variable was coded for 3 states: 1 = proximal texturing only, 2 = distal texturing only, 3 = continuous texturing, both proximal and distal, and 4 = no texture. The ancestral trait for male eremobatids in this study was the smooth texture condition with FG absent (Supplementary Fig. S10). The FG papillated texturing, on the other hand, was gained 9 independent times and many of those same gains were scored as being located on the distal end, toward the apex of the FF (Supplementary Fig. S10).

The flagellum complex (FC) is composed of morphologically distinct bristles that extend from the prolateral side of the FF. To explore the phylogenetic utility of the flagellum complex, we discretized flagellum complex length (FCL) and quantified this relative length with respect to FF. We coded 4 states for FCL: 1 = FC is ¼ the length of the FF, 2 = FC is ½ the length of FF, 3 = FC is ¾ the length of the FF, and 4 = FC starts from the manus (Supplementary Fig. S11). None of the ingroup taxa were observed to have an FC that extended the full length of the FF. Our ASR analysis of this discrete character suggests that the ½ length state was the ancestral state for all major clades with high probability (Supplementary Fig. S11). Taxa that coded for the FC starting at the manus only represented individual lineages, whereas the other character states were more inclusive, delineating multiple clades and species lineages. The quarter length was suggested to be the most prolific character state for extant taxa, whereas half-length was the second most widespread state.

We conducted 2 respective GMM analyses on our measurements and EF coefficients data set to determine if similar shapes and measures, determined by cluster membership, are phylogenetically informative. We treated each corresponding cluster information as discrete variables. The cluster designations using the EF coefficients for male chelicerae, on the other hand, were recovered as a phylogenetically informative trait that did not agree with a random model of evolution (Table 5; Fig. 8). The cluster memberships for male cheliceral measurements, on the other hand, was suggested to have weak phylogenetic signal compared to EFA cluster membership (2.869 vs. 11.761). The result of the ASR analysis for male measurements only supported a few extant clades with shared, derived cluster membership (Supplementary Fig. S12). Specifically, a clade comprising Hemerotrecha branchi species group members: Hemerotrecha macraMuma 1951,, H. marathoniMuma 1962,, H. sevilletaBrookhart & Cushing 2002,, H. cornutaBrookhart & Cushing 2002, a clade with Hemerotrecha denticulata species group members (Hemerotrecha denticulata, Hemerotrecha delicatula Muma 1989, Hemerotrecha parva Muma 1989, Hemerotrecha carsonana Muma 1989), and another clade with Hemerotrecha banksi species group members: H. prenticeiBrookhart & Cushing 2008, H. californica Banks 1899, H. pseudotruncata, and H. kaboomi.

Fig. 8.

Ancestral state estimation and mapping of Gaussian mixture modeling (GMM) cluster membership onto our phylogram for ingroup male eremobatid species inferred from Elliptical Fourier coefficients for chelicerae. Character states 1−4 refer to the 4 distinct clusters in morphospace (Fig. 5B).

ASR of Female Qualitative Characteristics

Four out of the 6 quantified discrete female variables demonstrated some degree of phylogenetic signal, all of which pertained to variables of female chelicerae. The ASR analysis of the EFA cluster association supported multiple independent lineages of each cluster membership, however, this result lacked a broad inclusion of multiple species (Supplementary Fig. S13), thus was weakly informative (Table 5). Of the 39 taxa considered, 11 taxa (28%) were coded as Cluster 1, 10 were coded for Cluster 2 (26%), and 18 were recorded as Cluster 3 (46%). For the cluster data set resulting from female measurements, the ancestral state was estimated to be Cluster 1 with the majority of the ingroup taxa pertained to this cluster designation (Supplementary Fig. S14). Only species, Eremorhax joshui, belonged to Cluster 3.

Coding for tooth count and tooth formula both resulted in 9 unique states (Supplementary Fig. S17). The most coded state for tooth count was 7 (52%) and most occurring state for tooth formula was FD-2FSD-FM-2FSM-FP (30%). Although these 2 characters were moderately phylogenetically informative and not in alignment with a random model of evolution (Table 5), there was 1 inclusive clade that supported many species with a single shared character state (Supplementary Fig. S17). This clade included Eremochelis kerni, E. giboi, E. flexacus, and E. branchi, both of which were unique by the most common character states describing female dentition mentioned above.

We objectively characterized opercular shape and opercular measures in morphospace using GMM. Opercular shape was classified into 2 shape clusters (Fig. 5F; Supplementary Fig. S15) and 3 were supported for traditional measures (Supplementary Figs. S6 and S16). However, although the ASR results for traditional opercular measures demonstrated several instances of shared cluster ancestry, both these characters were in accordance with a random model of evolution.

Discussion

Eremobatidae Phylogeny

Here we present a phylogenomic approach and a new phylogenetic hypothesis for the North American camel spider family Eremobatidae. Our results, we believe, give implications for new generic designations such that the genera are represented by clades of shared ancestry. Despite our limited taxon sampling for other Eremobates species groups, and lack of representation for Eremothera and Eremocosta, we expect that the well-supported clades recovered in our phylogenomic hypotheses will hold consistent as some of these clades are supported by morphological traits. The clade consisting of Horribates, Eremobates, Eremorhax, Eremochelis imperialis, Hemerotrecha serrata, and Eremochelis kerni to Eremochelis branchi all comprise taxa that are considerably large eremobatid species compared to the other taxa included in this study. Members of this clade tend to be 20 mm or more in body length, with cheliceral length being above 4.0 mm, on average. Although we did not sample some of the characteristically large species of Eremocosta in our analysis, we suspect that Eremocosta will fall within this large body size eremobatid clade, as well as Eremothera.

Several clades were generally supported by the morphological characters we examined here, and we believe that this is preliminary evidence that may support new generic-level placements within Eremobatidae. At this time, given our data, taxon sampling and results, we suggest that the Hemerotrecha branchi species and the Hemerotrecha denticulata species group merits elevation to a generic level rank in future taxonomic works. There are several implications that suggest that the Hemerotrecha banksi species group could be elevated to a generic rank, however, further evidence is required to determine whether Eremochelis striodorsalis and Eremochelis kastoni could be included in the same generic level rank.

We suggest the following synonymies for future taxonomic consideration: Hemerotrecha branchi and Hemerotrecha xena, distinguished by the presence or absence of papillae (Muma 1951), occur in various combinations within the same clade in our analysis, therefore are likely the same species with variable papillar morphologies. Second, Chanbria regalisMuma 1951 and Chanbria rectusMuma 1962, distinguished by the slight dorsal curvature of the FF in C. regalisMuma 1951, are also probable conspecifics (Fig. 4). Additionally, E. flexacus and E. acrilobatus are complementary species of the male and female, respectively. And finally, E. oregonensis a morphologically similar species to E. bidepressus and was found nested within multiple E. bidepressus representatives, thus should be synonymized.

Evolution of Historically Delimiting Characters

The absence of the FG has been historically used to place eremobatid species into the family Hemerotrecha (Muma 1951). The results of our ASR analysis of this trait suggest that the loss of the FG occurred multiple times throughout the evolutionary history of Eremobatidae. However, the FN and FG appear to be correlated traits. The FN, when present, is generally associated with the presence of an FG and vice versa (Fig. 7). Moreover, our results indicate that the FG is the most probably ancestral state. We speculate that the FG is a derived condition within Eremobatidae due to all males of the other extant solifuge families, including sister taxon Gylippidae (Kulkarni et al. 2023), are lacking an FG and instead possess a flagellum. The absence of an FN, on the other hand, resulted in being the most probable ancestral state, and thus is likely the solifuge plesiomorphic condition since all other non-eremobatid solifuge families lack an FN. The presence of the FN may have evolved to facilitate deep insertion of the FF into the female reproductive tract (Bird et al. 2015), during contact phase of eremobatid mating (Muma 1966a, Rowsell and Cushing 2020). However, suggested by the apparent homoplasy of this trait within Eremobatidae, and absence of the trait in the rest of Solifugae, the FN and the FG may not be primary structures for promoting the reproductive success within Eremobatidae. If the FG was analogous to the flagellum in solifuges, this observation would corroborate the assertion of several researchers who have acknowledged that the flagellum was likely not related to reproductive success (Heymons 1902,Cloudsley-Thompson 1961, Junqua 1962, Lawrence 1963, 1965). However, the presence of a flagellum is conserved throughout Solifugae, except for Eremobatidae, which implies a probable necessary function of the flagellum. Moreover, the evolutionary conservation of the FG and FN in Eremobatidae many clades would highlight some degree of importance of these traits for reproduction.

As an objective approach to establishing qualitative gross operculum morphologies, we assigned each sample operculum to a morphological cluster determined by the distribution of morphologies in morphospace using EFA cluster membership. However, this approach resulted in ambiguous phylogenetic patterns as multiple independent lineages recovered the same cluster designation (Supplementary Fig. S16). Therefore, we argue our attempt to discretize gross operculum morphologies through EFA cluster membership for operculum form is not an ideal variable to be used for taxonomic revision within Eremobatidae. This result may demonstrate that female operculum shape is evolving neutrally and likely has not played a strong role in the diversification of eremobatid taxa. Alternatively, this result may suggest that the female internal reproductive anatomy, rather than the external female operculum shape, could be a major factor in driving male chelicera variation. The muscle mechanism as to which sperm displacement happens has been suggested to depend on the structure of the female genital system by previous researchers, thus researchers have speculated that the female genital system morphology is important in cryptic female choice (Hellriegel and Ward 1998). Furthermore, the anatomical reconstruction of the female genital tract of 3 species of Galeodes suggests that the female genital tract morphology is highly variable (see Figs. 2−5, Peretti et al. 2021). Therefore, due to the known behavior of FF insertion into the reproductive tract, and the lack of phylogenetic informativeness for objective operculum shape, we hypothesize that the female internal anatomy, rather than the external genital plate morphology is strongly influencing male chelicerae morphology at the species level. In eremobatids, the transfer of sperm involves gonopore-to-gonopore contact (Muma 1966a), unlike other solifuge families in which transfer involves cheliceral sperm capture from a deposited sperm packet (Wharton 1986, Hrušková‐Martišová et al. 2010). Bird et al. (2015) suggested that, due to this lack of sperm capture in Eremobatidae, that perhaps the chelicerae evolved for other mating functions, such as for interacting with the female oviduct. We encourage future researchers to further investigate female genital tract morphology as this may uncover important traits that could be driving male chelicerae morphologies within Eremobatidae.

Considerations of New Qualitative Characters

In total, we identified 4 qualitative male characters used to describe cheliceral characteristics and 2 female characters associated with cheliceral tooth morphology. Of all the characters we considered and mapped onto our phylogeny, none of the exploratory characters were parsimonious and demonstrated a complex evolution of traits. Many of the characters examined were only consistently conserved for some clades with multiple independent gains and losses. Therefore, on a macroscale, none of these characters, nor the traditionally used characters should be treated as stand-alone for delimiting generic boundaries. Instead, we suggest that future taxonomic revisions should consider a consensus among a combination of morphological characteristics for well-supported clades.

The 4 newly proposed male qualitative characters that were investigated in a phylogenetic context were relative flagellar complex length (FCL); presence or absence of the medial tooth, moveable finger (MM); FG papillated texture; and flagellar groove papillated texture position. All investigated characters quantify male chelicerae variation, and our results suggest that all have some degree of phylogenetic signal. Of those characters, we noticed several consistent patterns for the ingroup taxa in our pruned phylogeny. These morphological patterns and combinations could be considered for future taxonomic revisions of this group. First, taxa that have an absent FG and FN character state, or simply no FG, generally lack the papillated texturing on the FF. The same pattern was observed for the opposite combination of character states in such a manner that, when the FG and FN are both present, there is a high probability that the species has a papillated texturing within the flagellar groove and texturing located distally. Regarding FCL, taxa with a present or reduced MM, an absent FG and FN, tend to have the flagellar complex emerging from half or more the length of the FF. Clades that represent this pattern are the Hemerotrecha branchi and banksi species groups and Chanbria. Notably, independent of MM, species with an FG and FN character state are disposed to having a flagellar complex starting from at least a quarter length of the FF. Although we observed these general patterns among our focal taxa, we acknowledge that there are exceptions to these patterns and perhaps quantifying other cheliceral characteristics, such as flagellar complex setae type, will help to elucidate the encompassing picture regarding male cheliceral trait evolution.

Previous studies on solifuge mating behavior give biological implications for the patterns mentioned above. It is known that males use their chelicerae to clasp, knead, and insert the FF into the female genital tract, pre- and post-sperm transfer (Muma 1966b, Rowsell and Cushing 2020). Despite this role in copulation, however, the functional significance of the morphological traits associated with the male chelicerae is poorly understood. The papillated texture within the flagellar groove, described for the first time here for eremobatids, may be associated with sperm packet adhesion, as a textured surface may assist in adhesion to the FG upon insertion of the FF into the female reproductive tract. In non-eremobatid taxa, it has been suggested that the male flagellum may have mechanoreceptive properties or may be affiliated with chemoreception (Lamoral 1975). Although eremobatids lack a flagellum, the presence of the papillated texture in some eremobatid taxa, and a similar texture observed on the flagellum of Branchia potensMuma 1951 and noted in Peretti et al. (2021) for Otlacola chacoensis, may be evidence for a sensory function beyond simply adhesion. However, to test this functional hypothesis, future studies are required to examine the ultrastructure of these papillae to verify if neurons are present to verify a sensory function.

For those species that lack an FG and papillated structuring, this could suggest that males lacking these traits are relying on other means for ensuring successful mating. For example, the Hemerotrecha branchi species group can be identified by dentition on the FF, lack of an FG, and enlarged, flattened apical flagellar complex setae (Muma 1951). These highly specialized setae are a synapomorphy for the group and exhibit a scale-like texture located dorsally. If this texture is analogous to the papillated texturing found in the FG in some species, then perhaps these setae are responsible for assuring sperm packet adhesion or are sensory modifications. Moreover, if the specialized flagellar complex setae play such function, then it further affirms their importance for copulation, especially for those species that lack an FG and texture.

The length of the flagellar complex (FC), a relative measure from which the setae begin to emerge from the FF, may also be a testament that emphasizes the importance of the specialized setae. Consistently, representatives of the Hemerotrecha branchi species group and the H. banksi species group had a flagellar complex length longer than a quarter length of the FF and, notably, have specialized setae that can be considered synapomorphies for the clade. Species that were placed in the H. banksi species group were described to have flagellar complex apical and subapical bristles that are flattened and broad (Muma 1951). Inspection of these bristles under SEM, however, shows that these prolateral bristles are tubular proximally, then distally taper to a semi-flattened shape. These large setae are nonetheless smooth; thus, it may de-emphasize the role of sperm packet adhesion to the FF in this group. In both these species groups, and other clades like Chanbria with a flagellar complex covering most of the FF, the limited FF area for insertion of the FF into the female could suggest that these highly modified setae are used to further perforate the genital opening to ensure a deep introduction of the sperm packet. Currently, however, it is largely unknown whether or how these specialized setae interact with the female, if they are also inserted into the female reproductive tract, or if they are solely used to interact with the sperm packet for insemination.

The lack of phylogenetic informativeness for tooth count and tooth formula of female chelicera further indicates many aspects of female chelicera are not useful in delimiting taxa at deeper phylogenetic levels within Eremobatidae. Although this has been suggested previously by other researchers (Bird et al. 2015), we decided to analytically test if the tooth variation exhibited among eremobatid female chelicera was phylogenetically informative. Our results suggest that female cheliceral tooth morphology is not under selection and female chelicerae tooth variation is phylogenetically unreliable. Given this lack of support for most of the discrete variables for females quantified in this study, and the obscure phylogenetic patterns resulting from EFA chelicera cluster membership, this suggests that attempts to discretize female variables can be subjective, thus ultimately misleading in a phylogenetic context. Female continuous variables on the other hand, we argue, can be informative at deeper levels (e.g., Supplementary Fig. S8). Our results suggest that continuous variables, variables that are specifically related to size, have the most weight in supporting major clades.

Evolutionary Implications of Measurements and Size

Character conceptualization by establishing discrete states is difficult in the context of solifuges as this approach omits intermediate states. Such states can have important implications for establishing a complete evolutionary picture and for communicating the extent of intraspecific variation. Our PCAs of the anatomical structures and the examination of phylogenetic statistics of the same continuous characters suggest that body size, as determined by the anatomical structures we examined here as a proxy, may have played a crucial role in the evolution of Eremobatidae. Additionally, we found an allometric signal between size and shape within our male cheliceral dataset.

From our morphospace plots, we recognized a general pattern regarding allometry. First, for the male chelicera measures, the representatives with the larger measures were oriented toward the negative end of PC1 (Supplementary Fig. S6A). Moreover, the cluster designations inferred by GMM for Cluster 1 and Cluster 2, principally comprised members with the biggest chelicera. Many of these individuals encompass representatives of Chanbria, Eremorhax, Eremochelis branchi, Hemerotrecha serrata, and other larger species of Eremochelis, all of which are larger than 20 mm in body size. Although comparatively smaller in sample size, we observed a similar pattern for female chelicera and operculum measures (Supplementary Fig. S6E and F), where allometrically larger structures were oriented toward the negative end of PC1 and such structures pertained to larger eremobatid species. Due to this consistent pattern, we suspect that allometry may play an important role within Eremobatidae at deeper taxonomic levels.

We quantitatively summarized eremobatid structures using an EFA and scaled each shape by actual size. Morphospace plots resulting from this approach were consistent with the trends observed with the traditional morphometrics data. There was considerable morphological overlap for male cheliceral shape and size for Eremochelis and Hemerotrecha (Fig. 6), which likely corroborates the polyphyletic result in our phylogenomic analysis. The overlap dispersed with an increasing PC1 with larger chelicera, positioned toward the positive end of PC1. We recognized a similar pattern for both female chelicera and opercula; however, the position of larger structures was opposite such that larger structures were oriented toward the negative end of PC1. Furthermore, there was also extensive morphological overlap for both female structures considered here.

The morphological overlap could suggest convergence of morphological shape and size among the ingroup taxa investigated here. This speculation is corroborated by the ASR analysis of PC1 (Fig. 6, Supplementary Fig. S6) and due to the consideration of size in cheliceral shapes, a similar pattern is reflected in EFA cluster membership for males (Fig. 8). Apart from Eremochelis andreasana, an early diverging species, closely related species among the clades of the Hemerotrecha branchi species group, H. banksi species group, H. simplex species group, and the clade that consisted of Eremochelis nudus and relatives, are more like each other in small body size than to other species with drastically contrasting allometry. The clade consisting of Eremochelis cochiseae, E. albaventralis, and E. bilobatus, the clade comprising Eremochelis imperialis to Hemerotrecha nr neotena, and the clade representing the Hemerotrecha denticulata species group (e.g., Hemerotrecha denticulata, H. parva, H. delicatula) all showed a probable ancestral state of having members with a moderate cheliceral length (Fig. 6A). The ASR analysis of PC1 for female chelicera and opercula reflected a similar pattern for the clades mentioned above (Fig. 6B, Supplementary Fig. S8); however, due to the unavailability of female conspecifics for some of the ingroup taxa in our study, therefore enhancing the effect of missing data, opercula PC1 values showed weak phylogenetic signal.

Several arachnid studies have hypothesized that large-bodied arachnids tend to be associated with arid environments. In a review of desert adaptation in spiders, Cloudsley-Thompson (1983) stated that desert-adapted spiders are typically large and size is hypothesized to be an evolutionary adaptation to prevent evaporative loss due to their small surface area to volume ratio (Cloudsley-Thompson 1983). Likewise, in a study of trapdoor spiders, Coyle and Icenogle (1994) observed that smaller-bodied Aliatypus specialize in mesic habitats and suggested that smaller individuals are prone to desiccation in dry environments due to their large surface area to volume ratio (Coyle and Icenogle 1994). Conversely, save 1 species examined in the empirical study, larger-bodied Aliatypus were affiliated with dry environments; therefore, Coyle and Icenogle (1994) postulated that the association of body size with environment type was an evolutionary adaptive shift. Main (1978) similarly mentioned that the large body size in Australian trapdoor spider, Anidiops villosus (Rainbow), was crucial for the success of this species in semiarid environments. Aside from these arachnid observations, the relationship of body size has been documented in other xeric adapted, arthropod taxa. In trogid beetles, arid-adapted Omorgus species are generally larger than mesic trogid beetles from the genus Trox (Le Lagadec et al. 1998). Given our results, we encourage future studies to investigate the mechanism behind the ecological influences of body size and the possible evolutionary relationship of size and habitat association not only within Eremobatidae but also in terrestrial arthropod systems.

Conclusion

In summary, our study tested the phylogenetic importance of traditionally utilized characters in eremobatid taxonomy and explored potentially new, phylogenetically informative characters to give insight into anticipated taxonomic revisions of this group of understudied arachnids. In this study, we found that many of the traditionally used traits initially used to establish generic boundaries recover most genera, especially the speciose eremobatid genera, as paraphyletic on the molecular level. Of all the traits examined here, none of the traits were conserved toward the root of the phylogeny and are consistently convergent, revealing a complex evolution of traits. We also found that cheliceral size is a good predictor of shape and attempts to objectively classify discrete shapes may not be useful for taxonomic delimitation in this group. However, we believe that our exploratory study informs new descriptions of generic boundaries beyond the traditional generic designations and facilitates the generation of new biological hypotheses regarding trait evolution within Eremobatidae. Moreover, in conjunction with the morphological patterns observed here and the phylogenetic relationships inferred from our analyses, we anticipate that this contribution will aid future generic-level revisions.

Acknowledgments

The authors would primarily like to thank John O. Brookhart for his mentorship, advice, comments, and for sharing his expertise of the group. Thanks to Dr. Carlos Santibañez-Lopez and Dr. Robert Kallal for helpful conversations and guidance on conducting the Elliptical Fourier analysis. We are grateful to Dr. Shahan Derkarabetian for his helpful advice and library prep training. Thanks to Dr. Lauren Esposito at CAS, former curator Dr. Oscar Francke and Dr. Edmundo Gonzalez-Santillán of the Colección National de Arácnidos at UNAM, Peter Oboyski at the Essig Museum, William Clark at the Orma J Smith Museum of Natural History, Dr. Maria Luisa Jiménez at CIBNOR, Michael Wall at the San Diego Natural History Museum, and Felipe Soto-Adames at the Florida State Collection of Arthropods for providing loans for solifuge specimens. Thanks to current members of the Cushing lab, Richard Ryan Jones and Goran Shikak, for their provided edits, suggestions, and advice, of which were invaluable to improving this manuscript. Thanks to Jaír Rojas Castillo, Diana Laura Batista Perales, Oscar Bejarano Mendoza, Edmundo Gonzalez-Santillán, Wendell Icenogle, Zachary Valois, and Marshal Hedin for assistance in specimen acquisition. Thanks to the BROADS, Chayotes, and the affiliated soccer community for their continuous support in the first author’s life. And finally, thanks to the three anonymous reviewers for providing valuable suggestions and edits to help strengthen this manuscript. This project was part of the first authors’ dissertation work and was supported by NSF Grant DEB-1754587 awarded to Dr. Paula Cushing and Dr. Matthew Graham.

Specimen Collection Statement

The authors attest that all legal and regulatory requirements, including export and import collection permits, have been followed for the collection of specimens from source populations at any international, national, regional, or other geographic level for all relevant field specimens collected as part of this study.

Author Contributions

Erika Garcia (Conceptualization [Lead], Data curation [Lead], Formal analysis [Lead], Investigation [Lead], Methodology [Lead], Project administration [Lead], Visualization [Lead], Writing—original draft [Lead]), Quincy Hansen (Data curation [Supporting], Investigation [Supporting], Writing—review & editing [Supporting]), and Paula Cushing (Data curation [Supporting], Funding acquisition [Lead], Investigation [Supporting], Project administration [Lead], Resources [Lead], Supervision [Supporting], Writing—review & editing [Supporting])

Data Availability

Raw sequences were deposited in the Sequence Read Archive (SRA: BioProject ID PRJNA982881). Concantenated 75% completeness alignment matrix, tree files, morphological data used for ASR, and complete phylogenetic signal statistics table are deposited in FigShare.