Asymmetrical rates of cladogenesis and extinction abound in the Tree of Life, resulting in numerous minute clades that are dwarfed by larger sister groups. Such taxa are commonly regarded as phylogenetic relicts or “living fossils” when they exhibit an ancient first appearance in the fossil record and prolonged external morphological stasis, particularly in comparison to their more diversified sister groups. Due to their special status, various phylogenetic relicts tend to be well-studied and prioritized for conservation. A notable exception to this trend is found within Amblypygi (“whip spiders”), a visually striking order of functionally hexapodous arachnids that are notable for their antenniform first walking leg pair (the eponymous “whips”). Paleoamblypygi, the putative sister group to the remaining Amblypygi, is known from Late Carboniferous and Eocene deposits, but is survived by a single living species, Paracharon caecus Hansen, 1921, that was last collected in 1899. Due to the absence of genomic sequence-grade tissue for this vital taxon, there is no global molecular phylogeny for Amblypygi to date, nor a fossil-calibrated estimation of divergences within the group. Here, we report several individuals of a previously unknown species of Paleoamblypygi from a cave site in Colombia. Capitalizing upon this discovery, we generated the first molecular phylogeny of Amblypygi, integrating ultraconserved element sequencing with legacy Sanger datasets and including described extant genera. To quantify the impact of sampling Paleoamblypygi on divergence time estimation, we performed in silico experiments with pruning of Paracharon. We demonstrate that the omission of relicts has a significant impact on the accuracy of node dating approaches that outweighs the impact of excluding ingroup fossils. Our results underscore the imperative for biodiversity discovery efforts in elucidating the phylogenetic relationships of “ dark taxa”, and especially phylogenetic relicts in tropical and subtropical habitats.

Species sampling

Specimens sequenced for this study were hand collected from field sites or contributed by collections of the National Museum of Natural History, the Smithsonian Institution, Washington, DC, USA (USNM); the Museum of Comparative Zoology, USA (MCZ); the Natural History Museum of Denmark, Copenhagen (NHMD); the Museu Nacional do Rio de Janeiro, Brazil (MNRJ); the California Academy of Sciences, San Francisco, USA (CAS); the State Museum of Natural History Stuttgart, Germany (SNMS); and the Universidade Federal do Piauí, Brazil (CHNUFPI). We selected exemplars of each of the five described extant families, and all 17 described genera: 20 specimens of Charinidae (133 described species), one Charontidae (15 described species), 1 Paracharontidae (1 described species), two Phrynichidae (35 described species), and ten Phrynidae (77 described species). Outgroup sampling leveraged previous phylogenomic works that have established Amblypygi as the sister group of vinegaroons (Thelyphonida), and short-tailed whip scorpions (Schizomida) (the clade Pedipalpi; Sharma et al. 2014; Ballesteros and Sharma 2019; Ballesteros et al. 2019). Pedipalpi in turn is understood to form a clade with spiders, scorpions, and pseudoscorpions, a relationship supported by rare genomic changes (e.g., Ontano et al. 2021, 2022; Ballesteros et al. 2022). We therefore included two representatives of spiders, one vinegaroon, three Schizomida, and rooted the tree with three scorpions. Accession data for all specimens in this study are provided in Supplementary Table S2.

Ultraconserved element sequencing

Ultraconserved element (UCE) sequences generated in this study were augmented with published UCE and RNASeq datasets (Supplementary Table S3). We analyzed 47 terminals with UCE data including 38 Amblypygi. For newly sequenced specimens, one leg was used for DNA extractions from one specimen using the DNeasy™ Tissue Kit (Qiagen Inc., Valencia, CA). Libraries were prepared and enriched following protocols in Faircloth et al. (2015), but following the modifications detailed in the Supplementary Material. All pools were enriched with the Arachnida probe set (Starrett et al. 2017) except Heterophrynus and Charon, which were enriched using the Spider2Kv1 probe set (Kulkarni et al. 2020) following the myBaits protocol 4.01 (Arbor Biosciences).

Inclusion of Sanger-sequenced terminals

While UCE datasets have demonstrated great potential for leveraging historical collections, the compatibility of these data with legacy datasets in Sanger-sequencing studies is often not evaluated. To integrate historical datasets in a comprehensive phylogenetic framework, we compiled another matrix comprised of 109 terminals for six publicly available Sanger-sequenced loci: 12S ribosomal RNA (12S), 16S rRNA (16S), cytochrome c oxidase subunit 1 (COI), histone H3 (H3), and the small and large subunits of nuclear ribosomal genes (18S and 28S, respectively). We integrated this data set by querying the raw files of the UCE data set for any potential match with 18S rRNA and 28S rRNA. Finally, COI and H3 markers were aligned with peptide translation using MACSE (Ranwez et al. 2011), with the invertebrate mitochondrial code implemented for COI. The remaining markers (12S, 16S, 18S and 28S) were aligned using MAFFT version 7 (Katoh and Standley 2013). Trimming was performed on all alignments using trimAL (Capella-Gutiérrez et al. 2009) with -gappyout.

Phylogenomic analyses

The assembly, alignment, trimming and concatenation of data were done using the PHYLUCE pipeline (publicly available at https://phyluce.readthedocs.io/en/latest/). UCE contigs were extracted using the Spider2Kv1 probe set (Kulkarni et al. 2020) and the Blended probe set (Maddison et al. 2020). We applied gene occupancies of 10%, 25% and 40% on the UCE data set. Additionally, we also analyzed 1% occupancy of the UCE data set to allow inclusion of all loci in the reconstruction of the phylogeny. We screened for orthologous and duplicate loci with the minimum identity and coverage of 65 and 65 matches.

To augment the UCE dataset with RNASeq datasets, we followed the assembly, sanitation and reading frame detection pipeline as in Fernández et al. (2018) for assembling the transcriptomes. Additionally, we ran the Perl script for Rcorrector (Song and Florea, 2015) for error correction and downstream efficiency prior to assembly. FASTA files of transcriptomes resulting from CD-HIT-EST were converted to 2bit format using faToTwoBit, (Kent et al. 2002). In the PHYLUCE environment (publicly available at https://phyluce.readthedocs.io/en/latest/tutorial-three.html), we created a temporary relational database to summarize probe to assembly match using: phyluce_probe_run_multiple_lastzs_sqlite function on the 2bit files. The ultraconserved loci were recovered by the phyluce_probe_slice_sequence_from_genomes command. The resulting FASTA files were treated as contigs and used to match the reads to the Blended probe set of Maddison et al. (2020).

Phylogenetic analyses were performed on two types of data sets: an unpartitioned UCE nucleotide data set alone (at difference thresholds of occupancy); and a matrix of unpartitioned UCE paired with the partitioned six Sanger loci matrix. Maximum likelihood analyses were performed using IQ-TREE (Nguyen et al. 2015) version 2. Model selection was allowed for each unpartitioned dataset using the TEST function (Kalyaanamoorthy et al. 2018, Hoang et al. 2018). Nodal support was estimated via 1,000 ultrafast bootstrap (UFBoot) replicates (Hoang et al. 2018) and Shimodaira-Hasegawa approximate likelihood ratio testing (SH-aLRT) with 1,000 iterations. We used the flag -bnni which reduces the risk of overestimating branch support with UFBoot due to model violations. This flag optimizes each bootstrap tree using a hill-climbing nearest neighbor interchange (NNI) search based on the corresponding bootstrap alignment (Hoang et al. 2018). We used gene (gCF) and site concordance factors (sCF) to evaluate the percentage of gene trees and decisive alignment sites containing a given branch in the maximum likelihood tree implemented in IQ-TREE (Minh et al. 2020).

To evaluate topological robustness of selected nodes, we evaluated signal for alternative placements using quartet likelihood mapping (Strimmer and Von Haeseler 1997) in IQ-TREE. Likelihood mapping was performed against the complete dataset, representing all UCE loci (i.e., not filtered for any occupancy threshold).

Phylogenomic dating

Divergence time estimation was performed using a Bayesian inference approach, as implemented in codeml and MCMCTree (both part of PAML v.4.8; Yang 2007; dos Reis and Yang 2019). We optimized the fossil information-based calibrations on a matrix comprising smallest UCE matrix (selected for highest taxon occupancy), plus the six Sanger loci. The maximum likelihood tree topology inferred for this dataset served as the basis for divergence time estimation, implementing uniform node age priors to accommodate the scarcity of terrestrial chelicerate fossils. The root age was set to a range of 545 Mya to 435 Mya, based on the age of Eramoscorpius brucensis, the oldest known arachnopulmonate. Weygoldtina anglica, the oldest known stem-Palaeoamblypygi, was used to constrain the most recent common ancestor (MRCA) of crown Amblypygi. The Cretaceous fossils Kronocharon prendinii and Britopygus weygoldtii were assigned to the MRCA of (Charinidae + Charontidae) and (Phrynidae + Phrynichidae), respectively. Four other outgroup nodes were calibrated using the oldest unambiguous fossils representing those clades. Justifications and references for node calibrations are provided in Supplementary Table S3. Both chains were run for 20,000 generations, with an additional 5% set for burnin. We used the independent rates clock model for all partitions and all analyses were run using the same seed value (arbitrarily set to 5) to make the results reproducible.

Separately, we inferred divergence times under a penalized likelihood approach, as implemented in LSD2, which uses a least-squares approach based on a Gaussian model and is robust to uncorrelated violations of the molecular clock (To et al. 2016). The root age was set to a maximum of 545 Ma (Supplementary Table S3). We used the commands prime and thorough to optimize the analyses, and cross validation was used to select the optimal smoothing parameter. Following Eberle et al. (2018), penalized likelihood optimization iterations were increased from the default of 2 to 5, and the number of penalized likelihood simulated annealing was doubled from 5,000 to 10,000. LSD2 requires at least one fixed date, so we used 314.6 Mya as a fixed date for the most recent common ancestor of Scorpiones.

To investigate the influence of the relict Paracharon on divergence time estimation (i.e., simulate the effect of not having discovered the Colombian specimens), we performed a separate family of analyses for both dating approaches (MCMCTree and LSD2), wherein we removed the terminal Paracharon sp.; the MRCA of the remaining Amblypygi (Euamblypygi) was therefore calibrated with the minimum age constraint reflecting the age of the Carboniferous fossil Weygoldtina anglica. We compared the median age and variance of nodes between runs to assess which parts of the chronogram were most prone to node age overestimation. In addition, to compare the significance of sampling Paracharon versus ingroup fossil calibrations, we performed a third family of analyses wherein we retained all terminals, but removed the three ingroup fossil node calibrations and keeping only the five calibration points for outgroups.