Mark A Phuong
email: markphuong@gmail.com

Michael E Alfaro

Gusti N Mahardika

Ristiyanti M Marwoto

Romanus Edy Prabowo

Thomas Rintelen

Philipp W H Vogt

Jonathan R Hendricks

Nicolas Puillandre

Thomas Von Rintelen

Philipp Wh Vogt

Lack of signal for the impact of venom gene diversity on speciation Mark A

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Why are some taxa more diverse than others? Species richness and phenotypic diversity are not distributed evenly across the tree of life [START_REF] Rabosky | Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation[END_REF]. For example, there exists over 10,000 species of birds, but their closest relatives (crocodiles and alligators) comprise only of 23 species. Differences in evolvability, or lineages' intrinsic capacity to adapt and diversify, is one reason commonly used to explain these disparities [START_REF] Wagner | Perspective: Complex Adaptations and the Evolution of Evolvability[END_REF][START_REF] Yang | Modularity, evolvability, and adaptive radiations: a comparison of the hemiand holometabolous insects[END_REF][START_REF] Jones | The mutation matrix and the evolution of evolvability[END_REF][START_REF] Pigliucci | Is evolvability evolvable? Nature reviews[END_REF][START_REF] Losos | Adaptive radiation, ecological opportunity, and evolutionary determinism[END_REF]. Evolvability is thought to be determined by the underlying genetic architecture of organismssome genomes of organisms have a greater propensity to generate variation that may be adaptive in the future [START_REF] Wagner | Perspective: Complex Adaptations and the Evolution of Evolvability[END_REF][START_REF] Jones | The mutation matrix and the evolution of evolvability[END_REF][START_REF] Pigliucci | Is evolvability evolvable? Nature reviews[END_REF]. For example, gene duplication increases evolvabilitythe copied gene is free from the selective pressures of the original gene [START_REF] Crow | What is the role of genome duplication in the evolution of complexity and diversity?[END_REF].

Mutation, selection, and drift can act on the copied gene, facilitating the possibility of new phenotypes to arise; this shapes the extent that taxa can diversify and exploit resources [START_REF] Crow | What is the role of genome duplication in the evolution of complexity and diversity?[END_REF]. Over long evolutionary time scales, clades that exhibit higher evolvability are predicted to have increased species richness and diversification rates [START_REF] Yang | Modularity, evolvability, and adaptive radiations: a comparison of the hemiand holometabolous insects[END_REF].

Despite the ubiquity of this concept in macroevolutionary theory, few studies explicitly test these predictions; this is possibly due to the difficulty of identifying genes responsible for phenotype [START_REF] Hoekstra | The Locus of Evolution: Evo Devo and the Genetics of Adaptation[END_REF]. Past studies that have attempted to test the impact of evolvability on diversification have produced mixed results [START_REF] Santini | Did genome duplication drive the origin of teleosts? A comparative study of diversification in ray-finned fishes[END_REF][START_REF] Soltis | Polyploidy and angiosperm diversification[END_REF][START_REF] Mayrose | Recently Formed Polyploid Plants Diversify at Lower Rates[END_REF][START_REF] Rabosky | Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation[END_REF][START_REF] Zhan | Comparative analysis reveals that polyploidy does not decelerate diversification in fish[END_REF][START_REF] Tank | Nested radiations and the pulse of angiosperm diversification : increased diversification rates often follow whole genome duplications[END_REF]Malmstrøm . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF]et al. 2016). For example, whole genome duplication events, which are hypothesized to increase the genomic potential of organisms, have been documented to increase [START_REF] Santini | Did genome duplication drive the origin of teleosts? A comparative study of diversification in ray-finned fishes[END_REF][START_REF] Soltis | Polyploidy and angiosperm diversification[END_REF][START_REF] Tank | Nested radiations and the pulse of angiosperm diversification : increased diversification rates often follow whole genome duplications[END_REF], decrease [START_REF] Mayrose | Recently Formed Polyploid Plants Diversify at Lower Rates[END_REF], and have no impact [START_REF] Zhan | Comparative analysis reveals that polyploidy does not decelerate diversification in fish[END_REF] on the long-term evolutionary success of clades. In another case, a positive correlation between evolvability and speciation rates exist when measuring evolvability through morphological proxies [START_REF] Rabosky | Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation[END_REF]. One limitation of past research on this hypothesis is the inability to tie genomic changes with ecological factors driving diversification patterns [START_REF] Robertson | Lineage-specific rediploidization is a mechanism to explain time-lags between genome duplication and evolutionary diversification[END_REF]. Although gene duplication and whole genome duplication events can increase the evolutionary capacity of organisms, genes that are ecologically relevant for adaptation may not be readily available for selection to drive divergence.

Here, we study the relationship between evolvability and diversification in cone snails (family, Conidae), a diverse group (> 900 spp.) of predatory marine gastropods. These snails feed on either worms, molluscs, or fish by paralyzing their prey with a cocktail of venomous neurotoxins (conotoxins, [START_REF] Duda | Molecular genetics of ecological diversification: duplication and rapid evolution of toxin genes of the venomous gastropod Conus[END_REF]. Cone snail provides a unique opportunity to test predictions of evolvability and diversification for the following reasons: first, cone snail species share an ecologically relevant trait, venom. Conidae species are globally distributed in tropical and subtropical regions, where >30 species can co-occur within the same habitat [START_REF] Kohn | Maximal species richness in Conus: diversity, diet and habitat on reefs of northeast Papua New Guinea[END_REF]. High numbers of species hypothesized to be able to co-occur because species have diversified to specialize on different prey using prey-specific conotoxins [START_REF] Duda | Molecular genetics of ecological diversification: duplication and rapid evolution of toxin genes of the venomous gastropod Conus[END_REF]. Second, venom genes are known and diversify through gene duplication [START_REF] Duda | Evolutionary diversification of multigene families: allelic selection of toxins in predatory cone snails[END_REF][START_REF] Kaas | Conopeptide characterization and classifications: an analysis using ConoServer[END_REF][START_REF] Kaas | ConoServer: updated content, knowledge, and discovery tools in the conopeptide database[END_REF][START_REF] Chang | Extensive and continuous duplication facilitates rapid evolution and diversification of gene families[END_REF]. Diet specialization is thought to be enabled by the rapid evolution of the genes that underlie conotoxinsestimated rates of gene duplication and nonsynonymous substitutions rates for conotoxin genes are the highest across . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] metazoans [START_REF] Duda | Evolutionary diversification of multigene families: allelic selection of toxins in predatory cone snails[END_REF][START_REF] Chang | Extensive and continuous duplication facilitates rapid evolution and diversification of gene families[END_REF]. Therefore, conotoxin genes provide a natural way to characterize differences in evolvability between clades.

We employ a sequence capture technique previously used in cone snails [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] to recover phylogenetic markers and conotoxin genes from 314 described species. We use the phylogenetic markers to reconstruct a time-calibrated phylogeny and perform trait-dependent diversification analyses to test the impact of evolvability on diversification patterns. We predict that clades with a greater number of conotoxin gene copies should have higher speciation rates. In addition, we test other traits that may have an impact on diversification patterns, including diet and larval dispersal strategy.

Methods

Bait design

We used a targeted sequencing approach to recover markers for phylogenetic inference and obtain an estimate of conotoxin gene diversity from Conidae species. For the phylogenetic markers, we identified loci using a previous Conidae targeted sequencing dataset [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]) and the Conidae transcriptome data from [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]). In the Conidae targeted sequencing dataset, the authors generated a phylogeny using 5883 loci across 32 species [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. For our sequencing experiment, we only retained loci that were >180bp and were present in at least 26 out of 32 taxa with at least 10X coverage. We chose to only include longer loci to increase confidence in identifying orthologous fragments in other Conidae species. To identify additional phylogenetic markers from the transcriptome data . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF], which consisted of venom duct transcriptomes from 12 species, we performed the following:

(1) identified reciprocal best blast hits between the assembled transcriptome and the Lottia gigantea protein reference (Simakov et al. 2013) using BLAST+ v2.2.31 (evalue = 1e-10). We also considered fragments that had their best hit to the protein reference, but to a non-overlapping portion (<20% overlapping).

(2) mapped reads using bowtie2 v2.2.7 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] (3) removed duplicates using picard-tools v.2.1.1 (http://broadinstitute.github.io/picard) (4) fixed assembly errors by calling single nucleotide polymorphisms (SNPs) using samtools v1.3 and bcftools v1.3 [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] (5) aligned sequences per locus using mafft v7.222 [START_REF] Katoh | MAFFT version 5: Improvement in accuracy of multiple sequence alignment[END_REF] (6) calculated uncorrected pairwise distances within each locus for all possible pairwise comparisons (7) removed sequences if the uncorrected pairwise distance was greater than the 90 th percentile for those pair of species (8) denoted exon boundaries by comparing the transcriptome sequences to the Lottia gigantea genome reference [START_REF] Simakov | Insights into bilaterian evolution from three spiralian genomes[END_REF], retaining exons >180bp For all retained phylogenetic markers, we also performed the following: (1) we generated an ancestral sequence using FastML v3.1 (Ashkenazy et al. 2012) between a Californiconus californicus sequence and another Conidae sequence that had the highest amount of overlap with the C. californicus sequence (we generated these ancestral sequences to decrease the genetic distances between the target sequence and the orthologous sequence from any Conidae species),

(2) removed sequences that had a GC content < 30% or > 70% because extreme GC contents can . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] reduce capture efficiency (Bi et al. 2012), (3) removed loci that contained repeats identified through the RepeatMasker v4.0.6 web server [START_REF] Smit | RepeatMasker Open-4[END_REF], and (4) performed a self-blast with the target sequences via blastn v2.2.31 (evalue = 1e-10) and removed loci that did not blast to itself with sequence identity >90%. The final set of target loci for phylogenetic inference included 1749 loci, with a total length of 470,435 bp.

To recover conotoxin loci, we targeted sequences generated from both the previous targeted sequencing dataset [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] and the transcriptome dataset [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]. For conotoxin sequences discovered from the targeted sequencing dataset [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF], we performed the following to generate our target sequences: (1) we trimmed each sequence to only retain the coding region and included 100bp flanking the exon,

(2) merged sequences using cd-hit v4.6.4 [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein[END_REF]) at 95% sequence similarity to reduce redundancy among conotoxin loci (3) masked repeats using the RepeatMasker v4.0.6 web server [START_REF] Smit | RepeatMasker Open-4[END_REF], and (4) retained loci >120bp to ensure that the locus was longer than our desired bait sequence length. We concatenated all sequences below 120bp to create a single, chimeric sequence for capture. The final set of target sequences from the previous targeted sequencing dataset consisted of 12,652 unique loci totaling 3,113,904 bp and a single concatenated sequence representing 351 merged loci with a total length of 37,936 bp. We also targeted conotoxin loci from the transcriptomes described in [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]) to obtain conotoxin loci from gene superfamilies that were not targeted in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] or performed poorly. We performed the following to generate a set of conotoxin loci from the transcriptome data: (1) we trimmed sequences from [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF] to only include the coding region and 100bp of the untranslated regions (UTRs), (2) merged sequences using cd-hit v4.6.4 [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein[END_REF]) at 97% sequence similarity to reduce redundancy among conotoxin . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; loci, and (3) masked repeats using the RepeatMasker v4.0.6 web server [START_REF] Smit | RepeatMasker Open-4[END_REF]. This filtered dataset contained 395 conotoxin loci with a total length of 171,317 bp.

We submitted the following datasets to MYcroarray (Ann Arbor, Michigan, USA) for bait synthesis: (1) 1749 loci for phylogenetic inference, (2) 12652 conotoxin loci using data from [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF], (3) a single concatenated sequence using data from [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF], and (4) 395 additional conotoxin loci using transcriptome data from [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]. We chose to synthesize a MYbaits-3 kit, which included 60,000 bait sequences to accommodate all the targeted loci. Because our aim was to recover sequences from species throughout Conidae, each bait sequence was 120bp in length, which increases the efficiency of recovering divergent fragments. We used a 2X tiling density strategy (a new probe every 60bp) across the sequences from datasets (1) and ( 2) and used a 4X tiling density strategy (a new probe every 30bp) across datasets (3) and (4). We chose to increase the tiling density for datasets (3) and (4) because the boundaries between exons were not denoted and we wanted to ensure effective capture of the conotoxin loci. The set of probe sequences will be made available on DRYAD following publication.

Genetic samples, library preparation, hybridization, and sequencing

We performed the targeted sequencing experiment across 362 samples representing both described Conidae species and unique lineages/potential new species identified during routine species verification using the mitochondrial locus (results not shown), CO1 (Table S1, [START_REF] Folmer | DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates[END_REF]. We also sequenced Bathyoma sp. as an outgroup based on a recent molecular phylogeny of the Conoideans, a clade of gastropods that includes Conidae (Table S1, [START_REF] Puillandre | The dragon tamed? A molecular phylogeny of the Conoidea (Gastropoda)[END_REF]. We obtained these genetic samples from two field expeditions in Indonesia and Australia and from five museum collections (Table S1). We extracted DNA from tissue using the EZNA Mollusc DNA kit (Omega Bio-Tek, Doraville, GA, USA). There was slight variation in tissue preservation strategy among samples, with most tissues preserved directly in 95% ethanol (Table S1). For 10 samples, tissue was not available but DNA was available from a previous extraction. For these samples, we ran the DNA through the EZNA Mollusc DNA kit to purify the DNA prior to library preparation. We extracted a minimum of 2000 ng per sample prior to library preparation, when possible. We sheared DNA using a Biorupter when necessary and used a 1X bead purification protocol to ensure that the DNA fragments per sample ranged from 250-600bp, centered on ~350bp. We aimed to generate libraries with longer fragment sizes to ensure that we could recover exons containing the mature toxin region, which are often only recoverable because they are flanking conserved regions that are targeted by our bait design [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF].

We prepared libraries for following the [START_REF] Meyer | Illumina sequencing library preparation for highly multiplexed target capture and sequencing[END_REF] protocol with the following modifications: (1) we started library preparation with at least 2000ng, rather than the 500ng suggested by the protocol to increase downstream capture efficiency, (2) we performed 1X bead clean-up for all enzymatic reactions and (3) generated dual-indexed libraries by incorporating adapters with unique 7bp barcodes. We were able to re-use libraries for the 32 species sequenced in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] and incorporated new indexes for these samples.

We generated equimolar pools of 8 samples and hybridized probes with 2000ng of the pooled DNA for ~24 hours. We substituted the adapter blocking oligonucleotides provided by MYcroarray with custome xGen blocking oligonucleotides (Integrated DNA technologies). We performed 3 independent post-capture amplifications using 12 PCR cycles and pooled these . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] products. We sequenced all samples across 5 Illumina HiSeq 4000 lanes with 100bp paired-end reads. We multiplexed 80 samples per lane for the first 4 lanes and multiplexed the remaining 43 samples on the last lane. Sequencing was carried out at the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley. We note that our third lane containing 80 samples was contaminated, with 65% of the reads belonging to corn DNA. We were able to resequence this entire lane, resulting in overall increased sequencing effort for samples belonging to our third lane (Table 1).

Data filtration and initial assembly

We filtered the raw read data as follows:

(1) we trimmed reads using Trimmomatic v0.36 under the following conditions: (a) we used the ILLUMINACLIP option to trim adapters with a seed mismatch threshold of 2, a palindrome clip threshold of 40, and a simple clip threshold of 15, (b) we performed quality trimming used the SLIDINGWINDOW option with a window size of 4 and a quality threshold of 20, (c) we removed reads below 36bp by setting the MINLEN option to 36, and (d) we removed leading and trailing bases under a quality threshold of 15.

(2) we merged reads using FLASH v1.2.11 [START_REF] Magoč | FLASH: fast length adjustment of short reads to improve genome assemblies[END_REF] with a min overlap parameter of 5, a max overlap parameter of 100, and a mismatch ratio of 0.05.

(3) we removed low complexity reads using prinseq v0.20.4 [START_REF] Schmieder | Quality control and preprocessing of metagenomic datasets[END_REF] using the entropy method with a conservative threshold of 60.

We assembled the filtered read data using SPAdes v3.8.1 using default parameters and reduced redundancy in the resultant assemblies with cap3 [START_REF] Huang | CAP3: A DNA sequence assembly program[END_REF] under default parameters and cd-hit v4.6 [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein[END_REF], sequence identity threshold = 99%).

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Phylogenetic data processing and filtering

To associate assembled contigs with the target sequences for phylogenetic inference, we used blastn v2.2.31 (word size = 11, evalue = 1e-10). For the set of target sequences that originated from the transcriptome dataset, we redefined exon/intron boundaries using EXONERATE v2.2.0 [START_REF] Slater | Automated generation of heuristics for biological sequence comparison[END_REF] using the est2genome model because we found that several predicted exons actually consisted of several smaller exons. For each sample, we mapped reads using bowtie2 (very sensitive local and no discordant options enabled) to a reference that contained only sequences associated with the targeted phylogenetic markers. We marked duplicates using picard-tools v2.0.1 and masked all regions below 4X coverage and removed the entire sequence if more than 30% of the sequence was below 4X coverage. We called SNPs using samtools v1.3 and bcftools v1.3 and estimated average heterozygosity across all contigs within a sample. We removed sequences if a contig had a heterozygosity value greater than two standard deviations away from the mean.

Conotoxin assembly, processing, and filtering

Commonly used assembly programs are known to poorly reconstruct all copies of multilocus gene families [START_REF] Lavergne | Optimized deep-targeted proteotranscriptomic profiling reveals unexplored Conus toxin diversity and novel cysteine frameworks[END_REF][START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]. To address this issue, we followed the conotoxin assembly workflow outlined in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. Briefly, we first mapped reads back to our assembled contigs using the 'very sensitive local' and no discordant' options. Then, we identified conotoxins within our dataset by using blastn v2.2.31 (word size = 11, evalue = 1e-10) to associate our assembled contigs (from SPAdes) with conotoxins we targeted in the bait design. We generated a set of unique conotoxin 'seed . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; sequences' (a short stretch [~100bp] of conotoxin-blasted sequence) using a combination of of the pysam module (https://github.com/pysam-developers/pysam), cd-hit v4.6 (percent identity = 98%), cap3 (overlap percent identity cutoff = 99%), blastn v2.2.31 (word size = 11, evalue=1e-10), and Tandem Repeats Finder v4.09 (Benson 1999, minscore = 12, maxperiod = 2). We mapped reads to these seed sequences using bowtie2 v2.2.6 (very sensitive local and no discordant options enabled) and built out the conotoxin sequences using the PRICE v1.2 algorithm, which uses an iterative mapping and extension strategy to build out contigs from initial seed sequences [START_REF] Ruby | PRICE: software for the targeted assembly of components of (Meta) genomic sequence data[END_REF]. We ran price on each seed sequence at 5 minimum percentage identity (MPI) values (90%, 92%, 94%, 96%, 98%) with a minimum overlap length value of 40 and a threshold value of 20 for scaling overlap for contig-edge assemblies. A reassembled sequence was retained if it shared 90% identity with the original seed sequence and we reduced redundancy by only retaining the longest sequence per seed sequence out of the 5 MPI assembly iterations. This approach is described in further detail in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. We note that the final conotoxin sequences per sample consisted of exon fragments, where each sequence represents a single conotoxin exon flanked by any adjacent noncoding region.

We updated our conotoxin reference database because we targeted additional conotoxin transcripts from [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF].We used blastn v2.2.31 (word size = 11, evalue =1e-10) and EXONERATE v2.2.0 to define exon/intron boundaries for these additional conotoxin transcripts and added them to our conotoxin reference database. The final conotoxin reference database consisted of conotoxin sequences with the coding regions denoted and gene superfamily annotated. We also annotated the conotoxin sequences for functional region (e.g., signal, pre, . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; mature, post) using blastn v2.2.31 (word size = 11, evalue = 1e-10) with a conotoxin reference database that was previously categorized by functional region [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF].

With the final conotoxin reference database, we performed blastn v2.2.31 (word size = 11, evalue = 1e-10) searches between the conotoxin reference and every sample's re-assembled conotoxin sequences. We retained sequences if they could align across the entire coding region of the reference sequence. We guessed the coding region for each retained sequence by aligning the query sequence with the reference conotoxin using mafft v7.222 and denoting the coding region as the region of overlap with the exon in the reference conotoxin. We fixed misassemblies by mapping reads with bowtie2 (very sensitive local and no discordant options enabled, score min = L, 70, 1) back to each conotoxin assembly and marked duplicates using picard-tools v2.0.1. We masked regions below 5X coverage and discarded sequences if coverage was below 5X across the entire predicted coding region. To generate the final set of conotoxin sequences per sample, we merged sequences using cd-hit v4.6.4 (percent identity = 98%, use local sequence identity, alignment coverage of longer sequence = 10%, alignment coverage of short sequence = 50%).

Targeted sequencing experiment evaluation

We generated the following statistics to evaluate the overall efficiency of the capture experiment: (1) we calculated the % reads mapped to our targets by mapping reads to a reference containing all targets (both phylogenetic markers and conotoxin sequences) using bowtie2 v2.2.7 (very sensitive local and no discordant options enabled, score min = L, 70, 1), (2) we calculated the % duplicates that were identified through the picard-tools, and (3) we calculated average coverage across the phylogenetic markers and conotoxin sequences. We also evaluated the effect . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] of tissue quality (measured by the maximum fragment length of the extracted DNA sample via gel electrophoresis) and genus (only on Conus, Profundiconus, and Conasprella, the three genera with more than 1 sample included in this study) on these capture efficiency metrics using an Analysis of Variance (ANOVA). To assess the effectiveness of conotoxin sequence recovery, we compared our capture results with conotoxin diversity estimates from [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] and calculated the average change in those estimates.

Phylogenetic inference

In addition to the 362 samples that we sequenced in this study, we obtained sequences for 10 other species (Table 1). For two of these species, we used data from another targeted sequencing study (Abdelkrim et al. unpublished). We used blastn (word size = 11, evalue = 1e-10) to identify loci that were present in our phylogenetic marker reference. These sequences were filtered under conditions similar to the filtering strategy applied to the phylogenetic markers in this study. For the other eight species, we used data from venom duct transcriptomes (Safavi-Hemami et al. unpublished). With these transcriptomes, we trimmed data using trimmomatic v0.36 and merged reads using flash using parameters previously described above.

We assembled each transcriptome using Trinity v2.1.1 [START_REF] Grabherr | Full-length transcriptome assembly from RNA-Seq data without a reference genome[END_REF]) reduced redundancy in these transcriptomes with cap3 and cd-hit (percent identity = 99%). We used blastn (word size = 11, evalue=1e-10) to associate contigs with the phylogenetic markers present in our dataset. We used bowtie2 v2.2.7 (very sensitive local and no discordant enabled), samtools v1.3, and bcftools 1.3 to map reads and call SNPs. We removed sequences if they were below 4X coverage for > 30% of the sequence and masked bases if they were below 4X coverage. We also removed sequences if they had a heterozygosity value two standard deviations away from . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] the mean heterozygosity within a sample. We used to mafft v7.222 to align loci across a total of 373 samples.

We inferred phylogenies under both maximum likelihood [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF]) and coalescent-based methods [START_REF] Mirarab | ASTRAL-II : coalescent-based species tree estimation with many hundreds of taxa and thousands of genes[END_REF]. We used RAxML v8.2.9 [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] to generate a maximum likelihood phylogeny using a concatenated alignment under a GTRGAMMA model of sequence evolution and estimated nodal support via bootstrapping. We generated the coalescent-based phylogeny using ASTRAL-II v5.5.9 [START_REF] Mirarab | ASTRAL-II : coalescent-based species tree estimation with many hundreds of taxa and thousands of genes[END_REF] with individual locus trees generated under default parameters in RAxML v8.2.9. We estimated local posterior probabilities as a measure of branch support [START_REF] Sayyari | Fast Coalescent-Based Computation of Local Branch Support from Quartet Frequencies Article Fast Track[END_REF]. Due to the underperformance of the capture experiment, we ran both phylogenetic analyses with loci that had 80% of the taxa, 50% of the taxa, and 20% of the taxa. For each iteration, we removed taxa that had > 90% missing data.

Time calibration

We estimated divergence times using a Bayesian approach with MCMCTree implemented in PAML v4.9g [START_REF] Yang | PAML 4: phylogenetic analysis by maximum likelihood[END_REF]. Given the size of our alignments, we first estimated branch lengths using baseml and then estimated divergence times using Markov chain Monte Carlo (MCMC). We used a HKY85 + Γ substitution model and used an independent rates clock model. We left all other settings on default. We performed two independent runs of the analysis and checked for convergence among the runs. To account for uncertainty in branching order in our phylogeny, we executed dating analyses across all trees generated from RAxML.

For time calibration, we applied a maximum constraint of 55 million years at the root of Conidae, which corresponds with the first confident appearance of Conidae in the fossil record . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; [START_REF] Kohn | Tempo and mode of evolution in conidae[END_REF]. We assigned 12 additional fossils (Table S2, Fig. S1 [START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF][START_REF] Hendricks | The Genus Conus (Mollusca:Neogastropoda) in the Plio-Pleistocene of the Southeastern United States[END_REF][START_REF] Hendricks | Glowing Seashells : Diversity of Fossilized Coloration Patterns on Coral Reef-Associated Cone Snail (Gastropoda : Conidae) Shells from the Neogene of the Dominican Republic[END_REF][START_REF] Hendricks | Diversity and preserved shell coloration patterns of Miocene Conidae (Neogastropoda) from an exposure of the Gatun Formation, Colón Province, Panama[END_REF]) to nodes throughout the phylogeny as minimum age constraints, which MCMCtree treats as soft bounds on the minimum age [START_REF] Yang | PAML 4: phylogenetic analysis by maximum likelihood[END_REF]. Further information on fossil placement on nodes can be found in the Supplementary. A recent paper showed that the number of species in Lautoconus may be overestimated (Abalde et al. 2017). To account for potential artificial inflation in the species richness of this clade, we artificially removed half the unique species in Lautoconus from our dataset and ran all dating analyses and downstream diversification analyses on this secondary dataset.

Characterizing diversification patterns

To visualize lineage accumulation patterns, we generated a log-lineage through time plot using the R package APE [START_REF] Paradis | APE : Analyses of Phylogenetics and Evolution in R language[END_REF]. We estimated diversification rates and identified rate shifts using BAMM (Bayesian Analysis of Macroevolutionary Mixtures) [START_REF] Rabosky | Automatic detection of key innovations, rate shifts, and diversitydependence on phylogenetic trees[END_REF], which uses reversible jump Markov chain Mone Carlo to explore potential lineage diversification models. To account for non-randomness in species sampling across Conidae genera, we applied generic-specific sampling fractions. Using the number of valid Conidae names on WoRMS as estimates of total species diversity in each genus (Worms Editorial Board 2017), we applied a sampling fraction of 32.1% to Profundiconus, 50% to Lilliconus, 100% to Californiconus, 16.7% to Pygmaeconus, 28% to Conasprella, and 33.7% to Conus. We ran BAMM for 100 000 000 generations and assessed convergence by calculating ESS values. We analyzed and visualized results using the R package BAMMtools [START_REF] Rabosky | BAMMtools : an R package for the analysis of evolutionary dynamics on phylogenetic trees[END_REF].

Trait dependent diversification

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

We tested for the impact of evolvability (measured as conotoxin gene diversity) on diversification patterns using two trait dependent diversification methods, focusing on the genus Conus. We focused our hypothesis testing on Conus because conotoxin diversity is wellcharacterized in this group [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF]) and the sequence capture approach used in this study likely represents uniform sampling in conotoxin gene diversity across the genus. This is in contrast to other genera in Conidae, such as Conasprella or Profundiconus, where low conotoxin diversity values are likely the result of poor knowledge of the venom repertoire of these genera (Fig. S2) First, we used BiSSE (binary state speciation and extinction, [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF] implemeneted in the R package diversitree (FitzJohn 2012), which employs a maximum likelihood approach to estimate the impact of a binary trait on speciation, extinction, and transition rates between character states. We coded the conotoxin gene diversity data as 'low' or 'high' across several thresholds (i.e., 250, 300, 350, 400, 500, 550, or 600 estimated conotoxin genes per species) and compared BiSSE models where speciation rates were allowed to vary or remain equal between traits. We applied a sampling fraction of 33.7%, taking the maximum number of Conus species to be the number of valid names on WoRMS (World Register of Marine Species, (Worms Editorial Board 2017)). We determined the best-fitting model using Akaike Information Criterion (AIC). Second, we used FiSSE (Fast, intuitive State-dependent Speciation-Extinction analysis), a non-parametric statistical test that assesses the effects of a binary character on lineage diversification rates [START_REF] Rabosky | FiSSE : A simple nonparametric test for the effects of a binary character on lineage diversification rates[END_REF]. We followed the same coding strategy as in the BiSSE analyses to convert conotoxin gene diversity counts to binary character states. Finally, we used STRAPP (Structured Rate Permutations on Phylogenies, [START_REF] Rabosky | A Robust Semi-Parametric Test for Detecting Trait-Dependent Diversification[END_REF]) implemented in the R package BAMMtools [START_REF] Rabosky | BAMMtools : an R package for the analysis of evolutionary dynamics on phylogenetic trees[END_REF].

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

STRAPP is a semi-parametric approach that tests for trait dependent diversification by comparing a test statistic with a null distribution generated by permutations of speciation rates across the tips of the phylogeny [START_REF] Rabosky | A Robust Semi-Parametric Test for Detecting Trait-Dependent Diversification[END_REF]. We generated the empirical correlation (method = Spearman's rank correlation) between speciation rates and conotoxin gene diversity and compared this test statistic with the null distribution of correlations generated by permutations of evolutionary rates across the tree. We performed a two-tailed test with the alternative hypothesis that there is a correlation between speciation rates and total conotoxin gene diversity.

We also tested the impact of diet and larval dispersal strategy on diversification patterns.

Both piscivory and molluscivory is known to have evolved from the ancestral vermivory condition in cone snails [START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF]) and these diet transitions may be associated with increased diversification rates due to access to new dietary niches. In addition, differing larval dispersal strategies including long-lived larval stages (planktotrophy) and short-lived and/or direct developing larvae (lecithotrophy) are hypothesized to impact long term diversification patterns [START_REF] Jablonski | Larval ecology and macroevolution in marine invertebrates[END_REF]). We coded diet as either vermivory, molluscivory, and piscivory using natural history information from (Jiménez-Tenorio & Tucker 2013). We tested the impact of speciation and extinction using MuSSE (multistate speciation and extinction, (FitzJohn 2012))

where speciation rates were allowed to vary or remained equal among traits. We excluded species that were documented to feed on multiple diet types from this analysis. For larval type, we used protoconch morphology from (Jiménez-Tenorio & Tucker 2013) to infer larval dispersal strategy, where multispiral protoconchs were indicative of planktotrophic larvae. We tested the impact of larval type on diversification patterns using BiSSE and FiSSE.

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Results

Targeted sequencing data

We sequenced an average of 9,548,342 reads (range: 1,693,918 -29,888,444) across the 363 samples (Table S1). After redefining exon/intron boundaries in the phylogenetic marker reference, we ultimately targeted 2210 loci. On average, we recovered 1388 of these loci per sample (range: 30 -1849, Table S1) at an average coverage of 12.39X (range: 3.08X -27.87X, Table S1). For the conotoxin dataset, each sequence we re-assembled contained a single conotoxin exon with any associated noncoding regions (referred to here as 'conotoxin fragments'). We recovered on average 3416 conotoxin fragments per sample (range: 74 -11535 fragments, Table S1) at an average coverage of 32.3X (range: 5.06X -65.77X, Table S1). When mapped to a reference containing both the phylogenetic markers and conotoxin genes, the % reads mapped to our targets was on average 14.86% (range: 0.7% -38.07%, Table S1) and the average level of duplication was 47.47% (range: 22.89% -89.06%, Table S1).

We found that genus had an impact on % mapped and % duplication, where non-Conus genera had lower % mapping and lower % duplication (Fig. S2). These differences likely occurred because conotoxin fragments were not easily recovered in these genera (ANOVA, p < 0.0001, Fig. S2). Genus did not have an impact on coverage or the number of phylogenetic markers recovered (ANOVA, p > 0.05, Fig. S2). We found that tissue quality, measured by the maximum fragment length visualized via gel electrophoresis, had a significant impact on the capture efficiency metrics (ANOVA, p < 0.0001, Fig. S3). DNA samples with strong genomic bands at the top of the gel tended to have higher % mapping, less % duplication, higher coverage, and a greater number of targets recovered (Fig. S3).

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Our final conotoxin sequence dataset consists of exon fragments and we do not have information on exon coherence (which exons pair together on the same gene). We were unable to assemble full conotoxin genes because conotoxin introns are long (>1 kilobases, [START_REF] Wu | Molecular evolution and diversity of Conus peptide toxins, as revealed by gene structure and intron sequence analyses[END_REF]) and exceed the average insert size of our sequencing experiment (~350bp). We recovered fragments from all 58 gene superfamilies we targeted and obtained 159,670 sequences containing some or all of the mature toxin region (Table S3). Total conotoxin gene diversity per species (estimated by summing across all signal region exon fragments and sequences containing the entire coding region) ranged from 5 to 1280 copies in Conus, 31 to 88 copies in Profundiconus, and 7 to 164 in Conasprella (Table S1). Total conotoxin diversity was 311 copies for Californiconus californicus, 12 copies for Pygmaeconus tralli, and 30 copies for the outgroup taxon, Bathyoma sp (Table S1). When compared to samples in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF], the average change (increase or decrease) in total conotoxin gene diversity was ~90 gene copies (Table S4). If samples performed poorly in the number of phylogenetic markers recovered, conotoxin gene diversity estimates tended to be lower in this study than in [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]) and vice versa (Fig. S4). The average absolute change in the number of fragments recovered per gene superfamily by region was 3.7 for sequences containing the signal region, 12.2 for the prepro region, 9.6 for the mature region, 48.9 for the post region, and 3.4 for sequences containing the entire coding region (Table S5, Fig. S5). We note several key outliers: the average absolute change in the number of fragments was 104.3 for the T gene superfamily containing the prepro region, 210.4 for the O1 gene superfamily prepro region, 57.4 for the O1 gene superfamily mature region, 219.9 for the O2 gene superfamily mature region, and 1417 for the T gene superfamily post region (Table S5, Fig. S5).

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Phylogeny

The amount of missing data from the alignments was 15.4% when a minimum of 80% of the taxa were present in each locus, 26.8% when 50% of the taxa were present, and 38.6% when 20% of the taxa were present. The number of loci retained in the alignment was 387 (107,011 bp) when a minimum of 80% of the taxa were present in each locus, 976 (237,027 bp) when 50% of the taxa were present, and 1476 loci (336,557 bp) when 20% of the loci were present. Across all methods and datasets, we recovered phylogenies with a moderate level of resolution (average number of nodes resolved = 71.1%, range = 61.4 -79.2%, Table S6). In general, as increased amounts of sequence data was given to the phylogenetic programs, more nodes became resolved (Table S6). While we recovered all 6 genera within Conidae with high confidence, relationships among subgenera were less supported (bootstrap and PP = 100%, Fig. 1, Fig. S6, S7, S8).

Divergence time estimation

We found evidence for three major branching events during the Eocene: (1) a branching

Diversification patterns

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

We found that most branching events within each genus began to occur in the Miocene and continued until the present (Fig. 1). When analyzing the entire dataset, we found support for diversification rate heterogeneity, where BAMM identified at least one rate shift across Conidae (Fig. 1,S10). Across the 95% credible set of distinct shift configurations, BAMM detected an increase in diversification rates on the branch leading to Lautoconus, a clade consisting mainly of species endemic to the Cape Verde islands (Fig. 1,S10). However, when examining an artificially reduced dataset consisting of half the species within Lautoconus, we detect no rate shift or a decrease in diversification rates leading to the Conus clade (Fig. S11).

Trait dependent diversification

Across all thresholds for the BiSSE analysis, we found that diversification rates were not influenced by conotoxin gene diversity. In all cases, the null model was either preferred (delta AIC > 2, Table S7) or was indistinguishable from a model where speciation and extinction were allowed to vary (delta AIC < 2, Table S7). Both the FiSSE and STRAPP analyses revealed that speciation rates were not correlated with conotoxin gene diversity (p > 0.05). These results were consistent across both the full dataset and the reduced dataset.

We found that diversification rates were not dependent on diet when analyzing the full dataset (Table S8). However, in the reduced dataset, we found a signal for diet-dependent speciation rates (delta AIC > 2, Table S8). We found that species with mollusk-feeding diets had the fastest speciation rates (0.33), followed by piscivory (0.24), and vermivory (0.16). For the larval dispersal trait, we found support for trait-dependent speciation rates in the full dataset (delta AIC > 2, Table S9), where species with short-lived larvae had higher speciation rates (0.27 . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; vs. 0.16). However, this result was not significant when examining the reduced dataset (Table S9).

Discussion

Capture results

Our targeted sequencing experiment underperformed initial testing of this sequencing method on cone snails [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. Although tissue quality impacted capture metrics (Fig. S3), the % of reads mapping to our targets for even our best samples was ~30% lower than expected [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. While it is difficult to determine the exact cause of this depression in our capture statistics, we hypothesized that changes made in the bait design between this study and [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF] may have led to poorer capture results. For example, we recovered an overabundance of conotoxin sequences containing the post region from the T gene superfamily that has no clear co-variation pattern with phylogenetic relatedness (Fig. S12), which likely indicates a large amount of non-specific binding due to conotoxin misclassification. In the future, we suggest re-designing the baits to only include sequences from only the most critical regions (signal region and mature region) to avoid nonspecific binding. Although overall capture efficiency statistics were low, the absolute change in conotoxin diversity estimates per gene superfamily was generally minor (Table S5). Therefore, we do not believe that total conotoxin diversity metrics were severely biased by the sequencing method.

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Phylogenetic relationships

Below, we discuss the results of our phylogenetic analyses, how the phylogenetic relationships compare to past work, and their implications for Conidae taxonomy. Unless otherwise noted, the results we highlight below have at least 90% bootstrap support in the RAxML analyses and 90% posterior probabilities from the ASTRAL-II analyses (Figure S7,S8).

When present results on subgeneric relationships starting from the top of the tree shown in Figure S6.

We recovered all six major deep lineages representing genera in Conidae that were previously described in recent molecular phylogenetic studies using mtDNA (Puillandre et al. 2014a;Uribe et al. 2017), Fig. 1, S6, S7, S8). Specifically, we find strong support for Profundiconus, Californiconus, Lilliconus, Pygmaeconus, Conasprella, and Conus, as separate and distinct lineages. We also confirm the branching order of these six genera that were recently described using mtDNA genomes [START_REF] Uribe | Phylogenetic relationships of Conidae based on complete mitochondrial genomes[END_REF], with Profundiconus being sister to all other genera, Pygmaeconus + Lilliconus sister to Californiconus, Californiconus + Lilliconus + Pygmaeconus sister to Conasprella, and these four genera sister to Conus.

Based on the molecular phylogeny from three mtDNA genes, monophyletic groupings of species from Conasprella were classified into several subgenera (Puillandre et al. 2014a;b). We note several differences between past results and our study in the relationships among these genera and their monophyly:

(1) Ximeniconus is sister to all other Conasprella in some trees, or we reconstructed a polytomy at the base of Conasprella, which contrasts with Conasprella (Kohniconus) arcuata recovered at the base of Conasprella in previous work (Puillandre et al. 2014a).

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

(2) Kohniconus is polyphyletic. In (Puillandre et al. 2014a), only a single species from Kohniconus was included and we find evidence for the non-monophyly of Kohniconus when we included the additional species, C. centurio. Given these results, we propose that C. emarginatus, C. delssertii, and C. centurio be placed in the subgenus Kohniconus and C. arcuata placed in a new subgenus.

(3) Endemoconus is paraphyletic. When including an additional species (C. somalica) not sequenced in (Puillandre et al. 2014a), we find that Endemoconus is not monophyletic.

Based on these results, C. somalica should be transferred to Conasprella.

Within Conus, our results largely confirm previous findings that C. distans is sister to all other Conus species and the relationships among subgenera remain tenuous and difficult to resolve (Puillandre et al. 2014a). We note the following differences in subgenera relationships and classification between our results and past work:

(1) We found support the sister relationship between Turriconus and Stephanoconus, which has not been recovered in a previous study (Puillandre et al. 2014a).

(2) We found support for the monophyly of Pyruconus across our RAxML analyses, but not our ASTRAL-II analyses. The monophyly of Pyruconus was not supported in (Puillandre et al. 2014a).

(3) C. trigonus and C. lozeti were classified into the subgenus (Plicaustraconus) based on morphological characters [START_REF] Jiménez-Tenorio | Illustrated catalog of the Living Cone snails[END_REF]Puillandre et al. 2014b). We found this subgenus to be polyphyletic when sequence data was obtained.

(4) Similar to (Puillandre et al. 2014a), we found that Textila + Afonsoconus is sister to Pionoconus. However, instead of the unsupported relationship of Asprella as sister to these three subgenera, we found support for Gastridium as the sister group.

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

(5) We found support for the sister relationship between Asprella and Phasmoconus, which conflicts with the unsupported relationship shown (Puillandre et al. 2014a), where these subgenera branch in different parts of the phylogeny.

(6) We find support for the following successional branch order: Tesselliconus, Plicaustraconus, Eugeniconus, and Conus. We found that Conus is sister to Leptoconus, Darioconus, and Cylinder, but the relationships among these three subgenera remained unresolved. This conflicts with (Puillandre et al. 2014a) as Cylinder was paraphyletic, whereas in our results with increased sampling of Eugeniconus, Cylinder became monophyletic.

(7) We did not find strong support for the subgenus Calibanus, contrasting with previous work (Puillandre et al. 2014a). In our results, we found that C. thalassiarchus and C. furvus were not sister to each other, or their relationship resulted in an unresolved polytomy. Additional investigation into the subgeneric status of these two species. Here, we found that it was nested within Dauciconus. No other species within this subgenus have been sequenced up until this point. Therefore, we synonymize Atlanticonus with Dauciconus because C. granulatus is the type species for Atlanticonus.

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

(10) Two species (C. pergrandis and C. moncuri) sequenced in this study were placed into the subgenus Elisaconus (Puillandre et al. 2014b). Our results do not support the monophyly of Elisaconus, as the sister relationship between C. moncuri and C. pergrandis was not supported in 5/6 trees. Additional data is required to classify C. moncuri and C. pergrandis into the appropriate subgenus.

(11)

C. cocceus was placed into Floraconus based on morphological characters in (Puillandre et al. 2014b). With sequence data, we found that it was actually nested within

Phasmoconus. Therefore, we transfer C. cocceus to the subgenus, Phasmoconus.

Classification within Conidae is known to be highly unstable (Jiménez-Tenorio & Tucker 2013; Puillandre et al. 2014a;b;Puillandre & Tenorio 2018). Although the phylogeny presented here improved understanding of subgeneric relationships and monophyly of subgenera, resolving relationships within Conidae still remains a significant challenge. Given the underperformance of our capture experiment (Table S1), it is unclear if the reason for the moderate power in resolving relationships is due to insufficient data/incomplete data or due to short internal branches during the origination of Conidae subgenera that are extremely difficult to resolve. Overall, our results suggest that both additional data and increased sampling of Conidae species are reasonable pursuits to continue attempting to resolve the phylogeny and classification of this family of marine snails.

Timing of diversification

The timing of splits between major are largely congruent with past estimates from a study using mtDNA genomes [START_REF] Uribe | Phylogenetic relationships of Conidae based on complete mitochondrial genomes[END_REF], Fig. 1,S9). However, our age estimates for the . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] branching events between Californiconus, Lilliconus, and Pygmaeconus are much younger (occurring across the Oligocene into the Miocene) than previous estimates (occurring across the Eocene into the Oligocene, [START_REF] Uribe | Phylogenetic relationships of Conidae based on complete mitochondrial genomes[END_REF], Fig. 1,S9). This discrepancy may have been caused by differences in fossil calibration, as we included many more fossils in this study compared to previous studies. The Conidae fossil record and analyses of several molecular phylogenetic studies suggest a major radiation of Conus during the Miocene [START_REF] Kohn | Tempo and mode of evolution in conidae[END_REF][START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF][START_REF] Uribe | Phylogenetic relationships of Conidae based on complete mitochondrial genomes[END_REF]. While we noted that many branching events within Conus ocurred during the Miocene into the present, we did not detect an increase in diversification on the branch leading to the origin of Conus (Fig. 1,S10,S11). This is congruent with diversification rates estimated from the fossil record [START_REF] Kohn | Tempo and mode of evolution in conidae[END_REF], suggesting that the accumulation of species during the Miocene may have been a function of an increased number of lineages present rather than an increase in diversification rates. The number of species we included in the subgenus Lautoconus had an impact on the BAMM diversification analyses. On the full dataset, BAMM detected an increase in diversification rates leading to Lautoconus (Fig. 1, S10), a known and documented radiation of cone snails [START_REF] Duda | Explosive radiation of Cape Verde Conus, a marine species flock[END_REF][START_REF] Cunha | Patterns of cladogenesis in the venomous marine gastropod genus Conus from the Cape Verde islands[END_REF]. However, when we remove half the species in response to recent work suggesting taxonomic inflation in this subgenus (Abalde et al. 2017), we do not detect the same shift.

Rather, there is partial support for no shift across Conidae, or a slight decrease in diversification rates leading to Conus (Fig. S11). These results suggest that the original diversification analyses and identified radiation of Lauotoconus may have been due to taxonomic inaccuracies biasing the diversification analyses results, rather than a true radiation. What is even more striking about these results is that we found minimal diversification rate heterogeneity across Conidae, despite the expansive species richness across this group. It is unclear whether this signal is real, or due to . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] other technical artifacts. For sample, although we included over 300 species in this study, this only represents ~30% of the total diversity in this group and may have hindered our ability to effectively estimate diversification rates. Similarly, new Conidae species are continually described, with over 100 species described over the last few years (Worms Editorial Board 2017). Therefore, our inability to estimate the number of living taxa may have weakened our ability to test the impact of diversification on this group.

Speciation rates and conotoxin gene diversity

Contrary to macroevolutionary expectations, we were unable to detect any relationship between speciation rates and conotoxin gene diversity across all trait dependent diversification analyses (Fig. 1, S11, Table S7). Even when performing the analyses with BiSSE, a method in recent years that has become the subject of criticism due to high false positive rates (Abosky 2017; Rabosky & Goldberg 2017), our analyses did not detect an impact of conotoxin gene diversity on diversification rates (Table S7). These results may have been expected, given that we found minimal levels of diversification rate heterogeneity in Conidae (at minimum, one shift, Fig. 1, S10, S11). As discussed previously, taxonomic instability in this group may have hindered our efforts to estimate past historical diversification patterns. However, we did find some signal for the impact of diet and larval dispersal strategy on diversification rates when using the BiSSe and MuSSE methods (Table S8,S9). Further work is needed to be fully confident in this signal given high false positive rates in these methods (Abosky 2017; Rabosky & Goldberg 2017) and given that our results depended on which dataset was used.

What is remarkable about these results is the lack of any signal on the impact of venom gene diversity on diversification rates in cone snails, even as we found some signal for trait-. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with[END_REF] dependent diversification in other Conidae characters. If this lack of signal is real, several biological factors may explain this decoupling between conotoxin gene diversity and speciation rates. A critical assumption in Conus biology is that ecological diversification driven by diet specialization is a major factor governing diversification dynamics in cone snails [START_REF] Duda | Molecular genetics of ecological diversification: duplication and rapid evolution of toxin genes of the venomous gastropod Conus[END_REF][START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF]. Past studies have shown that cone snail venom repertoires track their dietary breadth, providing a link between diet and venom evolution [START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF][START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. However, it is unclear whether or not the relationship between diet and venom evolution leads to ecological speciation due to divergence in prey preference. Ecological speciation is often difficult to detect in marine ecosystems and long-term diversification patterns may be better explained by traits that limit dispersal and promote isolation (Bowen et al.). Another possibility is that conotoxin phenotypic divergence may not be the rate-limiting factor in prey specialization and divergence [START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF]. Conotoxin genes are under continuous positive selection and gene duplication that allow venom components to change rapidly in response to the environment [START_REF] Duda | Molecular genetics of ecological diversification: duplication and rapid evolution of toxin genes of the venomous gastropod Conus[END_REF][START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF][START_REF] Chang | Extensive and continuous duplication facilitates rapid evolution and diversification of gene families[END_REF][START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. This persistent evolutionary change in the venom cocktail suggests that perhaps venom evolution is not necessarily the factor limiting dietary shifts among species and ultimately, speciation among taxa. Ecological opportunity is hypothesized as a necessary component for diversification [START_REF] Losos | Adaptive radiation, ecological opportunity, and evolutionary determinism[END_REF] and may be a more critical factor limiting Conidae diversification. Indeed, evidence from the fossil record and past Conidae molecular phylogenetic studies indicate a concentration of lineage formation during the Miocene [START_REF] Kohn | Tempo and mode of evolution in conidae[END_REF][START_REF] Duda | Origins of diverse feeding ecologies within Conus, a genus of venomous marine gastropods[END_REF][START_REF] Uribe | Phylogenetic relationships of Conidae based on complete mitochondrial genomes[END_REF], a period that is coincident with the formation of coral reefs in the Indo-Australian Archipelago [START_REF] Cowman | Coral reefs as drivers of cladogenesis: expanding coral reefs, cryptic extinction events, and the development of biodiversity hotspots[END_REF]. Our results also show a concentration of branching events during this period as well, though we do . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; not detect a shift in diversification rates (Fig. 1, S10, S11). Overall, our results point to increased taxonomic sampling and a holistic approach to investigating factors shaping diversification in Conidae for future work.

Venom evolution is assumed to be a key innovation that led to the evolutionary success of venomous animal lineages [START_REF] Pyron | EXTINCTION[END_REF][START_REF] Sunagar | Deadly Innovations: Unraveling the Molecular Evolution of Animal Venoms[END_REF]) and a large body of work is devoted towards understanding how venom evolves and responds to the environment over time [START_REF] Kordis | Adaptive evolution of animal toxin multigene families[END_REF][START_REF] Wong | Venom evolution through gene duplications[END_REF][START_REF] Bowen | Complex cocktails: the evolutionary novelty of venoms[END_REF]. However, the impact of venom evolution on higher-level diversification patterns is rarely tested. Here, we examined the effect of variation in the adaptive capacity of venom across Conidae species and found it had no influence on macroevolutionary diversification patterns. Although we do not detect a strong signal of conotoxin gene diversity shaping speciation rates in Conidae, it does not refute the importance of venom evolution in adaptation and prey specialization as venom may be necessary, but not sufficient, to promote speciation [START_REF] Duda | Geographic variation in venom allelic composition and diets of the widespread predatory marine gastropod Conus ebraeus[END_REF][START_REF] Safavi-Hemami | Specialized insulin is used for chemical warfare by fish-hunting cone snails[END_REF][START_REF] Chang | Age-related association of venom gene expression and diet of predatory gastropods[END_REF][START_REF] Phuong | Dietary breadth is positively correlated with venom complexity in cone snails[END_REF][START_REF] Phuong | Targeted Sequencing of Venom Genes from Cone Snail[END_REF]. Future work in other venomous animal systems may shed light on whether or not the ability to adapt to different prey through venom evolution translates to the long-term evolutionary success of taxa. Table S1. Sample information and capture efficiency metrics. We first list the species name listed in all phylogenetic analyses ("species") and the accepted taxonomic classification in the WoRMS database at the genus, subgenus, and species level ("WoRMS genus", "WoRMS subgenus", "WoRMS species"). We then list the specimen ID ("ID"), the collection source ("Collection"), the year the sample was collected ("Year Collected"), how the sample was preserved ("Preservation type"), and the country the sample originated from ("Country"). We then list "Gelsize", or the largest fragment size visualized via gel electrophoresis as a way to measure tissue quality. Values can be either "g" for genomic band, or 1500, 1000, or 500 for bands beginning at 1500bp, 1000bp, or 500bp. We list the data collection method ("Data . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; collection method"), the number of reads sequenced ("# of reads sequenced") and several capture efficiency metrics. Finally, we list the total estimated number of conotoxin genes per species ("Conotoxin gene diversity").

Table S2. Information on fossils used for calibration. We list the fossil taxon ("Fossil Species"), the genus ("Clade assignment"), extant species related to the fossil ("Compared with"), the formation ("Formation"), the age of the fossil ("Age"), and its citation ("Reference").

Table S3. Table showing the number of conotoxin sequences recovered per species for each gene superfamily. Within each gene superfamily, conotoxin sequences were categorized based on whether they contained the entire coding region or mostly the signal, prepro, mature, or post regions.

Table S4. Comparison of conotoxin gene diversity estimates between this study and [START_REF] Phuong | Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and venom diversity[END_REF]. These represent comparisons between technical replicates (capture experiment was performed on the same libraries in both studies).

Table S5. Comparison of conotoxin gene diversity estimates between this and (Phuong &

Mahardika 2017), broken down by gene superfamily. Within each gene superfamily, conotoxin sequences were categorized based on whether they contained the entire coding region or mostly the signal, prepro, mature, or post regions. These represent comparisons between technical replicates (capture experiment was performed on the same libraries in both studies).

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; Table S6. Number of nodes resolved depending on the amount of missing data and the tree inference method. Phylogenetic trees were inferred using either RAxML or ASTRAL-II. The "% taxa per locus was" the percent of samples needed per locus in order to retain the locus for phylogenetic inference.

Table S7. Venom gene diversity BiSSE AIC results. "Dataset" represents whether the full dataset was used or the reduced dataset. "Threshold" represents the conotoxin gene diversity value used to decide between "high" and "low" conotoxin diversity. Values above the threshold value were categorized as "high" and values below were categorized as "low". "AICvariable rates" shows AIC values for a model where speciation and extinction rates were allowed to vary depending on a trait. "AIC -equal rates" represents AIC values for the null model, where rates were not allowed to vary by trait. Table S8. Diet BiSSE AIC results. Model values were generated under a variable rates model (where speciation was allowed to vary) or under an equal rates model (speciation rates across trait states were equal).

Table S9. Larval dispersal type BiSSE AIC results. Model values were generated under a variable rates model (where speciation was allowed to vary) or under an equal rates model (speciation rates across trait states were equal). was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; of a concatenated alignment where loci were kept if at least 20% of species was present in the locus. Best tree is shown and erroneous and intraspecific tips were pruned. . CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018; represent faster speciation rates than cooler, blue colors. We note that in all shift configurations, there is a shift in diversification rates in the clade leading to Lautoconus.

  event leading to Profundiconus (56.5 mya, CI = 46.3 -65.3 mya, Fig. 1, S9), (2) a branching event leading to Conus (54.7 mya, CI = 42.5 -63.6 mya, Fig. 1, S9), and (3) a branching event separating Conasprella and Californiconus, Lilliconus, and Pygmaeconus (46.0 mya, CI = 36.5 -53.2 mya, Fig. 1, S9). The branching event leading to Californiconus occurred during the Oligocene (26.1 mya, CI = 13.8 -36.5 mya, Fig. 1, S9) and the split between Lilliconus and Pygmaeconus occurred during the Miocene (17.8 mya, CI = 9.25 -25.1 mya, Fig. 1, S9).

( 8 )

 8 C. sanderi was classified into its own subgenus (Sandericonus) based on morphological characters (Jiménez-Tenorio & Tucker 2013; Puillandre et al. 2014b). Here, when sequence data were obtained, we found it nested within Dauciconus. Therefore, we synonymize Sandericonus with Dauciconus because C. sanderi is the type species for Sandericonus. (9) C. granulatus was classified into its own subgenus (Atlanticonus) based on morphological characters (Jiménez-Tenorio & Tucker 2013; Puillandre et al. 2014b).

Figure 1 .

 1 Figure 1. Time calibrated maximum likelihood phylogeny of the cone snails. Phylogeny was

Figure S1 .

 S1 Figure S1. Node placement of fossils. Numbers correspond to node placement justification in

Figure S2 .

 S2 Figure S2. Boxplots showing impact of phylogeny (categorized by Conidae genus) on capture

Figure S3 .

 S3 Figure S3. Boxplots showing impact of tissue quality (estimated by maximum DNA fragment

Figure S4 .

 S4 Figure S4. Scatterplot showing relationship between the number of phylogenetic markers

Figure S5 .

 S5 Figure S5. Histograms showing absolute change in conotoxin sequence diversity per gene

Figure S6 .

 S6 Figure S6. Maximum likelihood phylogeny inferred using RAxML, where 20% of the taxa

Figure S7 .

 S7 Figure S7. Phylogenies inferred through the coalescent-based method, ASTRAL-II. Individual

Figure S8 .

 S8 Figure S8. Maximum likelihood phylogenies generated using a concatenated alignment. Nodes

Figure S9 .

 S9 Figure S9. Maximum likelihood phylogeny dated with 13 fossil node calibrations in MCMCtree.

Figure S10 .

 S10 Figure S10. 95% credible set of distinct shift configurations from BAMM for the full dataset.

Figure S11 .

 S11 Figure S11. 95% credible set of distinct shift configurations from BAMM for the reduced

Figure S12 .

 S12 Figure S12. Diversity estimates for the A gene superfamily signal region and the T gene

Acknowledgements

. CC-BY-NC 4.0 International license It is made available under a was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint (which . http://dx.doi.org/10.1101/359976 doi: bioRxiv preprint first posted online Jul. 1, 2018;

Data availability

Raw read data will be made available at the National Center for Biotechnology Information Sequence Read Archive. Bait sequences, conotoxin sequences, scripts, and final datasets used for analyses will be uploaded onto Dryad following publication.

A majority of the sampling material in this paper originates from numerous shore-based expeditions and deep sea cruises, conducted respectively by MNHN and Pro-Natura International (PNI) as part of the Our Planet Reviewed programme (SANTO 2006, ATIMO VATAE, MAINBAZA, INHACA 2011, GUYANE 2014, PAPUA NIUGINI, KAVIENG 2014), by MNHN