Spider phylogenomics: untangling the Spider Tree of Life

Sampling, extraction, assembly

Spider sequence data representing all major lineages were collected from previously published transcriptomic and genomic resources (N = 53) and supplemented with newly sequenced transcriptomes (N = 21) to form the target taxon set for the current study. Existing sequence data were acquired via the NCBI SRA database (http://www.ncbi.nlm.nih.gov/sra). Raw transcriptome sequences were downloaded, converted to fastq file format, and assembled using Trinity (Grabherr et al., 2011). Genomic data sets in the form of predicted proteins were downloaded directly from the literature (Sanggaard et al., 2014) for downstream use in our pipeline. Newly sequenced spiders were collected from a variety of sources, extracted using the TRIzol total RNA extraction method, purified with the RNeasy mini kit (Qiagen) and sequenced in-house at the Auburn University Core Genetics and Sequencing Laboratory using an Illumina Hi-Seq 2500. This produced 100bp paired end reads for each newly sequenced spider transcriptome, which were then assembled using Trinity. Proteins were predicted from each transcriptome using the program TransDecoder (Haas et al., 2013).

Core ortholog approach and data processing

We employed a core ortholog approach for putative ortholog selection and implicitly compared the effect of using a common arthropod core ortholog set and one compiled for spiders; the arthropod core ortholog set was deployed as described in Bond et al. (2014). To generate the spider core ortholog set, we used an all-versus-all BLASTP method (Altschul et al., 1990) to compare the transcripts of the amblypygid Damon variegatus, and the spiders Acanthoscurria geniculata, Dolomedes triton, Ero leonina, Hypochilus pococki, Leucauge venusta, Liphistius malayanus, Megahexura fulva, Neoscona arabesca, Stegodyphus mimosarum, and Uloborus sp.. Acanthoscurria geniculata and Stegodyphus mimosarum were represented by predicted transcripts from completely sequenced genomes while the other taxa were represented by our new Illumina transcriptomes. An e-value cut-off of 10–5 was used. Next, based on the BLASTP results, Markov clustering was conducted using OrthoMCL 2.0 (Li, Stoeckert & Roos, 2003) with an inflation parameter of 2.1.

The resulting putatively orthologous groups (OGs) were processed with a modified version of the bioinformatics pipeline employed by Kocot et al. (2011). First, sequences shorter than 100 amino acids in length were discarded. Next, each candidate OG was aligned with MAFFT (Katoh, 2005) using the automatic alignment strategy with a maxiterate value of 1,000. To screen OGs for evidence of paralogy, an “approximately maximum likelihood tree” was inferred for each remaining alignment using FastTree 2 (Price, Dehal & Arkin, 2010). Briefly, this program constructs an initial neighbor-joining tree and improves it using minimum evolution with nearest neighbor interchange (NNI) subtree rearrangement. FastTree subsequently uses minimum evolution with subtree pruning regrafting (SPR) and maximum likelihood using NNI to further improve the tree. We used the “slow” and “gamma” options; “slow” specifies a more exhaustive NNI search, while “gamma” reports the likelihood under a discrete gamma approximation with 20 categories, after the final round of optimizing branch lengths. PhyloTreePruner (Kocot, Citarella & Halanych, 2013) was then employed as a tree-based approach to screen each candidate OG for evidence of paralogy. First, nodes with support values below 0.95 were collapsed into polytomies. Next, the maximally inclusive subtree was selected where all taxa were represented by no more than one sequence or, in cases where more than one sequence was present for any taxon, all sequences from that taxon formed a monophyletic group or were part of the same polytomy. Putative paralogs (sequences falling outside of this maximally inclusive subtree) were then deleted from the input alignment. In cases where multiple sequences from the same taxon formed a clade or were part of the same polytomy, all sequences but the longest were deleted. Lastly, in order to eliminate orthology groups with poor taxon sampling, all groups sampled for fewer than 7 of the 11 taxa and all groups not sampled for Megahexura fulva (taxon with greatest number of identified OGs) were discarded. The remaining alignments were used to build profile hidden Markov models (pHMMs) for HaMStR with hmmbuild and hmmcalibrate from the HMMER package (Eddy, 2011).

For orthology inference, we employed HaMStR v13.2.3 (Ebersberger, Strauss & Von Haeseler, 2009), which infers orthology based on predefined sets of orthologs. Translated transcripts for all taxa were searched against the new set of 4,934 spider-specific pHMMs (available for download from the Dryad Data Repository) and an arthropod core ortholog set previously employed in Bond et al. (2014). In the spider core ortholog analysis, the genome-derived Acanthoscurria geniculata OGs were used as the reference protein set for reciprocal best hit scoring. Daphnia pulex was used as the reference species for putative ortholog detection in the arthropod core ortholog analysis. Orthologs sharing a core identification number were pooled together for all taxa and processed using a modified version of the pipeline used to generate the custom spider ortholog set. In both analyses, sequences shorter than 75 amino acids were deleted first. OGs sampled for fewer than 10 taxa were then discarded. Redundant identical sequences were removed with the perl script uniqhaplo.pl (available at http://raven.iab.alaska.edu/ ntakebay/) leaving only unique sequences for each taxon. Next, in cases where one of the first or last 20 characters of an amino acid sequence was an X (corresponding to a codon with an ambiguity, gap, or missing data), all characters between the X and that end of the sequence were deleted and treated as missing data. Each OG was then aligned with MAFFT (mafft—auto—localpair—maxiterate 1,000; Katoh (2005)). Alignments were then trimmed with ALISCORE (Misof & Misof, 2009) and ALICUT (Kück, 2009) to remove ambiguously aligned regions. Next, a consensus sequence was inferred for each alignment using the EMBOSS program infoalign (Rice, Longden & Bleasby, 2000). For each sequence in each single-gene amino acid alignment, the percentage of positions of that sequence that differed from the consensus of the alignment were calculated using infoalign’s “change” calculation. Any sequence with a “change” value greater than 75 was deleted. Subsequently, a custom script was used to delete any mistranslated sequence regions of 20 or fewer amino acids in length surrounded by ten or more gaps on either side. This step was important, as sequence ends were occasionally mistranslated or misaligned. Alignment columns with fewer than four non-gap characters were subsequently deleted. At this point, alignments shorter than 75 amino acids in length were discarded. Lastly, we deleted sequences that did not overlap with all other sequences in the alignment by at least 20 amino acids, starting with the shortest sequence not meeting this criterion. This step was necessary for downstream single-gene phylogenetic tree reconstruction. As a final filtering step, OGs sampled for fewer than 10 taxa were discarded.

In some cases, a taxon was represented in an OG by two or more sequences (splice variants, lineage-specific gene duplications (= inparalogs), overlooked paralogs, or exogenous contamination). In order to select the best sequence for each taxon and exclude any overlooked paralogs or exogenous contamination, we built trees in FastTree 2 (Price, Dehal & Arkin, 2010) and used PhyloTreePruner to select the best sequence for each taxon as described above. Remaining OGs were then concatenated using FASconCAT (Kück & Meusemann, 2010). The OGs selected by our bioinformatic pipeline were further screened in seven different ways (subsets listed in Table 2). OGs were first sorted based on amount of missing data; the half with the lowest levels was pulled out as matrix 2 (1,699 genes). From matrix 2, a smaller subset of OGs optimized for gene occupancy was extracted, resulting in matrix 3 (850 genes). The full supermatrix (matrix 1) was also optimized using the programs MARE (Meyer, Meusemann & Misof, 2011) and BaCoCa (Base Composition Calculator; Kück & Struck, 2014). MARE assesses the supermatrix by partition, providing a measure of tree-likeness for each gene and optimizes the supermatrix for information content. The full supermatrix was optimized with an alpha value of 5, to produce matrix 7 (1,488 genes, 58 taxa). From the MARE-reduced matrix, genes having no missing partitions for any of the remaining taxa (n = 50) were extracted to form a starting matrix for the BEAST analyses (details below). Matrix assessment was also conducted using BaCoCa, which provides a number of descriptive supermatrix statistics for evaluating bias in amino acid composition and patterns in missing data. This program was used to assess patterns of non-random clusters of sequences in the data, which could potentially mislead phylogenetic analyses. Matrix 4 represents a 50% reduction of the full supermatrix using BaCoCa derived values for phylogenetically informative sites as a guide; essentially reducing missing data from absent partitions and gaps. This matrix is similar, but not identical to matrix 2. Matrix 5 resulted from application of arthropod core OGs from Bond et al. (2014) to the extended taxon set. Matrix 6 represents the full spider core OG matrix (matrix 1) with Stegodyphus pruned from the tree. OGs for each matrix were concatenated using FASconCAT (Kück & Meusemann, 2010).

Phylogenetics

Table 2 summarizes run parameters of the seven individual maximum likelihood analyses conducted for each of the supermatrices. We selected the optimal tree for each supermatrix using the computer program ExaML ver. 3.0.1 (Kozlov, Aberer & Stamatakis, 2015). Models of amino acid substitution were selected using the AUTOF command in ExaML. Bootstrap data sets and starting parsimony trees for each matrix were generated using RAxML (Stamatakis, 2014) and each individually analyzed in ExaML. We generated 225–300 replicates for each matrix which were then used to construct a majority-rule bootstrap consensus tree; a custom python script was used to automate the process and write a bash script to execute the analyses on a high performance computing (HPC) cluster. The arthropod core OG bootstrap analysis was conducted using RAxML. All analyses were conducted on the Auburn University CASIC HPC and Atrax (Bond Lab, Auburn University).

A coalescent-based method as implemented in ASTRAL (Accurate Species TRee ALgorithm; Mirarab et al., 2014) was used to infer a species tree from a series of unrooted gene trees. The ASTRAL approach is thought to be more robust to incomplete lineage sorting, or deep coalescence, than maximum likelihood analysis of concatenated matrices and works quickly on genome-scale datasets (Mirarab et al., 2014). We first constructed individual gene trees for all partitions contained within matrix 1. Gene trees were generated using ML based on 100 RAxML random addition sequence replicates followed by 100 bootstrap replicates (Table 2). Subsequent species tree estimation was inferred using ASTRAL v4.7.6, from all individual unrooted gene trees (and bootstrap replicates), under the multi-species coalescent model.

A chronogram was inferred in a Bayesian framework under an uncorrelated lognormal relaxed clock model (Drummond et al., 2006; Drummond, 2007) using Beast v1.8.1 (Drummond et al., 2012). For this analysis we used 43 partitions of a matrix which included complete partitions for all taxa derived from the MARE-optimized matrix 7. The model of protein evolution for each partition was determined using the perl script ProteinModelSelection.pl in RAxML. BEAST analyses were run separately for each partition using eight calibration points based on fossil data. The most recent common ancestor (MRCA) of Mesothelae + all remaining spiders was given a lognormal prior of (mean in real space) 349 Ma (SD = 0.1) based on the Mesothelae fossil Palaeothele montceauensis (Selden, 1996). The MRCA of extant araneomorphs was given a lognormal prior of (mean in real space) 267 Ma (SD = 0.2) based on the fossil Triassaraneus andersonorum (Selden et al., 1999). The MRCA of extant mygalomorphs was given a lognormal prior of (mean in real space) 278 Ma (SD = 0.1) based on the fossil Rosamygale grauvogeli (Selden & Gall, 1992). The MRCA of Haplogynae + Hypochilidae was given a lognormal prior of (mean in real space) 278 Ma (SD = 0.1) based on the fossil Eoplectreurys gertschi(Selden & Penney, 2010). The MRCA of Deinopoidea (cribellate orb-weavers) was given a lognormal prior of (mean in real space) 195 Ma (SD = 0.3) based on the fossil Mongolarachne jurassica (Selden, Shih & Ren, 2013). The MRCA of ecribellate orb-weavers was given a lognormal prior of (mean in real space) 168 Ma (SD = 0.4) based on the fossil Mesozygiella dunlopi (Penney & Ortu, 2006). The MRCA of Nemesiidae, excluding Damarchus, was given a lognormal prior of (mean in real space) 168 Ma (SD = 0.4) based on the nemesiid fossil Cretamygale chasei (Selden, 2002). Finally, the MRCA of Antrodiaetidae was given a lognormal prior of (mean in real space) 168 Ma (SD = 0.4) based on the fossil Cretacattyma raveni (Eskov & Zonstein, 1990). Two or more independent Markov Chain Monte Carlo (MCMC) searches were performed until a parameter effective sample size (ESS) >200 was achieved. ESS values were examined in Tracer v1.5. Independent runs for each partition were assembled with LogCombiner v1.7.5 and 10% percent of generations were discarded as burn-in. Tree files for each partition where then uniformly sampled to obtain 10,000 trees. A total of 430,000 trees (10,000 trees from each partition) were assembled with LogCombiner v1.7.5 and a consensus tree was produced using TreeAnnotator v1.8.1. A chronogram containing all taxa was generated using a penalized likelihood method in r8s v1.8 (Sanderson, 2002). The 95% highest posterior density dates obtained for the BEAST analysis were incorporated as constraints for node ages of the eight fossil calibrated nodes. The analysis was performed using the TN algorithm, cross validation of branch-length variation and rate variation modeled as a gamma distribution with an alpha shape parameter.

To detect diversification rate shifts, we performed a Bayesian analysis of diversification in BAMM (Bayesian Analysis of Macroevolutionary Mixtures; Rabosky et al., 2014). For this analysis we used the chronogram obtained by the r8s analysis in order to maximize taxon sampling. To account for non-random missing speciation events, we quantified the percentage of taxa sampled per family (World Spider Catalog, 2015) and incorporated these into the analysis. We also accounted for missing families sampled. The MCMC chain was run for 100,000,000 generations, with sampling every 10,000 generations. Convergence diagnostics were examined using coda (Plummer et al., 2006) in R. Ten percent of the runs were discarded as burn-in. The 95% credible set of shift configurations was plotted in the R package BAMMtools (Rabosky et al., 2014).

Character state reconstructions of web type following Blackledge et al. (2009) were performed using a maximum likelihood approach. The ML approach was implemented using the rayDISC command in the package corHMM (Beaulieu, O’Meara & Donoghue, 2013) in R (Ihaka & Gentleman, 1996). This method allows for multistate characters, unresolved nodes, and ambiguities (polymorphic taxa or missing data). Three models of character evolution were evaluated under the ML method: equal rates (ER), symmetrical (SYM) and all rates different (ARD). A likelihood-ratio test was performed to select among these varying models of character evolution.

Summary of genomic data

Twenty-one novel spider transcriptomes were sequenced, with an average of 72,487 assembled contigs (contiguous sequences) ranging from 6,816 (Diguetia sp.) to 191,839 (Segestria sp.); specimen data and transcriptome statistics for each sample are summarized in Tables S1 and S2 respectively. Median contig length for the novel transcriptomes was 612 bp. The complete taxon set, including spider and outgroup transcriptomes from the SRA database, had an average contig number of 53,740 and a range of 5,158 (Paratropis sp.) to 202,311 (Amaurobius ferox) with a median contig length of 655. The newly constructed spider-specific core ortholog group (OG) set contained 4,934 OGs, more than three times the number of arthropod core orthologs used in prior spider analyses (Bond et al., 2014) and represents a significant step forward in generating a pool of reasonably well-vetted orthologs for spider phylogenomic analyses. The arthropod and spider core orthology sets had 749 groups in common; 4,185 OGs in the spider core were novel. Of the spider-core groups, 4,249 (86%) were present in the sequenced genome of our HaMSTR reference taxon of choice Acanthoscurria geniculata (Sanggaard et al., 2014) and were retained for use in downstream ortholog detection. The number of TransDecoder predicted proteins and ortholog detection success for each taxon is summarized in Table S2. Annotations for the arthropod set can be found in Bond et al. (2014); Table S3 summarizes gene annotations for the spider core ortholog set generated for this study. Our new HaMStR spider core ortholog set and Acanthoscurria geniculata BLAST database file can be downloaded from the Dryad Data Repository at doi:10.5061/dryad.6p072.

Phylogenetic analyses

Seven super matrices were generated for downstream non time-calibrated analyses (Fig. 2), one drawn from the arthropod core set and six using the spider core set. Data set sizes, summarized in Table 2, ranged from a maximum of 3,398 OGs with a higher percentage of missing cells (38.5%), 850 OGs with 19.6% missing, to 549 OGs (arthropod core set) with 33% missing data. Two matrices were generated using automated filtering approaches implemented by BaCoCa (Kück & Struck, 2014) and MARE (Meyer, Meusemann & Misof, 2011). In BaCoCa we sorted partitions using number of informative sites, capturing the top half (1,700 OGs) of the matrix containing the most informative sites. RCFV values generated by BaCoCa were <0.05 for all taxa in all partitions for each of the matrices, indicating homogeneity in base composition. Additionally, there was no perceptible taxonomic bias observed in shared missing data (Figs. S1–S6). The MARE optimized matrix comprised 58 taxa and 1,488 genes with 19.6% missing data. For graphical representations of gene occupancy for each matrix, see Figs. S7–S12. Blast2GO (Conesa et al., 2005) gene ontology distributions of molecular function for OGs recovered from both the spider and arthropod ortholog sets (Figs. S13 and S14) can be found in the Supplemental Information.

Our phylogenetic analyses (see Table 2 and ‘Discussion’), the results of which are summarized in Fig. 2, consistently recover many well-supported monophyletic groups: Araneae, Mygalomorphae, Araneomorphae, Synspermiata (i.e., Haplogynae excluding Filistatidae and Leptonetidae), Entelegynae, the RTA clade, Dionycha, and Lycosoidea. Within Mygalomorphae, Atypoidina and Avicularioidea are monophyletic; Nemesiidae is polyphyletic. Filistatidae (Kukulcania) emerges as the sister group to Hypochilus. Interestingly, Leptonetidae emerges as the sister group to Entelegynae. Eresidae is sister to Araneoidea, similar to findings of Miller et al. (2010). Deinopoidea is polyphyletic. Oecobiidae is sister to Uloboridae, which are together sister to Deinopidae plus the RTA clade. Homalonychidae and by implication the entire Zodarioidea (Miller et al., 2010), is sister to Dionycha plus Lycosoidea. Hahniidae, represented by the cryphoecine Calymmaria, is sister to Dictynidae. Thomisidae belongs in Lycosoidea as proposed by Homann (1971) and Polotow, Carmichael & Griswold (2015) (see also Ramírez, 2014).

Coalescent-based species-tree analysis in ASTRAL employed unrooted gene trees based on the 3,398 gene matrix as input and inferred a well-supported tree (most nodes >95% bs; Fig. 3). With few exceptions the topology recovered using this approach was congruent with the likelihood-based supermatrix analysis. Conflicting nodes, some corresponding to key araneomorph lineages, which were moderately to weakly supported in concatenated analyses, are summarized in Fig. 2.

A chronogram based on 43 partitions with no missing data (matrix 7, see Table 2) is shown in Fig. 4. MRCA divergence time estimates are summarized in Table 3: Mesothelae—Opisthothelae at 340 Ma (95% CI[287–398]); Mygalomorphae—Araneomorphae at 308 Ma (95% CI[258–365]); Synspermiata + Hypochilidae—Entelegynae at 276 Ma (95% CI[223–330]); RTA + Deinopoidea—Stegodyphus + Araneoidea at 214 Ma (95% CI[154–280]); RTA—Dionycha at 138.8 Ma (Fig. 4).

Diversification rate shift analysis estimated three instances of significant diversification shifts within spiders (95% credibility). The highest rate shift is within the RTA + Dionycha + Lycosoidea (Fig. 5) followed by Avicularioidea and within Araneoidea (f = 0.23; 0.21; Fig. 5).

Maximum likelihood ancestral state reconstruction of web type (Fig. 6) shows that the spider common ancestor likely foraged from a subterranean burrow, sometimes sealed by a trapdoor. The ancestral condition for araneomorphs may have been a stereotypical aerial sheet. Entelegynae ancestors probably spun orbs, which were subsequently lost at least three times. RTA taxa largely abandoned webs to become hunting spiders. Precise location of these character state shifts depends upon sufficient sampling; denser sampling reduces the number of unobserved evolutionary events. While this analysis contains only 47 of 114 spider families, the sequence and overall mapping to the spider backbone phylogeny is strongly supported.

Data characteristics and development of spider core orthologs

Transcriptome analyses are unquestionably data rich. Thousands of assembled sequences emerge from even modest RNA-seq experiments, providing, among other things, a basis for identifying phylogenetically informative orthologs. This bounty comes with a few caveats. Isoforms, paralogous sequences, and assembly artifacts (chimeric contigs) can mislead inference of single-copy orthologous genes. The data represent one snapshot—a specific organism, point in time, and combination of tissues—that can lead to gaps in downstream supermatrices due to stochastic sampling issues. Large amounts of missing data, due to missing loci and indels introduced during alignment, can arise post-assembly in the ortholog detection and filtering stages of phylogenomic analyses (Bond et al., 2014; Fernández, Hormiga & Giribet, 2014). Lemmon et al. (2009) and a number of other authors (Roure, Baurain & Philippe, 2013; Dell’Ampio et al., 2014; Xia, 2014) have discussed the potential negative effects of such missing data in large phylogenomic (transcriptome-based) datasets. Recent studies argue that the phylogenetic signal from transcriptomes can conflict with alternative reduced representation approaches like targeted sequence capture (Jarvis et al., 2014; Brandley et al., 2015; Prum et al., 2015). From vast amounts of bird genome protein-coding data, Jarvis et al. (2014) concluded that these loci were not only insufficient (low support values), but also misleading due to convergence and high levels of incomplete lineage sorting during rapid radiations.

Simulation studies now predict that 10s–100s of loci will resolve most phylogenies, albeit this is sensitive to factors such as population size or speciation tempos (Knowles & Kubatko, 2011; Leache & Rannala, 2011; Liu & Yu, 2011). To mitigate the impacts of paralogy, incomplete lineage sorting, and missing data, we developed a priori a set of spider core orthologs that comprise a database consisting of over 4,500 genes that are expected to be recovered from most whole spider RNA extractions and are likely orthologous. We summarize the annotations for each of the genes in the HaMStR pHMM file in Table S3.

Our approach enhances repeatability, downstream assessment, scalability (taxon addition), and data quality. Studies that employ pure clustering approaches like OMA stand-alone (Altenhoff et al., 2013) may produce more data (i.e., more “genes”) on the front end; however, they present some problems in terms of ease of scalability. Although adding more genes is one strategy, it is increasingly clear that taxon sampling and data quality are very important (Lemmon & Lemmon, 2013; Bond et al., 2014).

A modified view of spider evolution and key innovations

Once considered the “crowning achievement of aerial spiders” (Gertsch, 1979), the orb web and consequent adaptive radiation of araneoid spiders (ecribellate orb weavers and their relatives) has captured the imagination of spider researchers for over a century. The evolution of adhesive threads and the vertical orientation of the orb web, positioned to intercept and retain flying insects, has been long considered a “key innovation” that allowed spiders to inhabit a new adaptive zone (Bond & Opell, 1998). It is important to note that several prior authors speculated about orb web adaptive value, such as Levi (1980), Opell (1979), Opell (1982) and Coddington (1986) although Bond & Opell (1998) quantified the pattern in a formal phylogenetic framework. Over 25% of all described spider species are araneoids. Given orb weaver monophyly on quantitative phylogenies (Griswold et al., 1998; Blackledge et al., 2009), rigorous empirical studies tended to confirm the orb as a prime cause of spider diversification (Bond & Opell, 1998). Nevertheless, a lack of correlation of the orb web and species richness has been apparent for some time. Griswold et al. (1998) noted that over 50% of Araneoidea no longer build recognizable orb webs and suggested that “the orb web has been an evolutionary base camp rather than a summit.”

Bond et al. (2014) tested two alternative evolutionary scenarios for orb web evolution, reflecting different analytical results; parsimony implied multiple independent origins, and maximum likelihood implied one origin and subsequent multiple losses. The current study (Fig. 6) favors the latter: the orb evolves at the base of the araneoid + deinopoid + RTA clade, but is lost at least three times independently. Large amounts of morphological and behavioral data (albeit often correlated with features essential to the orb) still support the single origin hypothesis (Coddington, 1986; Coddington, 1991; Scharff & Coddington, 1997; Griswold et al., 1998; Agnarsson, Coddington & Kuntner, 2013). Our results suggest both that the orb web originated earlier than previously supposed, and that heretofore-unsuspected clades of spiders descend from orb weavers. In a sense, this ancient origin hypothesis reconciles the implications of genomic data with the classical evidence for multiple, homologous, complex, co-adapted character systems.

Recent discoveries of large, cribellate orb web-weaving taxa from the late Triassic agree with our molecular dates. Diverse Mesozoic deinopoids (Selden, Ren & Shih, 2015) are consistent with the “orb web node” at 213 Ma (Fig. 4 and Table 3). Under this view, modern uloborids and deinopids are distinct remnants of this diverse group. Selden, Ren & Shih (2015) previously noted that if other extant taxa “emerged from the deinopoid stem or crown group it would render the whole-group Deinopoidea paraphyletic”; we discuss this scenario in detail below.

Contrary to the contemporary paradigm that the evolution of the orb web and adhesive sticky threads elevated rates of diversification among the araneoid spiders, our BAMM analysis (Fig. 5) indicates that the highest rates of diversification likely occurred among the RTA spiders followed by mygalomorphs and then araneoids as a distant third, the latter driven–in part–by the secondarily non-orb weaving theridiids and linyphiids. These results imply that other foraging strategies (e.g., cursorial hunting and irregular sheets) were a more “successful” strategy than the orb. Indeed, the point estimate for the RTA node during the early Cretaceous (138.8 Ma; Fig. 4 and Table 3) precedes the subsequent diversification of the RTA clade at 125–100 Ma.

This date coincides with the Cretaceous Terrestrial Revolution (KTR). Angiosperms radiated extensively at 125–90 Ma (Crane & Friis, 1987; Wang, Zhang & Jarzembowski, 2013), as did various plant-dependent insect lineages, including beetles (McKenna et al., 2009; Mckenna et al., 2015), lepidopterans (Wahlberg, Wheat & Peña, 2013), ants (Moreau et al., 2006) and holometabolous insects in general (Misof et al., 2014), although some insect lineages do not show a pulse (e.g., darkling beetles; Kergoat et al., 2014). Spiders, as important insect predators, may also have diversified rapidly along with their prey (e.g., Penney, Wheater & Selden, 2003; Peñalver, 2006; Selden & Penney, 2010). The fossil and phylogenomic data presented here show that most spider lineages predate the KTR (Selden & Penney, 2010; Bond et al., 2014). Among these, the RTA clade especially, but also mygalomorphs and araneoids, diversified in response to the KTR insect pulse. That aerial web spinners specialized on rapidly radiating clades of flying insects is hardly surprising. Similarly, if forest litter habitats became more complex and spurred insect diversification (Moreau et al., 2006) ground-dwelling spiders may also have diversified at unusual rates. Perhaps the most dramatic change in insect abundances occurred with the origin and early diversification of social insects (Hölldobler & Wilson, 1990) and beetles (Mckenna et al., 2015). Both groups date back to 150–125 Ma and diversified during the KTR (LaPolla, Dlussky & Perrichot, 2013; Ward, 2014; Legendre et al., 2015). A major increase in these insect groups may have favoured spiders that feed on cursorial prey and thus could help explain the concurrent increase in diversification in the RTA clade, mygalomorphs, and non-orb weaving araneoids such as cobweb weavers (Dziki et al., 2015).

Taken together, this new evidence on character evolution, divergence estimates, and rates of diversification indicates that previous conclusions regarding the timing and rate of spider evolution were imprecise. Our data support an ancient orb web hypothesis that is further bolstered by a wealth of fossil data showing that a cribellate deinopoid stem group likely diversified during the early Mesozoic. Molecular divergence clock estimates are consistent with the placement of the orb web further down the tree as well as suggesting that some of the greatest rates of species diversification coincided with the KTR. The latter suggests that spiders took advantage of increased abundance of cursorial prey.

These findings likely diminish the hypothesis proposed by Bond & Opell (1998) that the vertically oriented orb web represented a key innovation, particularly in light of the fact that over half of araneoid species do not build an orb web (e.g. Theridiidae and Linyphiidae; noted by Griswold et al., 1998; Fernández, Hormiga & Giribet, 2014). We already knew that major orb web-weaving groups are very successful in spite of abandoning the orb (Blackledge et al., 2009).

Spider systematics

Although our results show that many classical ideas in spider systematics require revision (e.g., mygalomorph families, Haplogynae, paleocribellates, higher araneoids, and RTA + dionychan lineages), they also robustly support many classical taxonomic concepts.