DOI: 10.5061/dryad.m63xsj449
title: The rediscovery of a relict unlocks the first global phylogeny of whip spiders (Amblypygi)

description: 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.
description: <em>Species
 sampling</em> 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. <em>Ultraconserved element
 sequencing</em> 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
 <em>Heterophrynus</em> and <em>Charon</em>, which
 were enriched using the Spider2Kv1 probe set (Kulkarni et al. 2020)
 following the myBaits protocol 4.01 (Arbor Biosciences). <em>Inclusion of Sanger-sequenced
 terminals</em> 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 -<em>gappyout</em>. <em>Phylogenomic
 analyses</em> 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: <em>phyluce_probe_run_multiple_lastzs_sqlite function</em> on the 2bit files. The ultraconserved loci were recovered by the <em>phyluce_probe_slice_sequence_from_genomes</em> 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). <em>Phylogenomic dating</em> 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 <em>Eramoscorpius brucensis</em>, the oldest known arachnopulmonate. <em>Weygoldtina anglica</em>, the oldest known stem-Palaeoamblypygi, was used to constrain the most recent common ancestor (MRCA) of crown Amblypygi. The Cretaceous fossils <em>Kronocharon prendinii</em> and <em>Britopygus weygoldtii </em>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 <em>prime</em> and <em>thorough</em> 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 <em>Paracharon</em> 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 <em>Paracharon</em> sp.; the MRCA of the remaining Amblypygi (Euamblypygi) was therefore calibrated with the minimum age constraint reflecting the age of the Carboniferous fossil <em>Weygoldtina anglica</em>. 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 <em>Paracharon</em> 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.