GigaScience, 8, 2019, 1–9 doi: 10.1093/gigascience/giz099 Data Note DATA NOTE The draft genome sequence of the spider Dysdera silvatica (Araneae, Dysderidae): A valuable resource for functional and evolutionary genomic studies in chelicerates Jose Francisco Sánchez-Herrero 1 , Cristina Frı́as-López 1 , Paula Escuer 1 , Silvia Hinojosa-Alvarez 1,2 , Miquel A. Arnedo 3 , Alejandro Sánchez-Gracia 1,* and Julio Rozas 1,* 1 Departament de Genètica, Microbiologia i Estadı́stica, Universitat de Barcelona (UB) and Institut de Recerca de la Biodiversitat (IRBio), Diagonal 643, 08028 Barcelona, Spain ; 2 Jardı́n Botánico, Instituto de Biologı́a, Universidad Nacional Autónoma de México, Tercer Circuito Exterior S/N, Ciudad Universitaria Coyoacán, 04510 México DF, México and 3 Departament de Biologia Evolutiva, Ecologia i Ciències Ambientals, Universitat de Barcelona (UB) and Institut de Recerca de la Biodiversitat (IRBio), Diagonal 643, 08028 Barcelona, Spain ∗ Correspondence address. Julio Rozas, Departament de Genètica, Microbiologia i Estadı́stica, Universitat de Barcelona (UB) and Institut de Recerca de la Biodiversitat (IRBio), Diagonal 643, 08028 Barcelona, Spain. E-mail: jrozas@ub.edu http://orcid.org/0000-0003-4543-4577; Alejandro Sánchez-Gracia, Departament de Genètica, Microbiologia i Estadı́stica, Universitat de Barcelona (UB) and Institut de Recerca de la Biodiversitat (IRBio), Diagonal 643, 08028 Barcelona, Spain; . E-mail: elsanchez@ub.edu http://orcid.org/0000-0002-6839-9148 Abstract Background: We present the draft genome sequence of Dysdera silvatica, a nocturnal ground-dwelling spider from a genus that has undergone a remarkable adaptive radiation in the Canary Islands. Results: The draft assembly was obtained using short (Illumina) and long (PaciBio and Nanopore) sequencing reads. Our de novo assembly (1.36 Gb), which represents 80% of the genome size estimated by flow cytometry (1.7 Gb), is constituted by a high fraction of interspersed repetitive elements (53.8%). The assembly completeness, using BUSCO and core eukaryotic genes, ranges from 90% to 96%. Functional annotations based on both ab initio and evidence-based information (including D. silvatica RNA sequencing) yielded a total of 48,619 protein-coding sequences, of which 36,398 (74.9%) have the molecular hallmark of known protein domains, or sequence similarity with Swiss-Prot sequences. The D. silvatica assembly is the first representative of the superfamily Dysderoidea, and just the second available genome of Synspermiata, one of the major evolutionary lineages of the “true spiders” (Araneomorphae). Conclusions: Dysderoids, which are known for their numerous instances of adaptation to underground environments, include some of the few examples of trophic specialization within spiders and are excellent models for the study of cryptic female choice. This resource will be therefore useful as a starting point to study fundamental evolutionary and functional questions, including the molecular bases of the adaptation to extreme environments and ecological shifts, as well of the origin and evolution of relevant spider traits, such as the venom and silk. Received: 6 May 2019; Revised: 27 June 2019; Accepted: 30 July 2019  C The Author(s) 2019. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 1 2 The draft genome sequence of the spider Dysdera silvatica (Araneae, Dysderidae) Keywords: Araneomorphae; hybrid genome assembly; genome annotation; Canary Islands recent review on arachnid genomics identified the superfam- ily Dysderoidea (namely, Dysderidae, Orsolobidae, Oonopidae, and Segestriidae) as one of the priority candidates for genome sequencing [14]. The new genome, intended to be a reference genome for genomic studies on trophic specialization, will also be a valuable source for the ongoing studies on the molecu- lar components of the chemosensory system in chelicerates [15]. Besides, because of the numerous instances of indepen- dent adaptation to caves [16], the peculiar holocentric chromo- somes [17], and the evidence for cryptic female choice mecha- nisms [18,19] within the family, the new genome will be a use- ful reference for the study of the molecular basis of adaptation to extreme environments, karyotype evolution, and sexual se- lection. Additionally, a new fully annotated spider genome will greatly improve our understanding of key features, such as the venom and silk. The availability of new genomic information Figure 1 Male of Dysdera silvatica from Teselinde (La Gomera, Canary Islands). in a sparsely sampled section of the tree of life of spiders [14] Photo credit: Miquel Arnedo. will further provide valuable knowledge about relevant scientific questions, such as gene content evolution across main arthro- Data Description pod groups, including the consequences of whole-genome du- plications, or the phylogenetic relationships with Araneae. Spiders are a highly diverse and abundant group of predatory arthropods, found in virtually all terrestrial ecosystems. Approx- imately 45,000 spider species have been recorded to date [1]. Sampling and DNA extraction The nocturnal ground family Dysderidae ranks 17th out of 118 currently accepted spider families in number of species. The We sampled adult individuals of D. silvatica in different locali- type genus of the family, Dysdera Latreille, 1804, includes half ties of La Gomera (Canary Islands) in March 2012 and June 2013 of the family diversity (282 species). This genus is remarkable (Supplementary Table S1-1). The species was confirmed in the in several aspects. First, it represents one of the few cases of laboratory, and samples were stored at −80o C until its use. For stenophagy, i.e., prey specialization, across spiders [2]. Many Illumina and PaciBio libraries (see below), we extracted genomic species in the genus have evolved special morphological, behav- DNA using Qiagen DNeasy Blood & Tissue Kit (Qiagen, Hilden, ioral, and physiological adaptations to feed on woodlice, includ- Germany, 74104) ) according to the manufacturer’s protocol. For ing modifications of mouthparts, unique hunting strategies, and the Oxford Nanopore libraries, we used a modified version of effective restriction to assimilation of metals into its tissues [3– the Blood & Cell Culture DNA Mini Kit (Qiagen). Due to the high 7]. Because of their chemical defenses and ability to accumulate amount of chitin present in spiders we incubated fresh origi- heavy metals from the soil, woodlice are usually avoided as prey nal samples 48 h at 32o C, avoiding a centrifugation step prior by most spiders, including generalist Dysdera [2,4, 5,7]. Although to sample loading to Qiagen Genomic tips, permitting the so- mostly circumscribed to the Mediterranean region, Dysdera has lution to precipitate by gravity. We also added an extra wash colonized all the Macaronesian archipelagoes and has under- with 70% ethanol and centrifuged the solution at >5,000g for gone a remarkable species diversification in the Canary Islands 10 min at 4o C. We quantified the genomic DNA in a Qubit flu- [8]. As many as 55 species have been recorded across the 7 main orometer (Life Technologies, Thermo Fisher Scientific Inc., USA) islands and islets of this archipelago, being most of them single- using the dsDNA BR (double stranded DNA Broad Range) Assay island endemics [9]. Although multiple colonization events may Kit and checked its purity in a NanoDrop 2000 spectrophotome- account for the initial origin of species diversity the bulk of this ter (Thermo Fisher Scientific Inc.). diversity is the result of in situ diversification [8]. Dysdera spi- ders have adapted to a broad range of terrestrial habitats within DNA sequencing the Canary Islands [9]. Interestingly, many co-occurring species significantly differ in mouthpart sizes and shapes, presumably We sequenced the genome of D. silvatica using 4 different se- owing to adaptations to a specialized diet [6,7], suggesting that quencing platforms (Table 1; Supplementary Table S1-2). First, stenophagy has evolved multiple times independently in these we used the Illumina HiSeq2000 to obtain the genome sequence islands [10]. Although behavioral and physiological experiments of a single male (100 bp, paired-end [PE] reads, 100 PE; TruSeq have revealed a close correlation between morphological traits library). The flow-cell lane generated ∼51 Gb of sequence, rep- and prey preference in Dysdera, little is known about the molec- resenting a genome coverage of 30× (assuming a genome size ular basis of trophic adaptations in this genus. of ∼1.7 Gb; see below). The genome of a female was sequenced Here we present the draft assembly and functional annota- using a mate pair (MP) approach; for that we used Nextera 5 kb- tion of the genome of the Canary Island endemic spider Dys- insert 100 PE libraries and the HiSeq2000 to generate ∼40 Gb of dera silvatica Schmidt, 1981 (NCBI:txid477319; Fig. 1). This study sequence (∼23× of coverage). A third individual (male) was used is the first genomic initiative within its family and just the sec- for single-molecule real-time (SMRT) sequencing (PacBio long ond within the Synspermiata [11], a clade that includes most reads). We used 8 SMRT libraries (20 kb SMRT bell templates), of the families formerly included in Haplogynae, which was re- which were sequenced using the P6-C4 chemistry in a PacBio cently shown to be paraphyletic [12,13] (Fig. 2). Remarkably, a RSII platform. We obtained a yield of ∼9.6 Gb (raw coverage of Sánchez-Herrero et al. 3 Dysdera silvatica SYNSPERMIATA ARANEOMORPHAE Loxosceles reclusa ARANEAE Parasteatoda tepidariorum Latrodectus hesperus ENTELEGYNAE ARACHNIDA Stegodyphus mimosarum CHELICERATA MYGALOMORPHAE ARTHROPODA Acanthoscurria geniculata Mesobuthus martensii SCORPIONES ECDYSOZOA Centruroides sculpturatus Ixodes scapularis PARASITIFORMES Metaseiulus occidentalis ACARI Tetranychus urticae ACARIFORMES Euroglyphus maynei Limulus polyphemus XIPHOSURA Drosophila melanogaster Bombyx mori HEXAPODA Pediculus humanus Daphnia pulex CRUSTACEA Strigamia maritima MYRIAPODA Hypsibius dujardini TARDIGRADA Caenorhabditis elegans NEMATODA Neo-Proterozoic Paleozoic Mesozoic Cz Figure 2 Phylogenetic relationships of the species used for the D. silvatica genome annotation (see Supplementary Table S1-11 for further details) and completeness analysis. Because the chelicerata phylogeny is controversial (e.g., [20], [21]), we set the most conflictive clades as polytomies. Divergence times were obtained from Carlson et al. (2017) [22] and the TimeTree web server (http://www.timetree.org/). Cz, cretaceous period. Table 1. Sequencing data and library information Coverage Run ID Library Insert size Read lengths Lanes Total bases Raw read pairs (×)a PE Illumina HiSeq200 - Truseq 370 bp 100×100 PE 1 51,202,445,102 506,954,902 30 MP Illumina HiSeq200 - Nextera 5 kb 100×100 PE 1 39,609,522,995 392,173,495 23 Nanopore Nanopore 1D Libraries - Nanopore 5 23,193,357,481 20,534,058 14 PacBio PacBio RSII 20 Kb SMRTbell - SMRT 8 9,652,844,880 1,455,288 6 a Based on the genome size estimated by flow cytometry ∼1.7 Gb. ∼6×). Finally, 2 additional females were used for the 5 runs of ative G0/G1 peak positions of the stained nuclei using a Gallios Nanopore sequencing (Nanopore 1D libraries). We got a yield of flow cytometer (Beckman Coulter, Inc, Fullerton, CA); the results ∼23.2 Gb (∼14× coverage) (Table 1; Supplementary Table S1-2). were based on the average of 3 spider replicates, counting a min- imum of 5,000 cells per individual. In addition, we also estimated the D. silvatica genome size D. silvatica chromosome and genome size from the distribution of k-mers (from short reads) with Jelly- D. silvatica has a diploid chromosome set of 6 pairs of autosomes fish v.2.2.3 (Jellyfish, RRID:SCR 005491) [28]. The distribution of k- and 2 (females are XX; 2n = 14) or 1 (males are X0) sex chromo- mers of size 17, 21, and 41 (GenomeScope (GenomeScope, RRID: somes (M. A. Arnedo, unpublished results). Using flow cytometry SCR 017014) [29]) resulted in a haploid genome size of ∼1.23 Gb and the genome of the German cockroach Blattella germanica (1C (Supplementary Fig. S1). The discrepancy between k-mer– and = 2.025 Gb, J. S. Johnston, personal communication; see also [23]) cytometry-based estimates may be caused by the presence of as reference, we determined that the haploid genome size of D. repetitive elements [30], which can affect k-mer estimates. silvatica is ∼1.7 Gb. For the analysis, we adapted the Hare and Johnston [24] protocol for spiders species, without using male Read preprocessing palps and chelicers to avoid analyzing haploid or endoreplicated cells, respectively [25,26]. Shortly, we isolated cells from the head To avoid including contaminants in the assembly step, we of the male cockroach, and legs and palps from female spiders. searched the raw reads for mitochondrial, bacterial, archaeal, We incubated the cells in LB0.1 with 2% of tween [27], propidium and virus sequences. We downloaded all genomes of all these iodide (50 μg/mL), and RNAse (40 μg/mL). After 10 minutes, the kinds available in the GenBank database (Supplementary Ta- processed tissue was filtered using a nylon mesh of 20 μm. We ble S1-3) and used BLASTN v2.4.0 (BLASTN, RRID:SCR 001598) determined the DNA content of the diploid cells through the rel- [31] to detect and filter all contaminant reads (E-value <10−5 ; 4 The draft genome sequence of the spider Dysdera silvatica (Araneae, Dysderidae) >90% alignment length; >90% identity). We preprocessed raw Table 2. Dysdera silvatica nuclear genome assembly and annotation reads using PRINSEQ v.0.20.3 (PRINSEQ, RRID:SCR 005454) [32]. statistics We estimated some descriptive statistics, such as read length Genome assemblya Value and k-mer representation, and calculated the amount of adapter sequences and exact duplicates. Assembly size (bp) 1,359,336,805 Quality-based trimming and filtering was performed accord- % AT/CG/N 64.91%/34.83%/0.26% ing to the chemistry, technology, and library used (Supplemen- Number of scaffolds 65,205 tary Table S1-4). For the short-insert 100 PE library, we used Trim- Longest scaffold 340,047 momatic v0.36 (Trimmomatic, RRID:SCR 011848) [33] with spe- N50 38,017 cific lists of adapters of the TruSeq v3 libraries to filter all reads L50 10,436 shorter than 36 bp or with minimum quality scores < 30 along Repeat statisticsb 4-bp sliding windows. We also filtered trailing and leading bases Number of elements 3,284,969 with a quality score < 10. Long-insert MP libraries were prepro- Length (bp) [% Genome] 731,540,381 [53.81%] cessed using NxTrim v0.4.1 [34] with default parameters (Sup- Genome annotationa plementary Table S1-4a and b). We preprocessed the raw PacBio Protein-coding genes 48,619 reads using the SMRT Analysis Software (SMRT Analysis Soft- Functionally annotated 36,398 (74.86%) ware, RRID:SCR 002942) [35], by generating circularized consen- Without functional 12,221 (25.14%) sus sequence to further perform a polishing analysis with Pilon annotation v1.22 (Pilon, RRID:SCR 014731) [36] based on short reads (Supple- tRNA genes 33,934 mentary Table S1-4c). a See also Supplementary S1-7. b Summary of the RepeatMasker analysis (See also Supplementary Table S1-9). De novo genome assembly of repetitive sequences, (i) D. silvatica–specific repetitive el- We used MaSuRCA v3.2.9 (MaSuRCA, RRID:SCR 010691) [37] for a ements generated with RepeatModeler v1.0.11 [46], (ii) the hybrid de novo assembly of the D. silvatica genome (Supplemen- Dfam Consensus [48] (version 20170127), and (ii) the RepBase tary Fig. S2). Additionally, we performed a scaffolding phase us- (version 20170127) [49,50]. We identified 2,604 families of repeti- ing AGOUTI (minimum number of joining reads pairs support, k tive elements, where 1,629 of them (62.6%) were completely un- = 3) [38], and the raw reads from a D. silvatica RNA sequencing known. Repetitive sequences accounted for ∼732 Mb, which rep- (RNAseq) experiment [39] (Supplementary Table S1-5 and S1-6). resent 53.8% of the total assembly size (Table 2; Supplementary During the assembly phase, we chose for each software the pa- Table S1-9a). Remarkably, most abundant repeats are from un- rameter values that generated the best assembly (Supplemen- known families, 22.6% of the assembled genome. The repeti- tary Table S1-7) in terms of (i) continuity and contig size statis- tive fraction of the genome also include DNA elements (16.8%), tics, such as the N50, L50, and the total number of sequences and LINEs (10.7%), and SINEs (1.85%), and a small fraction of other el- bases assembled; and (ii) completeness measures, obtained as ements, including LTR elements, satellites, simple repeats, and the fraction (and length) of a series of highly conserved proteins low-complexity sequences. We found that the 10 most abundant present in the draft genome. Particularly, we used 5 datasets, repeat families among the 2,604 identified in D. silvatica account BUSCO v3 (BUSCO, RRID:SCR 015008) with genome option [40] for ∼7% of the genome and encode 5 unknown, 3 SINEs, and 2 using (i) the Arthropoda or (ii) the Metazoa dataset, (iii) the 457 LINEs, with an average length of ∼193, ∼161, and ∼1,040 bp, re- core eukaryotic genes (CEGs) of Drosophila melanogaster [41], (iv) spectively (Supplementary Table S1-9b). the 58,966 transcripts in the D. silvatica transcriptome [39], and We also studied the distribution of the high-covered genome (v) the 9,473 1:1 orthologs across 5 Dysdera species, D. silvatica; regions to describe the spacing pattern among repetitive se- D. gomerensis Strand, 1911; D. verneaui Simon, 1883; D. tilosensis quences. In particular, we searched for genomic regions that Wunderlich, 1992; and D. bandamae Schmidt, 1973 obtained from have a higher than average sequencing coverage above a par- the comparative transcriptomics analysis of these species [42]. ticular threshold. Because repetitive regions are more prone to Finally, we performed an additional search to identify and re- form chimeric contigs in the assembly step, we only used Ma- move possible contaminants in the generated scaffolds (Supple- SuRCA super reads, and longer than 10 kb and free of Ns (34,937 mentary Table S1-7). We discarded 16 contaminant sequences > contigs; 1.12 Gb). We estimated the coverage after mapping the 5 kb. The final assembly size of the D. silvatica genome (Dsil v1.2) short reads (from the 100PE library) to those contigs. We defined was ∼1.36 Gb, with an N50 of ∼38 kb (Table 2). as high-coverage regions (HCRs) those with a coverage ≥2.5× or We determined the average genome coverage for each se- 5× the genome-wide average (∼30×), in a region of ≥150, ≥500, quencing library with SAMtools v1.3.1 (SAMtools, RRID:SCR 002 ≥1,000, or ≥5,000 bp (Supplementary Fig. S4a; Supplementary 105) [43], by mapping short reads (using bowtie2 v2.2.9 [bowtie2, Table S2). We found a large number of contigs encompassing RRID:SCR 005476] [44]) or long reads (using minimap2 [45]) to the ≥1 HCR. For instance, 21,614 contigs (∼61.9%) include ≥1 HCR of final draft assembly (Table 1; Supplementary Table S1-8; Supple- 150 bp with >2.5× coverage (an average of 2.48 HCRs per contig; mentary Fig. S3). 77.7 HCR per Mb) (Supplementary Table S2-2a). For HCRs of >5× coverage, the results are also remarkable (10,604 contigs have ≥1 HCR of 150 bp, corresponding to 25.6 HCR per Mb). As expected, Repetitive DNA sequences the longer the HCR the smaller the fraction in the genome; in- We analyzed the distribution of repetitive sequences in the deed, we found that the genome is encompassing ∼5 HCR per genome of D. silvatica, using either a de novo with RepeatModeler Mb (HCR, longer than 1 kb at 2.5×). The distances between con- v1.0.11 (RepeatModeler, RRID:SCR 015027) [46], or a database- secutive HCRs doenot show clear differences between the 2.5× guided search strategy with RepeatMasker v.4.0.7 (Repeat- and 5× thresholds (Supplementary Fig. 4b and S5; Supplemen- Masker, RRID:SCR 012954) [47]. We used 3 different databases tary Table S2-2b). Sánchez-Herrero et al. 5 Table 3. Completeness analysisa Number Identified BLAST analysisb (%) Parasteatoda genes (n = 30,041) 19,580 (65.2) Single-copy Dysdera (n = 9,473) 8,420 (88.9) Single-copy spiders (n = 2,198) 2,141 (97.4) CEG (n = 457) 438 (95.8) BUSCO analysisc Metazoa (n = 978) Identified BUSCO 882 (90.2) Complete (C) 689 (70.5) Single copy (S) 662 (67.7) Duplicated (D) 27 (2.8) Fragmented (F) 193 (19.7) Missing (M) 96 (9.8) Figure 3 Bar plot of the annotation of the repetitive elements within the HCRs (2.5× threshold) at different intra-HCR length cutoffs (150, 500, 1,000, and Artrhopoda (n = 1,066) 5,000 bp) (Supplementary Table S2-2a). Colors represent the type of repeat ele- Identified BUSCO 959 (89.9) ment identified by RepeatMasker. ”Other types” class includes the LTR elements, Complete (C) 736 (69.1) small RNA, and satellite information that represent a small fraction. Single copy (S) 702 (65.9) Duplicated (D) 34 (3.2) We found a strong relationship between the length of the Fragmented (F) 223 (20.9) Missing (M) 107 (10.0) HCR and the type of the included repetitive elements (Fig. 3; Supplementary Table S2-3). For instance, while LINEs represent a Completeness analysis of the 36,398 functional annotated proteins of D. silvat- 8.62% of the repetitive elements in the whole genome, they are ica. clearly enriched in the HCRs (36.12% in HCRs longer than 150 bp; b − BLASTP searches against different datasets. E-value cutoff < 10 3 , alignment 12.08% in HCRs longer than 5,000 bp) (Fig. 3; Supplementary Ta- length cutoff > 30%, and identity cutoff > 30%. c ble S2-3a); the same was found for the small RNA fraction (ribo- BUSCO analysis using default parameters against different datasets (BUSCO, RRID:SCR 015008). somal RNA). In contrast, the fraction of low-complexity repet- itive sequences is much less represented in small HCRs than in the whole genome (∼1.3%). We also found that the cover- tified in 5 Dysdera species and the 1:1 orthologs among spiders age threshold has little effect on the results (Supplementary Ta- available at OrthoDB v10 (OrthoDB, RRID:SCR 011980) [56] (8,792). ble S2-3; Supplementary Fig. S6), either for the main families or After several iterative training rounds, we applied MAKER2, Au- across subfamilies (Supplementary Table S2-4 and S2-5). gustus, and SNAP, adding other sources of evidence: (i) transcript Given that the HCR analysis covers an important fraction of evidence (Dsil-RefGuided transcriptome), (ii) RNAseq reads exon the assembled bases (∼82%), the present results can likely be junctions generated with HISAT2 [51] and regtools [57], and (iii) extrapolated to the whole genome. Therefore, the relatively low proteins annotated in other arthropods, especially chelicerates N50 of the D. silvatica genome draft is very likely to be caused (Fig. 2; Supplementary Table S1-11). The annotation process re- by abundant interspersed repeats preventing genome continu- sulted in 48,619 protein-coding and 33,934 transfer RNA (tRNA) ity. Despite the low N50 we estimated that the draft presented genes. The mean annotation edit distance (AED) upon protein- here is mostly complete in terms of functional regions (see be- coding genes was 0.32 (Supplementary Fig. S6), which is typical low). of a well-annotated genome [58, 59]. After each training and it- erative annotation round, we checked the improvement of the annotation by means of the cumulative fraction of AED (Supple- Transcriptome assembly and genome mentary Table S1-12a; Supplementary Fig. S7). annotation We searched for the presence of protein domain signatures We used the newly generated genome sequence to obtain a in annotated protein-coding genes using InterProScan v5.15-54 reference-guided assembly of the D. silvatica transcriptome with (InterProScan, RRID:SCR 005829) [60,61], which includes infor- the RNAseq data from Vizueta et al. [39]. We used HISAT2 v2.1.0 mation from public databases (see additional details in Sup- (HISAT2, RRID:SCR 015530) [51] to map the RNAseq reads to plementary Table S1-7). Additionally, we used NCBI BLASTP the reference and Trinity v2.4.0. (Trinity, RRID:SCR 013048) [52] v.2.4.0 (BLASTP, RRID:SCR 001010) [31] (E-value cutoff <10−5 ; (genome-guided bam, max intron = 50 kb, min coverage = 3) >75% alignment length) against the Swiss-Prot database to an- to assemble the transcriptome (named ”Dsil-RefGuided tran- notate D. silvatica genes. We found that 74.9% (36,398 genes) scriptome”; Supplementary Table S1-10). We used the MAKER2 of the predicted protein-coding genes have hits with records v2.31.9 (MAKER2, RRID:SCR 005309) [53] genome annotation of either InterPro (32,322 genes) (InterPro, RRID:SCR 006695) or pipeline for the structural annotation of D. silvatica genes (Sup- Swiss-Prot (17,225 cases) (Table 2; Supplementary Table S1-7). plementary Fig. S2), using both ab initio gene predictions and an- notation evidences from D. silvatica and other sources. For the ab Completeness initio gene predictions we initially trained Augustus v3.1.0 (Au- gustus, RRID:SCR 008417) [54] and SNAP (SNAP, RRID:SCR 00212 We determined the completeness of the D. silvatica genome 7) [55] softwares using scaffolds longer than 20 kb, and BUSCO assembly (Table 3) using BLASTP (E-value cutoff <10−3 ; >30% gene models generated from completeness searches. Then we of alignment length and identity > 50%). We searched for ho- iteratively included a reliable set of proteins for a further train- mologs of the functionally annotated peptides (36,398) (i) among ing. This dataset was composed of the 9,473 orthologs 1:1 iden- CEG genes of Drosophila melanogaster [41]; (ii) among the pre- 6 The draft genome sequence of the spider Dysdera silvatica (Araneae, Dysderidae) dicted peptides of Parasteatoda tepidariorum, a spider with a well- precluded obtaining a high-continuity draft. The characteristic annotated genome [62]; (iii) among the 9,473 1:1 orthologs across holocentric chromosomes of Dysderidae [17] may also explain 5 Dysdera species; and (iv) among the 2,198 single-copy genes the observed genome fragmentation; indeed, it has been re- identified in all spiders and available in OrthoDB v10 [56]. We cently shown that genome-wide centromere-specific repeat ar- found in D. silvatica a high fraction of putative homologs (95.8% rays are interspersed among euchromatin in holocentric plants of CEG genes, and 97.4% spider-specific single-copy genes; Ta- (Rhynchospora, Cyperceae) [66]. ble 3). Furthermore, the analysis based on the putative homologs Nevertheless, the completeness and the extensive annota- of the single-copy genes included in the BUSCO dataset (BUSCO, tions achieved for this genome, as well as the new reference- RRID:SCR 015008) [40], applying the default parameters for the guided transcriptome, make this draft an excellent source tool genome and protein mode, also demonstrated the high com- for further functional and evolutionary analyses in this and pleteness of the genome draft. Indeed the analysis recovered the other related species, including the origin and evolution of rel- ∼90% of Metazoa or Arthropoda genes (v9), and nearly 70% of evant spider traits, such as venom and silk. Moreover, the avail- them are complete in D. silvatica. ability of new genomic information in a lineage with remarkable We extended the search for D. silvatica homologs to a broader evolutionary features such as recurrent colonizations of the un- taxonomic range (Fig. 2; Supplementary Table S1-11) by includ- derground environment or complex reproductive anatomies in- ing other metazoan lineages and performing a series of local dicative of cryptic female choice, to cite 2 examples, will further − BLASTP searches (E-value cutoff < 10 3 ; >30% alignment length). provide valuable knowledge about relevant scientific questions, We found that a great majority of D. silvatica genes are shared such as the molecular basis of adaptation to extreme habitats or among arthropods (57.9%), 11,995 of them (32.95%) also being the genetic drivers of sexual selection, along with more general present in Ecdysozoa (Fig. 4a). Remarkably, 9,560 genes appears aspects related to gene content across main arthropod groups, to be spider-specific, 4,077 of them being specific (unique) of D. the consequences of whole-genome duplications, or phyloge- silvatica. Despite almost all these species-specific genes having netic relationships with the Araneae. Additionally, because this interproscan signatures, the annotation metrics are poor com- genus experienced a spectacular adaptive radiation in the Ca- pared with genes having homologs in other species (Supplemen- nary Islands, the present genome draft could be useful to further tary Table S1-12b; Supplementary Figs S7 and S9); indeed, they studies investigating the genomic basis of island radiations. have an average number of exons (2.8) and gene length (∼168aa), which may reflect their partial nature. They could be part of very large genes interspersed by repeats or complex sequences diffi- Availability of supporting data and materials cult to assemble. The analysis using OrthoDB (v10) [56] across The whole-genome shotgun project has been deposited at 5 chelicerates (including D. silvatica) identified 1,798 genes, with DDBJ/ENA/GenBank under accession number QLNU00000000 1:1 orthologous relationships (Fig. 4b), while 12,101 D. silvatica and project ID PRJNA475203. The version described in this article genes showed other more complex orthologous/homologous re- is version QLNU01000000. This project repository includes raw lationships (Fig. 4b, Supplementary Table S1-12c and S3-1). The data, sequencing libraries information, and assemblies of the analysis across the genome annotations of some representative mitochondrial and nuclear genomes. Other relevant datasets arthropods identified 950 genes with 1:1 orthologous relation- such as annotation, reference-guide assembled transcripts, re- ships (Supplementary Fig. S8, Supplementary Table S1-12c and peat, and HCR data, as well as other data relevant for the repro- S3-2). ducibility of results, are available in the GigaDB dataset [67]. Mitochondrial genome assembly and Additional file annotation File S1. Supplemental Material Summary We assembled the mitochondrial genome of D. silvatica (mtDsil) SanchezHerrero Dsilvatica SupMaterial Summary.pdf from 126,758 reads identified in the 100PE library by the soft- ware NOVOPlasty [63]. Our de novo assembly yielded a unique contig of 14,440 bp (coverage of 878×) (Supplementary Table S1- Availability of supporting source code and 13). CGVIEW (CGVIEW, RRID:SCR 011779) [64] was used to gen- requirements erate a genome visualization of the annotated mtDsil genome The scripts employed and developed in this project are available (Supplementary Fig. S10). We identified 2 ribosomal RNAs, 13 under the github repository: protein-coding genes, and 15 tRNAs (out of the putative 22 tR- Project name: Genome assembly of Dysdera silvatica NAs). Based on the contig length and the inability of standard Project home page: https://github.com/molevol-ub/Dysdera sil automatic annotation algorithms to identify tRNA with missing vatica genome arms, as reported for spiders [65], the complete set of tRNAs is Operating system(s): Platform independent most likely present for this species. Programming language: Bash, Perl, Python, R License: MIT Conclusion Abbreviations We have reported the assembly and annotation of the nuclear and mitochondrial genomes of the first representative of the AED: annotation edit distance; AGOUTI: Annotated Genome Op- spider superfamily Dysderoidea and the second genome of a timization Using Transcriptome Information; BLAST: Basic Local Synspermiata, one of the main evolutionary lineages within the Alignment Tool; bp: base pair; BUSCO: Benchmarking Univer- “true spiders” (Araneomorphae) and still sparsely sampled at sal Single Copy Orthologs; CEG: core eukaryotic gene; Cz: Cre- the genomic level [14]. Despite the high coverage and the hybrid taceous period; Dsil: Dysdera silvatica; Gb: gigabase pairs; GC: assembly strategy, the repetitive nature of the D. silvatica genome guanine cytosine; GO: Gene Ontology; HCR: high-coverage re- Sánchez-Herrero et al. 7 (a) (b) Figure 4 (a) Pie chart illustrating the taxonomic distribution of positive BLAST hits of the D. silvatica protein-coding genes against the sequence data of species included in Fig. 2. (b) Homology relationships among D. silvatica (Dsil) and different chelicerates genomes available in OrthoDB v10 [56], Parasteatoda tepidariorum (Ptep), Stegody- phus mimosarum (Smim), Ixodes scapularis (Isca), and Tetranychus urticae (Turt). Red and orange bars indicate the fraction of single-copy genes (1:1 orthologs) identified in all species, and in all but 1 (e.g., missing in 1 species), respectively. The dark and light green bar indicate the fraction of orthologs present in all species and in all but 1, respectively, that are not included in previous categories. The blue bar (other orthology/homology) shows other more complex homologous relationships. The results were generated by uploading D. silvatica proteins to the OrthoDB web server. gions; Isca: Ixodes scapularis; kb: kilobase pairs; LINE: long in- Acknowledgments terspersed nuclear element; LTR: long terminal repeats; Ma- We acknowledge the Garajonay National Parks for granting col- SuRCA: Maryland Super-Read Celera Assembler; Mb: megabase lection permits and helping with lodging and logistics during pairs; MP: mate pair; Mya: million years ago; NCBI: National Cen- fieldwork. We also thank CNAG (Centro Nacional de Análisis ter for Biotechnology Information; PacBio: Pacific Biosciences; Genómico) for the Nanopore sequencing facilities. PE: paired-end; PRINSEQ: PReprocessing and INformation of SE- Quence data; Ptep: Parasteatoda tepidariorum; RNAseq: RNA se- quencing; SINE: short interspersed nuclear element; Smim: Ste- References godyphus mimosarum; SMRT: Single-Molecule Real Time; tRNA: 1. World Spider Catalog (2018). 2018. http://wsc.nmbe.ch. Ac- transfer RNA; Turt: Tetranychus urticae. cessed on April 2019. 2. Pekár S, Toft S. Trophic specialisation in a predatory group: Competing interests the case of prey-specialised spiders (Araneae). Biol Rev 2015;90(3):744–61. The authors declare that they have no competing interests. 3. Hopkin SP, Martin MH. Assimilation of zinc, cadmium, lead, copper, and iron by the spider Dysdera crocata, a predator of Funding woodlice. Bull Environ Contam Toxicol 1985;34:183–87. 4. Pekár S, Lı́znarová E, Řezác M. Suitability of woodlice prey This study was supported by the Ministerio de Economı́a y for generalist and specialist spider predators: a comparative Competitividad of Spain (CGL2012-36863, CGL2013-45211, and study. Ecol Entomol 2016;41(2):123–30. CGL2016-75255), and by the Comissió Interdepartamental de Re- 5. Toft S, Macı́as-Hernández N. Metabolic adaptations for iso- cerca I Innovació Tecnològica of Catalonia, Spain (2014SGR-1055 pod specialization in three species of Dysdera spiders from and 2014SGR1604). J.F.S.-H. was supported by a Formación del the Canary Islands. Physiol Entomol 2017;42(2):191–98. Profesor Universitario (FPU) grant (Ministerio de Educación of 6. Řezáč M, Pekár S. Evidence for woodlice-specialization Spain, FPU13/0206); C.F.-L. by an IRBio PhD grant; S.H-A by Be- in Dysdera spiders: behavioural versus developmental ap- cas Postdoctorales en el Extranjero CONACyT; A.S.-G. by a Beat- proaches. Physiol Entomol 2007;32(4):367–71. riu de Pinós grant (Generalitat de Catalunya, 2010-BP-B 00175); 7. Řezác M, Pekár S, Lubin Y. How oniscophagous spiders over- and J.R. and M.A.A. were partially supported by ICREA Academia come woodlouse armour. J Zool 2008;275(1):64–71. (Generalitat de Catalunya). 8. Arnedo MA, Oromı́ P, Ribera C. Radiation of the spider genus Dysdera (Araneae, Dysderidae) in the Canary Islands: cladistic assessment based on multiple data sets. Cladistics Authors’ contributions 2001;17:313–353. J.R., A.S.-G., and M.A.A designed the study. C.F.-L., J.F.S.-H., 9. Macias-Hernandez N, de la Cruz López S, Roca-Cusachs M, P.E., and S.H-A. processed the samples and extracted DNA. et al. A geographical distribution database of the genus Dys- J.F.S.-H. performed the bioinformatics analysis and drafted the dera in the Canary Islands (Araneae, Dysderidae). Zookeys manuscript. J.F.S.-H., A.S.-G., and J.R. interpreted the data. All au- 2016;625(625):11–23. thors revised and approved the final manuscript. 10. Arnedo MA, Oromı́ P, Múrria C, et al. The dark side of an is- 8 The draft genome sequence of the spider Dysdera silvatica (Araneae, Dysderidae) land radiation: systematics and evolution of troglobitic spi- formatics 2017;33(14):2202–4. ders of the genus Dysdera Latreille (Araneae: Dysderidae) in 30. Austin CM, Tan MH, Harrisson KA, et al. De novo genome the Canary Islands. Invertebr Syst 2007;21(6):623. assembly and annotation of Australia’s largest freshwater 11. Michalik P, Ramı́rez MJ. Evolutionary morphology of the male fish, the Murray cod (Maccullochella peelii), from Illumina and reproductive system, spermatozoa and seminal fluid of spi- Nanopore sequencing read. GigaScience 2017;6(8):1–6. ders (Araneae, Arachnida) - current knowledge and future 31. Altschul SF, Gish W, Miller W, et al. Basic Local Alignment directions. Arthropod Struct Dev 2014;43(4):291–322. Search Tool. J Mol Biol 1990;215(3):403–10. 12. Wheeler WC, Coddington JA, Crowley LM, et al. The spi- 32. Schmieder R, Edwards R. Fast identification and removal of der tree of life: phylogeny of Araneae based on target- sequence contamination from genomic and metagenomic gene analyses from an extensive taxon sampling. Cladistics datasets. PLoS One 2011;6(3):e17288. 2017;33(6):574–616. 33. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexi- 13. Fernández R, Kallal RJ, Dimitrov D, et al. Phylogenomics, ble trimmer for Illumina sequence data. Bioinformatics diversification dynamics, and comparative transcriptomics 2014;30(15):2114–20. across the spider tree of life. Curr Biol 2018;28(9):1489–97. 34. O’Connell J, Schulz-Trieglaff O, Carlson E, et al. NxTrim: opti- 14. Garb JE, Sharma PP, Ayoub NA. Recent progress and mized trimming of Illumina mate pair reads. Bioinformatics prospects for advancing arachnid genomics. Curr Opin In- 2015;31(12):2035–7. sect Sci 2018;25:51–7. 35. PacBio. Single Molecule Real Time (SMRT). https: 15. Vizueta J, Rozas J, Sánchez-Gracia A. Comparative genomics //www.pacb.com/products-and-services/analytical-soft reveals thousands of novel chemosensory genes and mas- ware/smrt-analysis/. sive changes in chemoreceptor repertories across chelicer- 36. Walker BJ, Abeel T, Shea T, et al. Pilon: an integrated tool ates. Genome Biol Evol 2018;10(5):1221–36. for comprehensive microbial variant detection and genome 16. Deeleman-Reinhold CL. The genus Rhode and the assembly improvement. PLoS One 2014;9(11):e112963. harpacteine genera Stalagtia, Folkia, Minotauria, and Kaemis 37. Zimin AV, Marçais G, Puiu D, et al. The MaSuRCA genome (Araneae, Dysderidae) of Yugoslavia and Crete, with remarks assembler. Bioinformatics 2013;29(21):2669–77. on the genus Harpactea. Rev Arachnol 1993;10(6):105–35. 38. Zhang SV, Zhuo L, Hahn MW. AGOUTI: improving genome 17. Diaz MO, Maynard R, Brum-Zorrilla N. Diffuse centromere assembly and annotation using transcriptome data. Giga- and chromosome polymorphism in haplogyne spiders of the Science 2016;5(1):31. families dysderidae and segestriidae. Cytogenet Genome Res 39. Vizueta J, Frias-López C, Macı́as-Hernández N, et al. Evolu- 2010;128(1-3):131–8. tion of chemosensory gene families in arthropods: insight 18. Uhl G. Two distinctly different sperm storage organs in from the first inclusive comparative transcriptome analysis female Dysdera erythrina (Araneae: Dysderidae). Arthropod across spider appendages. Genome Biol Evol 2017;9(1):178– Struct Dev 2000;29(2):163–9. 96. 19. Burger M, Kropf C. Genital morphology of the haplogyne 40. Simão FA, Waterhouse RM, Ioannidis P, et al. BUSCO: assess- spider Harpactea lepida (Arachnida, Araneae, Dysderidae). ing genome assembly and annotation completeness with Zoomorphology 2007;126(1):45–52. single-copy orthologs. Bioinformatics 2015;31(19):3210–2. 20. Ballesteros JA, Sharma PP. A critical appraisal of the 41. Parra G, Bradnam K, Ning Z, et al. Assessing the gene space placement of Xiphosura (chelicerata) with account of in draft genomes. Nucleic Acids Res 2009;37(1):289–97. known sources of phylogenetic error. Syst Biol 2019, 42. Vizueta, J., Macı́as-Hernández, N., Arnedo, MA., Rozas, doi:10.1093/sysbio/syz011. J.and Sánchez-Gracia, A.(2019) Chance and predictabil- 21. Lozano-Fernandez J, Tanner AR, Giacomelli M, et al. Increas- ity in evolution: the genomic basis of convergent di- ing species sampling in chelicerate genomic-scale datasets etary specializations in an adaptive radiation. Mol. Ecol. provides support for monophyly of Acari and Arachnida. Nat doi:10.1111/mec.1519931359512 Commun 2019;10:2295. 43. Li H, Handsaker B, Wysoker A, et al. The Sequence 22. Carlson DE, Hedin M. Comparative transcriptomics of Alignment/Map format and SAMtools. Bioinformatics Entelegyne spiders (Araneae, Entelegynae), with empha- 2009;25(16):2078–9. sis on molecular evolution of orphan genes. PLoS One 44. Langmead B, Salzberg SL. Fast gapped-read alignment with 2017;12(4):e0174102. Bowtie 2. Nat Methods 2012;9(4):357–9. 23. Gregory TR. Animal Genome Size Database. 2018. http://ww 45. Li H. Minimap2: pairwise alignment for nucleotide se- w.genomesize.com. quences. Bioinformatics 2018;34(18):3094–100. 24. Hare EE, Johnston JS. Genome size determination using flow 46. Smit AF, Hubley R. RepeatModeler Open-1.0. 2008–2015. http: cytometry of propidium iodide-stained nuclei. Methods Mol //www.repeatmasker.org. Biol 2011;772:3–12. 47. Smit AF, Hubley R, Green P. RepeatMasker Open-3.0. 1996– 25. Rasch EM, Connelly BA. Genome size and endonuclear DNA 2010. http://www.repeatmasker.org. replication in spiders. J Morphol 2005;265(2):209–14. 48. Wheeler TJ, Clements J, Eddy SR, et al. Dfam: a database of 26. Gregory TR, Shorthouse DP. Genome sizes of spiders. J Hered repetitive DNA based on profile hidden Markov models. Nu- 2003;94(4):285–90. cleic Acids Res 2012;41(D1):D70–D82. 27. Dpooležel J, Binarová P, Lcretti S. Analysis of nuclear 49. Bao W, Kojima KK, Kohany O. Repbase Update, a database DNA content in plant cells by flow cytometry. Biol Plant of repetitive elements in eukaryotic genomes. Mobile DNA 1989;31(2):113–20. 2015;6(1):11. 28. Marçais G, Kingsford C. A fast, lock-free approach for effi- 50. Jurka J, Kapitonov VV, Pavlicek A, et al. Repbase Update, cient parallel counting of occurrences of k-mers. Bioinfor- a database of eukaryotic repetitive elements. Cytogenet matics 2011;27(6):764–70. Genome Res 2005;110(1-4):462–7. 29. Vurture GW, Sedlazeck FJ, Nattestad M, et al. GenomeScope: 51. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced fast reference-free genome profiling from short reads. Bioin- aligner with low memory requirements. Nat Methods Sánchez-Herrero et al. 9 2015;12(4):357–60. improving coverage, classification and access to protein se- 52. Haas BJ, Papanicolaou A, Yassour M, et al. De novo tran- quence annotations. Nucleic Acids Res 2019;47(D1):D351–60. script sequence reconstruction from RNA-seq using the Trin- 61. Jones P, Binns D, Chang HY, et al. InterProScan 5: ity platform for reference generation and analysis. Nat Pro- genome-scale protein function classification. Bioinformatics toc 2013;8(8):1494–512. 2014;30(9):1236–40. 53. Holt C, Yandell M. MAKER2: an annotation pipeline and 62. Schwager EE, Sharma PP, Clarke T, et al. The house spider genome-database management tool for second-generation genome reveals an ancient whole-genome duplication dur- genome projects. BMC Bioinformatics 2011;12(1):491. ing arachnid evolution. BMC Biol 2017;15(1):62. 54. Stanke M, Steinkamp R, Waack S, et al. AUGUSTUS: a web 63. Dierckxsens N, Mardulyn P, Smits G. NOVOPlasty: de novo server for gene finding in eukaryotes. Nucleic Acids Res assembly of organelle genomes from whole genome data. 2004;32(Web Server issue):W309–12. Nucleic Acids Res 2016;45(4):gkw955. 55. Korf I. Gene finding in novel genomes. BMC Bioinformatics 64. Stothard P, Wishart DS. Circular genome visualization and 2004;5(1):59. exploration using CGView. Bioinformatics 2005;21(4):537–9. 56. Kriventseva EV, Kuznetsov D, Tegenfeldt F, et al. Or- 65. Masta SE, Boore JL. The complete mitochondrial genome thoDB v10: sampling the diversity of animal, plant, fun- sequence of the spider Habronattus oregonensis reveals re- gal, protist, bacterial and viral genomes for evolutionary arranged and extremely truncated tRNAs. Molec Biol Evol and functional annotations of orthologs. Nucleic Acids Res 2004;21(5):893–902. 2019;47(D1):D807–D811. 66. Marques A, Ribeiro T, Neumann P, et al. Holocentromeres in 57. Feng YY, Ramu A, Cotto KC, et al. RegTools: integrated analy- Rhynchospora are associated with genome-wide centromere- sis of genomic and transcriptomic data for discovery of splic- specific repeat arrays interspersed among euchromatin. Proc ing variants in cancer. bioRxiv 2018, doi:10.1101/436634. Natl Acad Sci U S A 2015;112(44):13633–8. 58. Eilbeck K, Moore B, Holt C, et al. Quantitative measures for 67. Sánchez-Herrero JF, Frı́as-López C, Escuer P, et al. Supporting the management and comparison of annotated genomes. data for “The draft genome sequence of the spider Dysdera BMC Bioinformatics 2009;10(1):67. silvatica (Araneae, Dysderidae): a valuable resource for func- 59. Yandell M, Ence D. A beginner’s guide to eukaryotic genome tional and evolutionary genomic studies in chelicerates.” Gi- annotation. Nature Rev Genet 2012;13(5):329–42. gaScience Database 2019; http://dx.doi.org/10.5524/100628. 60. Mitchell AL, , Attwood TK, Babbitt PC, et al. InterPro in 2019: