Taxon Set and taxonomic data
The class Chondrichthyes is composed of two subclasses, the Holocephali (chimaeras) and the Elasmobranchii (sharks and rays), and includes 14 orders, 60 families, 198 genera and 1192 species (Supplementary table 1). Systematic relationships within Chondrichthyes, as with other taxa, are in flux[1–4], so we used the most recent combination of taxonomy and phylogeny to identify our taxon set (Chondrichthyan Tree of Life; downloaded October 15, 2015). Due to the inherent prolonged time required for analyses of this nature, taxonomic revisions often occur prior to completion and publication[5–9]. To aid readers in navigating recent changes that we could not incorporate into the analysis, we provide annotations to our master taxonomy highlight recent taxonomic revisions (Supplementary Table 1, column G) and include recently described species that have not been assessed and could not be included in this study (Supplementary Table 6). The subclass Holocephali includes one superorder (Holocephalimorpha) containing one order of Chimaeriformes (chimaeras); three families, six genera, and 49 species. The subclass Elasmobranchii includes three superorders: Batoidea, Galeomorphii, and Squalomorphii[10]. The superorder Batoidea includes four orders Myliobatiformes (stingrays), Rajiformes (skates), Rhinopristiformes (guitarfishes, wedgefishes, and sawfishes), and Torpediniformes (electric and thornback rays); 23 families, 86 genera, and 639 species. The sharks comprise two superorders: Galeomorphii and Squalomorphii. The superorder Galeomorphii includes four orders: Carcharhiniformes (ground sharks), Heterodontiformes (bullhead sharks), Lamniformes (mackerel sharks), and Orectolobiformes (carpet sharks); 23 families, 75 genera and 347 species. The superorder Squalomorphii includes five orders Hexanchiformes (cow sharks), Pristiophoriformes (saw sharks), Echinorhiniformes (bramble sharks), Squaliformes (dogfish sharks), and Squatiniformes (angel sharks); 11 families, 31 genera, and 157 species. The taxonomic hierarchy described above comprises the taxonomic data that we used to place and constrain those taxa without DNA sequence data (see Supplementary Table 1 and “Taxon-complete Analyses” subsection below).
DNA data matrix
We assembled a DNA data supermatrix from pre-existing GenBank and Barcode of Life Data System records (downloaded on or before September 15, 2014) as well as 54 novel sequences generated for this study. Data from GenBank can present particular challenges, outlined in Naylor et al. 2012[1], thus all sequence and species validity was checked prior to analysis. Of particular concern with GenBank sequence data is the potential for misidentified specimens leading to erroneous placement. As a check for this an initial set of trees were generated using RAxML[11] and topology was hand checked to verify reasonable placement of species included in our matrix. The matrix is composed of a novel set of 15 coding and non-coding regions as follows: 2 non-protein-coding mitochondrial loci (12S and 16S rDNA, 2,037 bp), 11 protein-coding mitochondrial loci (CO1, CO2, CO3, Cyt b, ND1, ND2, ND3, ND4, ND4L, ND5, and ND6; 10,341 bp), and 2 nuclear protein-coding loci (RAG1, 2,538 bp; SCFD2, 582 bp; Supplementary table 2). The novel sequence data generated for this study include 8 CO1, 1 Cyt b, 1 ND2, 1 ND4, and 42 RAG1 sequences. The supermatrix included representatives from all 14 orders, 59/60 families (98%), 173/198 genera (88%), and 642 species (645 originally, of which three were subsequently synonymized with valid names and removed) out of 1192 species (54%). At its maximum extent, the DNA data matrix comprises 15,498 bp; however, the alignment is sparse and taxonomic coverage averages 30% across loci (mitochondrial loci: 13–81%; nuclear loci 12–27%; Supplementary Table 2).
We used MAFFT v.7.221[12–14] to conduct local alignments for each locus. Nuclear and mitochondrial protein-coding sequences are straightforward to align, but mitochondrial non-protein-coding sequences are subject to high frequencies of insertions and deletions (indels). Therefore, we aligned the indel-rich mitochondrial 12S and 16S non-protein coding sequences in two stages: first, we aligned the sequences by taxonomic order, and second, we combined the resulting order-specific alignments and realigned the entire set together. We removed start and stop codons from aligned protein-coding sequences prior to testing nucleotide substitution models. We used JMODELTEST 2.0[15,16]and AIC criterion to identify the best-fit nucleotide substitution model for each locus. Three closely related best-fit models were identified: GTR + Γ (nuclear RAG1), GTR + Γ + Ι (mitochondrial 12S, 16S, CO1, CO2, Cyt b, ND1, ND2, ND3, ND4, ND4L, ND5 and ND6) and SYM + Γ + Ι (mitochondrial CO3 and nuclear SCFD2). Invariant site models, such as GTR + Γ + Ι and SYM + Γ + Ι, have been criticized because the proportion of invariant sites and the gamma shape parameter cannot be optimized independently. As a consequence, it is impossible to obtain reliable estimates for these parameters simultaneously[17]. The SYM model is a constrained, nested version of the GTR model.
We used RAxML[11] to infer individual gene trees and species trees based on partitioned, concatenated analyses. RAxML uses a computationally efficient version of GTR (GTRCAT) that accommodates rate heterogeneity and GTRCAT is a good fit for our small set of best-fit models. In the final step, the GTRCAT model optimizes parameters and calculates likelihood under GTR + Γ[11]. We conducted 1,000 bootstrap replicates in each analysis and inspected the resulting topologies for consistency and bootstrap support before proceeding with partitioned, concatenated analyses. We used a pragmatic and iterative approach to partitioning, which was informed by constraints on a small, sparse dataset and trade-offs associated with variation in nucleotide substitution rate and process. Nuclear sequences typically evolve slowly relative to mitochondrial sequences, and the two nuclear protein-coding sequences (RAG1 and SCFD2) were assigned to separate partitions. Although the mitochondrion is inherited as a single locus, its protein-coding and nonprotein-coding sequences exhibit different nucleotide substitution rates and indel frequencies. As a consequence, mitochondrial nonprotein-coding and protein-coding loci were assigned to two separate partitions. Mitochondrial protein-coding sequences are subject to codon-position-specific rate variation. One consequence of this rate variation is that the faster evolving mitochondrial 3rd codon position nucleotides (M3CPN) can become saturated over long periods of evolutionary time, and there is precedence for excluding M3CPN from phylogenetic reconstructions for ancient clades[18,19]. As part of our iterative approach to partitioning, we first conducted RAxML analyses with 1,000 bootstrap replicates and then used ROGUENAROK[20] to identify rogue taxa that erode bootstrap support. In ROGUENAROK analyses we specified the parameters as follows: Threshold: 50% majority rule consensus, Optimize: support, and Max Drop Set: 3. Additionally, we employed a raw-improvement-score threshold (0.5), which corresponds to a relatively large improvement of overall support, to identify and exclude rogue taxa[20].
When we excluded the M3CPN, RAxML analyses generated trees with low-overall bootstrap support, and two iterations of rogue identification yielded 60 rogue taxa (9.3%). Despite concerns over saturation, we tried including the M3CPN in a standard, uniform block alignment. This resulted in a substantial increase in overall bootstrap support, and two iterations of rogue identification yielded just 22 rogue taxa (3.4%). However, bootstrap support remained low at relatively deep nodes (among orders and among families within orders) within the superorder Batoidea. Batoidea comprises 4 orders and 23 families: Myliobatiformes (10 families), Rajiformes (3 families), Rhinopristiformes (5 families), and Torpediniformes (5 families). As a consequence, the impact of low bootstrap support on topology is large. In an attempt to reduce the potential negative impacts of saturation and generate increased bootstrap support for the deeper nodes within Batoidea, we implemented a novel, staggered-by-order alignment approach for the M3CPN. Instead of using a single uniform block alignment, we extracted all of the M3CPN and placed them in a separate partition. This resulted in two partitions for the mitochondrial protein-coding loci, one with first- and second-codon position nucleotides and another with third-codon position nucleotides. The partition containing first- and second-codon position nucleotides remained as a standard uniform block alignment. For the partition containing the 3rd codon position nucleotides, we staggered the alignment by taxonomic order. The consequence of this novel partitioning/alignment approach is that the M3CPN sequence data can only speak to affinities within, but not between, orders. This approach resulted in increased bootstrap support within Batoidea without loss of support in other parts of the phylogeny. Two iterations of rogue identification yielded 21 rogue taxa (3.3%). Importantly, 11 taxa, including the worst offenders in both analyses including M3CPN, were included in the shared subset of rogue taxa. Using the staggered-by-order approach for the M3CPN alignment and two iterations of rogue identification, we identified and excluded 21 rogue taxa from 15 genera: Bathyraja (n=2), Carcharhinus (2), Centrophorus, Dipturus (2), Discopyge, Isogomphodon, Leptocharias, Mustelus, Narke (2), Orectolobus, Potamotrygon (3), Raja, Spiniraja, Squalus, and Squatina (Supplementary Table 1). There were still many nodes with bootstrap support < 70% and these were collapsed prior to incorporating unresolved taxa, when the identified rogue taxa were also reintroduced. This means that there was likely very little effect of rogue selection on overall topology. While the novel staggered-by-order approach employed here with the M3CPN partition appears promising, it warrants further study.
Temporal Calibration
Ideally, calibration fossils should be subjected to a formal phylogenetic analysis or exhibit diagnostic apomorphies[21]; unfortunately, relatively few fossils assigned to Chondrichthyes meet these criteria[22]. There are three exceptions: (1) Chondrenchelys problematicus, which has affinities with stem Holocephalii (Chimaeriformes)[23], (2) Tingitanius tenuimandibulus, which has affinities with stem Platyrhinidae (thornback rays)[24], and (3) Protospinax annectans, which has affinities with the stem of the superorder Squalomorphii, a clade including Hexanchiformes, Pristiophoriformes, Squaliformes, and Squatiniformes[25]. We identified 7 additional calibration fossils that are distributed across the Chondrichthyan phylogeny; several of these are represented by substantial articulated remains (Supplementary Table 3). Chondrenchelys problematicus is the only formally-vetted calibration fossil and, importantly, it provides a hard minimum bound for the chondrichthyan root node (333.56 MY)[22]. The root node of Chondrichthyes is further characterized by a soft maximum bound of 422.4 MY (see Benton et al., 2015 for justification). Following Ho and Phillips[26], we used these hard minimum and soft maximum bounds to specify a lognormal calibration density and selected a mean that bounded 95% of the probability density within the 88.84 MY interval between the hard minimum and soft maximum bounds and a standard deviation that split the probability density evenly across the midpoint (377.98 MY) of this interval. While the ages of the 9 other calibration fossils provide hard minimum bounds for their calibrated nodes, there is not sufficient information to generate a calibration density for any of them.
We used treePL 1.0 in Ubuntu 14.04[27], which implements a flexible rate-smoothing algorithm, to assign a timeline of diversification to the phylogeny. Given the rate-smoothing behavior of treePL and the reported low substitution rate in chimaeras[28], we expect the actual crown Chimaeriformes to be older, and the Elasmobranchii crown age to be younger, than those reported here. Given uncertainty in the precise phylogenetic affinities of at least 7 of the 9 additional fossils, we chose to conduct two sets of dating analyses, one that included C. problematicus only, and another that included all 10 calibration fossils (Supplementary Table 3). For both of these calibration scenarios we generated a random sample of 500 root-node ages from the lognormal calibration density that we constructed for the root node (together measuring the root age we report), and then conducted separate treePL analyses using each of the 500 root-node ages. Our treePL analyses proceeded in two stages. In the first, we performed cross validation analyses (“cv” and “randomcv” commands) and tested performance of the available optimization routines (“prime” command). In the second, we incorporated control options (“thorough” command) to ensure that the preferred optimization routine converged. The 500 resulting treePL-dated topologies for each of the two calibration scenarios were subsequently used in the taxon complete analyses described below. Importantly, the topology recovered from our RAxML analyses of the concatenated data matrix, remained fixed across all of these treePL analyses; only the timeline of diversification changed between the two calibration scenarios and across the 500 root-node ages. The treePL output was converted into a single 500-tree Newick file, processed with TreeAnnotator 1.7.5 [29] using the default settings (no burnin; posterior probability limit 0.0; maximum clade credibility tree; median node heights).
Taxon-complete trees We added taxa without DNA sequence data to each of the 500 treePL-dated, molecular trees (“stage 1” trees) and then used a taxon-addition and polytomy-resolver algorithm (modified from Kuhn et al. 2011; details below) to generate a large distribution of fully resolved, taxon-complete candidate phylogenies. We used two taxonomic sources: the Chondrichthyan Tree of Life and the International Union for the Conservation of Nature Red List. Our distribution of taxon-complete trees includes three types of species[31]. Type 1 species have genetic data and are represented in the stage 1 trees. Type 2 species have no genetic data (or were identified as rogue taxa), but have at least one congener represented in the stage 1 trees. Type 3 species have no genetic data (or were identified as rogue taxa), and have no congeners in the stage 1 trees. Type 2 and Type 3 species were allowed to populate particular clades using taxonomic information and the topology (via node identities) of the stage 1 trees. We outline the rules we used to include taxa without any sequence data below.
Type 1 species were anchored relative to one another as resolved in the stage 1 trees, and we used 70% bootstrap-support (BS) as a threshold to topologically constrain inferred nodes during taxon addition. There were four scenarios in which the stage 1 trees needed to be modified to enforce genus, family, or order monophyly (Supplementary Table 4).
• First, there were nine genera with relatively weak support (BS < 70%; range: 7-68%) for genus monophyly within a highly supported clade (BS ≥ 70%; Supplementary Table 4). For these nine genera instances we pruned ten type 1 taxa from the stage 1 trees, reducing its size from 620 to 610 sp. We reincorporated these ten pruned taxa subsequently as type 2 species by constraining them to their named clades.
• Second, there were four families (Anacanthobatidae, Hemigaliidae, Somniosidae, and Triakidae) where there was weak support (BS < 70%; range: 33-62%) for family non-monophyly. In these four instances we collapsed the weakly supported nodes and subsequently reconstituted clades to enforce family monophyly[31] (Supplementary Table 4).
• Third, there were 24 mixed-genus and/or mixed-family clades with strong evidence (BS ≥ 70%) against monophyly for at least one genus or family. These “mixed clades” took on a variety of forms, from simple paraphyly to complex interdigitation of sub-genera or sub-families. We enforced genus monophyly for any genus or family within these mixed clades, unless there was strong evidence (BS ≥ 70%) against monophyly.
• Finally, for consistency between taxa with and without genetic data, we assumed that all genera, families, and orders were monophyletic unless there was strong evidence (BS ≥ 70%) against monophyly in the stage 1 trees. This rule affected one subgenus (Galeus minor clade), nine genera (Atlantoraja, Centrophorus, Chiloscyllium, Halaelurus, Mobula, Rajella, Rhinobatos, Sphyrna and Squalus), one mixed-genus clade (Dentiraja, Dipturus, Spiniraja together with Zearaja), one family (Pristidae) and one order (Orectolobiformes) in our stage 1 trees – each of these had only weak evidence (BS < 70%; range 32-69%) for monophyly.
After using these rules to modify the stage 1 trees, we imposed topological constraints on the placement of the remaining type 2 and type 3 species, including the 21 rogue taxa. Each type 2 species was restricted to its genus or its mixed-genus clade. There were four genera (Dasyatis, Galeus, Himantura, and Triakis) with strong evidence (BS ≥ 70%) against monophyly in the stage 1 trees that also required the addition of type 2 species. For these four genera, type 2 congeners were added to the largest candidate sub-clade for the genus. Each type 3 species was restricted to its named genus, and the entire genus was constrained in its placement among other genera according to higher-level (supergenus, family or order) taxonomic information (Supplementary Table 1). Twenty-three of the 198 recognized genera were not represented in the stage 1 trees, and these 23 genera were restricted to 6 of the 14 orders and 18 of the 60 families (Supplementary Table 1). There were 9 type 3 species not assigned to a family on http://sharksrays.org/. In these instances, we referred to the IUCN http://iucn.org/ for the original family-level designation (Supplementary Table 1).
The 500 dated taxon-complete trees were then each resolved using Polytomy Resolver[30], modified to allow for the partial constraints enumerated above. The Polytomy Resolver algorithm uses a customized R-script to generate an input file for BEAST 1.x, based on the original dated phylogeny, and the taxonomy additions outlined above. The generated input file leverages BEAST’s ability for “prior only sampling”, combined with a series of hierarchical topology constraints—both time-based and monophyly-based—that define the original tree topology and allowable taxon additions, to sample taxon-complete trees across tree-space. For both fossil calibration scenarios mean Growth Rate (birth – death) was unconstrained, while we set uniform flat priors for relative Death Rate (death/birth; one fossil calibration: 0.4 – 0.85; ten fossil calibration: 0.25 – 0.75). Each taxon infilling scenario was run in BEAST 1.7.5[29] for 3 million generations including a 1 million generation burn-in. Samples were drawn every 100,000 generations to avoid temporal autocorrelation between draws, and the resulting 20 trees from each of 500 scenarios were collated into a pseudo-posterior distribution of 10,000 fully resolved trees.