Distinct Responses of Elasmobranchs and Ray-Finned Fishes to Long-Term Global Change

Both biotic and abiotic factors likely played a role in influencing the diversification patterns of clades. Although the role of environmental forcing on the long-term evolution of biodiversity has been explored for invertebrate clades, little is known about how vertebrate groups responded to environmental changes. Among vertebrates, fishes (ray-finned fishes and elasmobranchs) have a long, rich, and complex evolutionary history comprising numerous diversification and extinction events. Yet, knowledge on the causes for the diversity fluctuations of these most speciose aquatic vertebrate clades in modern marine and continental ecosystems were restricted to qualitative interpretations. Here we use multiple regression methods to quantitatively examine the role of six abiotic parameters over the long-term variations of elasmobranch and actinopterygian genus-level diversity. We find that marine actinopterygian diversity is mainly controlled by temperature while continental fragmentation is the primary driver of the diversity fluctuations of elasmobranchs. Sea-level variations correlate positively with the diversity variations of both marine groups, whereas none of the tested proxies explains the diversity variation of freshwater ray-finned fishes. Our results indicate that such contrasting responses are mainly due to ecological and life-history trait differences between these groups.


INTRODUCTION
Understanding the relationships between biodiversity and global environmental variations through geological time represents a major focus in macroevolution and global change studies. However, this challenge has long been hampered by controversies over how much reliable the fossil record is (Erwin, 2009). It has been demonstrated that taking the fossil record at face value to estimate diversity variations produces inaccurate palaeodiversity reconstructions due to several geological and sampling biases that affect its direct reading (Raup, 1976;Peters and Foote, 2001;Smith, 2001;Peters, 2005Peters, , 2006Benton and Emerson, 2007). Over the past decade however, a number of analytical approaches have been developed to circumvent these biases and produce more reliable estimates of deep-time diversity variations based on fossil data (Norell, 1993;Alroy et al., 2008;Alroy, 2014;Silvestro et al., 2014). Meanwhile, geochemical studies have gathered a tremendous amount of data (Prokoph et al., 2008;Veizer and Prokoph, 2015) on a number of proxies that reflect climate and environmental parameters (Cárdenas and Harries, 2010;Hannisdal and Peters, 2011;Mayhew et al., 2012). These recent advances have permitted to test the hypothesis of a long-term control of various environmental factors over the fluctuations of taxonomic richness throughout the Phanerozoic. These studies mainly focused on marine invertebrate data (Mayhew et al., 2008(Mayhew et al., , 2012Cárdenas and Harries, 2010;Ezard et al., 2011;Hannisdal and Peters, 2011;Liow et al., 2015), whereas those that explored the role of some environmental variations over marine or continental vertebrate clades are comparatively scarce and concern low-diversity tetrapod groups (Marx and Uhen, 2010;Benson and Butler, 2011;Butler et al., 2011;Martin et al., 2014;Mannion et al., 2015). Multiple lines of evidence indicate that temperature acted as a major control over the evolutionary history of clades (Aguirre and Riding, 2005;Hunt et al., 2005;Marx and Uhen, 2010;Ezard et al., 2011;Figueirido et al., 2012;Mayhew et al., 2012;Condamine et al., 2013Condamine et al., , 2019Martin et al., 2014;De Blasio et al., 2015;Mannion et al., 2015), but it has been demonstrated that other factors such as sea-level (Peters, 2005;Benson and Butler, 2011;Hannisdal and Peters, 2011), tectonics (Zaffos et al., 2017), and nutrient availability (Cárdenas and Harries, 2010) likely acted over the evolution of biodiversity.
We focus on elasmobranchs (sharks, skates, and rays) and actinopterygians (ray-finned fishes)-the most speciose aquatic vertebrate clades in modern marine and continental ecosystems-to test the hypothesis of a long-term control of abiotic factors over vertebrates' diversity. These clades encompass a tremendous range of ecological, biological, morphological and environmental adaptations and their evolutionary history includes three of the biggest diversification events among jawed vertebrates (Alfaro et al., 2009). Previous family-level analyses indicated similar diversification patterns for marine actinopterygians and elasmobranchs, whereas the diversity variations of freshwater actinopterygians show contrasting patterns (Guinot and Cavin, 2015). However, quantitative analyses of the determinants of the diversification patterns of these clades are lacking at lower taxonomic levels. Here we use multiple regression methods to examine the role of six abiotic parameters over the long-term diversity variations of elasmobranch, marine and freshwater ray-finned fish clades based on bias-corrected diversity data at genus level.

Taxic Diversity Data
The diversity datasets used in the present work (Figures 1A,B, Supplementary Data 1) encompass all extant and extinct elasmobranch genera (541 taxa) over most of the time interval of the group (Upper Triassic to Recent) and all extinct actinopterygian genera (780) over the Upper Jurassic-Paleocene interval. Post-Paleocene phylogenetic and fossil record data for actinopterygians could not be included in the analyses for a two-fold reason: integrating extinct taxa into the still fluctuating molecular phylogenies of recent acanthomorphs alongside the necessary re-assessment of the systematic affinities of much of the enormous Eocene to Recent fossil record are tasks beyond the scope of this study.
Data were extracted from a phylogeny-based analysis of elasmobranch and actinopterygian fossil records and diversification patterns (Guinot and Cavin, 2016). This analysis produced a genus-level supertree for each clade and provided phylogenetic diversity estimates (PDE). For each timescaled phylogeny, polytomies were treated in two ways following Boyd et al. (2011): (i) resolved in chronological order [pectinate arrangement placing the taxon having the oldest First Appearance Datum (FAD) at the base of the resolved clade and the one with the youngest FAD in the most derived position] and (ii) resolved in reverse-chronological order (pectinate arrangement placing the taxon having the youngest FAD at the base of the resolved clade and the one with the oldest FAD in the most derived position). This produced two topologies (chronological and reverse-chronological polytomy resolutions) for each of the supertrees considered. Uncertainty in the age of FAD were considered in each tree topologies by randomly picking an absolute age within the FAD range of each taxon, for each 1,000,000 replicates as in Boyd et al. (2011). Two tree topologies were selected for each polytomy resolution method based on their congruence scores (Wills, 1999;Pol and Norell, 2001), which assess the fit of phylogenetic relationships with stratigraphic ranges of terminal taxa. This resulted in four phylogenetic diversity estimates (PDE) for each clade: most congruent chronological, least congruent chronological, most congruent reverse-chronological and least congruent reverse-chronological (see Guinot and Cavin, 2016). Actinopterygian PDE data were further divided into three datasets corresponding to three ecological preferences: (i) exclusively marine actinopterygian genera, (ii) exclusively freshwater actinopterygian genera, and (iii) genera from mixed environments (i.e., genera either including species from freshwater and marine environments, or including euryhaline taxa). The latter were excluded from our analyses because of equivocal environmental distributions. For each dataset, analyses were performed on each four PDE data. We also selected the median diversity value for each time bin in order to sum up all PDE in one curve (Figures 1A,B).

Selection of Explanatory Variables
Data on abiotic variables (Figures 1C-H, Supplementary Data 1) were taken from the literature and averaged for each time bin to comply with the diversity data. Two sets of data on sea level fluctuations were selected. Data from Haq et al. (1987) are based on sequence stratigraphy and span the Triassic-Recent interval. Sea level data from Miller et al. (2005) reflect eustatic changes due to variations in continental ice sheet volumes, which are recorded by δ 18 O from foraminifera over the Early Jurassic-Recent interval. Data of Miller et al. (2005), although more recent, do not span the entire time interval of the elasmobranch evolutionary history. We therefore created a composite sea level dataset by standardizing the original data at minimal and maximal values of 0 and 100, respectively, following Cárdenas and Harries (2010) and using the function scale () of the R-package "Scales" (Wickham, 2017). We then selected the scaled data of Miller et al. (2005) for the Aalenian-Recent interval and those of Haq et al. (1987) for the Carnian-Toarcian interval. Oxygen isotope records (expressed as δ 18 O values) were selected as representatives of trends in global climate change. Although climate change is multifactorial, δ 18 O records are considered as an inverse proxy for global climate change that mainly reflects temperature fluctuations (also global ice volume during icehouse periods) (Zachos et al., 2008). The δ 18 O data used here were taken from Veizer and Prokoph (2015). We selected the surface-water (mixed surface layer) data for the temperate and tropical realms excluding deep-sea and highlatitude data to preserve a homogenous signal over time bins and to comply with the dominant geographic and environmental distributions of the taxa included in our dataset. The δ 13 C data of Prokoph et al. (2008) were selected to represent a proxy of productivity and preservation of organic matter (Cárdenas and Harries, 2010). Proxies representing nutrient input from continental weathering (Martin, 1996;Shields, 2007) ( 87 Sr/ 89 Sr) and from recycling of organic material in ocean sediments (Cárdenas and Harries, 2010) (δ 34 S) over the Triassic-Recent interval were extracted form Prokoph et al. (2008). Missing value for two time bins (Tithonian and Hettangian) in the δ 34 S data were extrapolated using the spline interpolation method with the function aspline () of the R-package "akima" (Akima and Gebhardt, 2016) to perform an univariate akima interpolation at 0.5 years intervals. These were then averaged for each of the time bins used here. The index of continental fragmentation produced by Zaffos et al. (2017) that ranges from 0 (maximum tectonic plate aggregation) to 1 (complete plate fragmentation) was selected as another abiotic variable to represent temporal changes in the geographic arrangement of continental crust.
Proxies representing geological sampling were selected in order to account for the possible effect of this parameter in biasing diversity estimates. The elasmobranch and actinopterygian fossil records have been sampled since the end of the 18th century with uneven information about historic collections, which precludes a direct measure of the number of fish-bearing collections. Consequently, sampling effort is represented here by counts of fossiliferous geological stratae and we also considered counts of fossil collections, which were tested alternatively in our models. Collection or strata ages that span several time bins due to uncertainty were considered present in their whole age range. Data representing this proxy were downloaded from the Paleobiology Database (http:// paleobiodb.org, accessed on July 2018) and divided into marine and continental subsets.

Generalized Least Squares
This approach is a multiple regression method that offers the possibility of comparing a single dependent variable (here the PDE) to multiple explanatory variables (proxies). Whereas ordinary least square (OLS) regressions assume no serial correlation, generalized least squares (GLS) do not assume that data series or values within data series are independent and include autoregressive models in the analyses to account for autocorrelation issues. The set of models tested here comprises all possible combinations among the six explanatory variables considered as well as a null (intercept) model, which resulted in 64 models. The marine-related explanatory variables represented by δ 13 C and δ 34 S were not considered in analyses of the freshwater actinopterygian PDE because of their irrelevance as proxies for continental macroevolutionary patterns. However, we chose to keep the 87 Sr/ 89 Sr variable as it also represents variations in weathering intensity of landmasses and orogenic uplifts (Raymo and Ruddiman, 1992;Shields, 2007). This resulted in 33 models (including the null model) for the freshwater actinopterygian analyses. All models were duplicated to include a variable representing the rock record sampling, either expressed as the global count of fossiliferous stratae or as the number of fossil collections. This was made with the aim of considering the potential increased fit of adding such variables to our models, which would suggest incomplete sampling bias correction of the diversity variables. Following the method used in some previous macroevolutionary studies (e.g., Hunt et al., 2005;Marx and Uhen, 2010), we first fitted our 128 models (66 models for the freshwater data) using OLS and then fitted autoregressive models [function ar () of the R-package "Stats"] of order 1 (first-order autoregressive covariance model) to the residuals of the OLS. We then ran our models using GLS with the correlation structure of the residuals provided for the first order of autoregressive models. We subsequently compared the fit of autoregressive orders 0 (no autocorrelation) and 1 for each model and selected the best model for each combination of variables using the modified version of Akaike's information criterion for small sample sizes (AICc). The conditional probabilities for each model (Akaike weights, wAICc) as well as the generalized coefficient of determination (R 2 ) were computed to inform on the proportion of variance in the dependent variable that is explained by each model.
Non-normality of the residuals and heteroskedasticity (unequal variance of residuals) can cause issues in biasing estimates of the standard errors and of the variance of the coefficient estimators (incorrect testing of confidence intervals and hypothesis testing). We tested for normality and homoskedasticity of residuals using the Jarque and Bera (1980) and Breusch and Pagan (1979) tests, respectively. Most models for the elasmobranch dataset and some models in the marine and freshwater actinopterygian datasets contained non-normally distributed residuals that could not be removed by transforming the data (either Log or Square root transformation). For elasmobranch data, non-normality was avoided by deleting the data point for the Pleistocene-Recent interval and, in most instances, models that contained non-normality in all datasets were not ranked among the best fit models and did not affect the interpretation of our results. One exception is the freshwater actinopterygian analysis on the least congruent chronological PDE data, which contains non-normally distributed residuals. However, results were comparable to the same models for other freshwater actinopterygian PDE datasets with normally distributed residuals (see below). Heteroskedasticity was absent from the models of the actinopterygian datasets but present in some models of the elasmobranch dataset and we used heteroskedasticity-consistent standard error (and p-value) estimators (White, 1980) to correct estimates. We used the function vcovHC () of the R-package "Sandwich" to produce heteroskedasticity-consistent estimations of the variancecovariance matrix of coefficient estimates using the H3 type, as recommended by Long and Ervin (2000). However, such method might overcorrect for heteroskedasticity since significance is sometimes lost after H3 type correction for variables in models that contain heteroskedasticity, whereas similar homoskedastic models used for other PDE estimates provided significant correlations (see below).

Ecology Data
We compiled the main ecologies of living marine actinopterygian and elasmobranch species (Supplementary Data 2) with the aim of estimating the proportion of species living on (benthic) or close to (demersal, bathydemersal) the bottom among these two clades. Species with a range of ecologies that are not strictly benthic/demersal (e.g., benthopelagic) were not counted. Data for marine actinopterygians were extracted from FishBase (Froese and Pauly, 2019) through the R package "rfishbase" (Boettiger et al., 2012) and those for elasmobranchs were mainly taken from Ebert et al. (2013) and Last et al. (2016) with some complements from FishBase. For elasmobranchs, among a total of 1,203 species, 1,157 are marine and 1,003 marine species are benthic, demersal, or bathydemersal. Among the 14,324 strictly marine actinopterygian species, 6,560 species are benthic, demersal or bathydemersal.

RESULTS
Results of our analyses on the median diversity data are summarized in Table 1 and those of all PDE datasets are detailed in Supplementary Data 3-5. Regardless of the taxonomic dataset considered, none of the models including the sampling bias proxies as an additional explanatory variable is ranked among the best fit models (Supplementary Data 3-6). This tends to suggest that PDE provide robust estimations of diversity and efficiently correct for sampling biases. Alternatively, the sampling biases affecting the elasmobranch and actinopterygian datasets might not be accurately represented by the available global sampling bias proxies. However, in the absence of a direct measure of the number of "fish"-bearing collections (see section Selection of Explanatory Variables), there are no other available data that quantify sampling biases.
Results for the elasmobranch dataset indicate close AICc values for the first four best fit models as well as low wAICc (Supplementary Data 3), suggesting a nearly equivalent fit. However, all four best fit models, which altogether represent a cumulative wAICc of 0.807, indicate a strong (from 587 to 597 according to models) and statistically significant (p < 9.98 e−10 ) positive correlation coefficient for the continental fragmentation index. Sea level correlates weakly (0.909-0.92) but significantly in all first four models. Similarly, sulfur correlates moderately (3.14-3.28) in all four models with statistical significance for all of them, although significance is lost after the use of heteroskedasticity-consistent estimates for the two models that contain heteroskedasticity.
Results for the marine actinopterygian dataset indicate that the best fit model (wAICc = 0.497, R 2 = 0.986) is the one that comprises all six abiotic variables. Among them, surfacewater (mixed surface layer) δ 18 O values have a strong negative correlation coefficient (−37.57) with strong statistical support (p = 0.015). In addition, a positive correlation (4.35) is found significant (p = 0.019) for the variable representing sea level fluctuations. Similar results are found for all four PDE datasets taken independently (Supplementary Data 4). Analyses on the freshwater actinopterygian data identify the model including δ 18 O, continental fragmentation and 87 Sr/ 89 Sr data as explanatory variables (wAICc = 0.565, R 2 = 0.90). However, none of the correlation coefficients is statistically supported, and similar results are obtained for each PDE hypothesis (Supplementary Data 5).
Our results indicate that sea-level fluctuations correlate positively with the diversity variations of both marine actinopterygians and elasmobranchs. Furthermore, elasmobranch diversity is strongly positively correlated with continental fragmentation whereas diversity fluctuations of marine actinopterygians are predominantly affected by sea surface temperature variations. It could be expected that differences between both marine clades be due to the different time interval covered by the elasmobranch (Upper Triassic-Recent) and actinopterygian (Upper Jurassic-Paleocene) datasets. However, analyses ran with the elasmobranch dataset over the Upper Jurassic-Paleocene interval still indicate that continental Significance at alpha = 0.05 is indicated with asterisks (*). R 2 is the generalized coefficient of determination. AICc is the modified version of Akaike's information criterion for small sample sizes. Akaike weights (wAICc) represent the probability that the model is the best one given the observed data and considering the set of candidate models. AR order indicates the order of the autoregressive model. JB and BP indicate the p-value of the Jarque-Bera (non-normality) and Breusch-Pagan (heteroscedasticity) tests, respectively. All models were duplicated with the addition of a sampling proxy to account for the potentially biasing effect of sampling on PDE trends. Results shown here comprise the number of geological stratae as sampling proxy, those using the number of fossil collections are provided in Supplementary Data 6. Complete model rankings and parameter values for each PDE are included in Supplementary Data 3-5. Cont. Frag., continental fragmentation index.
fragmentation has a strong and significant positive correlation coefficient in all best fit models (Supplementary Data 7).

DISCUSSION
A positive relationship between diversity and continental flooding, through eustatic sea level change, has been observed for both marine invertebrates (Peters, 2005;Hannisdal and Peters, 2011) and marine tetrapods Butler et al., 2011;Mannion et al., 2015). Continental flooding has been considered as a potential driver of marine diversity through a species-area relationship (Peters, 2005) by controlling the extent of epicontinental seas [although the relationship might be more complex, see Holland (2012)], which host most of the marine vertebrate and invertebrate diversity today (Tittensor et al., 2010).
A previous macroevolutionary analysis on marine vertebrates indicated that sea level is positively correlated with the diversity of shallow marine tetrapods, but showed a weaker relationship with open sea taxa . It is likely that the positive correlation found with marine actinopterygians and elasmobranchs reflects the dominance of shallow water taxa in the clades analyzed here. Numerous studies have suggested or found a link between sea level (or continental flooding) and sampling bias metrics (Sepkoski, 1976;Peters and Foote, 2001;Smith, 2001;Peters, 2005Peters, , 2006Benton and Emerson, 2007), which led to establishing the "common cause" hypothesis (Peters, 2005). However, our results indicate a lower fit for models that include a sampling proxy, suggesting that sea level has a genuine biological effect and represents an important driver of the longterm fluctuations of marine biodiversity. Although sea level likely played a role in shaping the evolutionary history of marine fish clades, our results indicate that other abiotic variables had stronger effects on their paleobiodiversity variations. We found that marine actinopterygian diversity correlates strongly negatively with δ 18 O values, which are considered as an inverse proxy for global climate change that mainly reflects temperature fluctuations (also global ice volume during icehouse periods) (Zachos et al., 2008). The positive link between temperature and diversity has probably been the most widely recovered relationship in macroevolutionary studies on both marine and terrestrial groups and across many lineages (Aguirre and Riding, 2005;Hunt et al., 2005;Marx and Uhen, 2010;Ezard et al., 2011;Figueirido et al., 2012;Mayhew et al., 2012;Condamine et al., 2013Condamine et al., , 2019Martin et al., 2014;De Blasio et al., 2015;Mannion et al., 2015). Temperature is also recognized as acting on the spatial distribution of biodiversity in modern ecosystems through the latitudinal diversity gradient (Willig et al., 2003;Hillebrand, 2004;Tittensor et al., 2010). Although considering a direct relationship between climate and available energy might be an oversimplified view (Clarke and Gaston, 2006;Erwin, 2009), empirical data indicate that increasing energetic resources provided by warmer climates have increased mutation rates and shortened generation times (Rohde, 1992;Gillooly et al., 2002) and hence supported faster speciation and/or lower extinction rates (Allen et al., 2006). The correlation recovered here provides additional evidence for a positive effect of temperatures on clades' diversification and quantitatively confirms previous hypotheses that marine ray-finned fish diversity variations were linked with sea surface temperatures . This is exemplified by the tremendous Cenomanian and Palaeocene-Eocene actinopterygian diversifications (Guinot and Cavin, 2016), which occurred during hyperthermal events (Zachos et al., 2008;O'Brien et al., 2017). In addition, it has been demonstrated that higher seawater temperatures tend to shorten both planktonic larval and egg duration times in marine organisms (Hirst and Lopez-Urrutia, 2006;Duarte, 2007;O'Connor et al., 2007). Shorter planktonic larval and egg phases imply decreasing dispersal capabilities, which promotes isolation and speciation (Duarte, 2007;O'Connor et al., 2007). The vast majority of ray-finned fish groups include larval phase and most comprise a free/pelagic egg phase in their development (Kendall et al., 1984). This temperature-sensitive trait may partly explain the stronger effect of climate on actinopterygian diversity than on elasmobranchs, which produce less egg cases (or youngs per liters) and have no larval or free egg phases. However, a recent species-level analysis provided evidence that lamniform extinction rates were negatively correlated with temperature (Condamine et al., 2019), which indicates that this parameter also played a role in some clades of mostly pelagic elasmobranchs.
Nonetheless, our analyses show that elasmobranch diversity variations most strongly correlate with continental fragmentation, suggesting a first-order effect of this environmental parameter on elasmobranch diversification patterns. Plate tectonics has long been hypothesized to have impacted global biodiversity (Valentine and Moores, 1970) by reducing biodiversity during periods of continental coalescence and increasing biodiversity during breakup of continental masses. Under this hypothesis, continental fragmentation produces geographic barriers that lead to the evolution of new lineages via vicariance, which contributes to the global biodiversity increase. In the marine realm, continental fragmentation also produces more coastlines and neritic environments that are known to correlate with the taxic diversity of coastal groups in modern ecosystems (Tittensor et al., 2010). However, this plate tectonic regulation hypothesis lacked quantitative testing in a macroevolutionary context until a recent analysis demonstrated that the long-term trends in biodiversity of Phanerozoic skeletonized marine invertebrates are well predicted by the supercontinent coalescence-breakup cycle (Zaffos et al., 2017). The temporal similarities between the fluctuations of elasmobranch biodiversity and the continental coalescence-breakup cycles (Figures 1A,D) is exemplified by the almost continuous Jurassic-Late Cretaceous diversity increase that matches a steady rise in continental fragmentation, which reflects the breakup of Pangaea and particularly of Gondwana (Zaffos et al., 2017). This result on a marine vertebrate clade supports previous invertebrate-based evidence that the shifting distribution of continental landmasses has been a major driver of long-term global Phanerozoic biodiversity patterns Moores, 1970, 1972;Zaffos et al., 2017). However, we found no significant effect of continental fragmentation on the biodiversity variations of marine actinopterygians, which, added to our results for temperature, indicates contrasting responses to long-term environmental changes between these two clades. Living elasmobranch species are dominated by benthic and demersal taxa (86.7%), while these ecologies represent minor components (45.8%) of marine actinopterygian diversity (Supplementary Data 2). Fossil-based estimates of the ecological distribution of Mesozoic ray-finned fishes remain to be completed. However, the large post-Mesozoic diversification of reef-associated demersal fish lineages (Alfaro et al., 2007;Cowman and Bellwood, 2011) and the absence during the Mesozoic of typical demersal clades such as the flatfishes (Pleuronectiformes) (Friedman, 2008), or of extinct clades with similar ecology, indicate that the demersal component of rayfinned fishes was even lower during the Mesozoic than during the Cenozoic. These results are consistent with the expectation that habitat area and geographic barriers more strongly influence species richness of clades dominated by benthic taxa, which are more prone to vicariance due to limited dispersal capabilities (Hart and Pearson, 2011). Finally, we note that δ 34 S correlates positively with elasmobranch diversity, but statistical significance is not present in all datasets. Although regarded as a proxy for organic nutrient inputs or shelf redox conditions (Cárdenas and Harries, 2010;Mayhew et al., 2012), the ecological significance of δ 34 S is contentious (Cárdenas and Harries, 2010;Hannisdal and Peters, 2011) and a full understanding of the direct link between this parameter and biodiversity variations requires further study. Apart from sea level, which appears to be a common driver in the marine groups investigated here, our results indicate that different abiotic parameters played a role between the macroevolutionary history of elasmobranchs (continental fragmentation) and marine actinopterygians (climate change).
An analysis of the fish family-level diversity dynamics demonstrated that both marine actinopterygians and elasmobranchs follow a similar equilibrium model of diversity variation (Guinot and Cavin, 2015). This apparent contradiction with the present analyses may relate to the use of different taxonomic levels as it has been proposed that logistic diversity curves prevail for higher taxonomic ranks and gradually change toward an exponential distribution when lower levels are considered (Benton, 1997;Lane and Benton, 2003;Benton and Emerson, 2007). Although post-Paleocene genus-level diversity data are lacking for actinopterygians, our pre-Eocene data ( Figure 1B) suggest that such contrasting taxonomic-level patterns might occur in marine ray-finned fishes. Yet, the genus-level elasmobranch diversification pattern also fits an equilibrium model (Guinot and Cavin, 2015), and our results indicate that this equilibrium model is mainly tectonically driven. This suggests that continental fragmentation probably drove higher-level taxonomic richness, which can be regarded as reflecting key morphological or physiological adaptations (Guinot and Cavin, 2015), and that diversification patterns at lower taxonomic levels are influenced by different parameters according to the ecology of the clades.
Another noticeable result is the absence of correlation between the tested variables (continental fragmentation, sea level, δ 18 O and 87 Sr/ 89 Sr) and the long-term freshwater ray-finned fish diversity fluctuations. This result should be parallelized with previous family-level analysis showing that the diversity dynamics of freshwater ray-finned fishes do not follow that of their marine relatives, i.e., the variation is better described with an exponential function rather than with a logistic one (Guinot and Cavin, 2015). Potential causes at the origin of this pattern are proposed and discussed by Guinot and Cavin (2015). The present study highlights the fact that continental fish diversity may be under control of abiotic parameters that are unquantifiable at the macro scale (e.g., modifications of regional hydrographic systems, orogenesis). In addition, biotic factors such as competition and predation, which are known to enhance speciation and to be more intense in the freshwater realm (Vermeij and Grosberg, 2010), may play a preponderant role in shaping continental fish diversity. It has been proposed that sea level/continental flooding may have exerted a control on the diversity variations of continental clades on the basis of the well-established species-area relationship (Barnosky et al., 2005), but none of the analyses that quantitatively tested this hypothesis provided significant results Mannion et al., 2011). Our results on continental actinopterygians lend support to the absence of a long-term effect of sea level on diversity in the continental realm, although some large variations probably punctually affected diversity changes (Hallam, 1989;Hallam and Wignall, 1999).
Our analyses cover all available abiotic parameters that are so far quantifiable and proposes potential environmental drivers and their outcomes on the two most speciose vertebrate clades in today's aquatic realms. The contrasting response of different "fish" clades to extrinsic parameters indicates that biological factors such as life-history traits and ecology likely control the response of clades to environmental change. Many macroevolutionary studies have analyzed the drivers of the Phanerozoic marine invertebrate clades as a whole (Cárdenas and Harries, 2010;Mayhew et al., 2012), with no distinction between groups (taxonomic, biological, physiological, ecological). Based on our results, it is expected that results found in such Phanerozoic marine invertebrate analyses are unbalanced by dominant ecologies among marine invertebrates and that the real patterns of abiotic control on deep-time biodiversity variations are more complex than previously expected. Although this is a first step in the understanding of what shapes the evolutionary history of clades, biotic interactions including competition Silvestro et al., 2015;Voje et al., 2015) and predation (Huntley and Kowalewski, 2007;Marx and Uhen, 2010) have been more scarcely tested but have been demonstrated to play a role in the diversification patterns of clades, including sharks (Condamine et al., 2019). Exploring these intrinsic biotic effects represent another challenge for future studies, as is the understanding of the long-term drivers of continental fish diversity.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
GG and LC designed the study and wrote the manuscript. GG compiled the data and performed the analyses.

FUNDING
This paper is a contribution to the project Fish response to longterm global change supported by the Swiss National Science Foundation .

ACKNOWLEDGMENTS
We wish to thank Emmanuel Paradis and Roger Benson for helpful discussions on the methodology. Fabien Condamine is thanked for his support and fruitful discussions as well as the two reviewers for providing useful comments on an earlier version of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.

2019.00513/full#supplementary-material
Supplementary Data 1 | Abiotic and diversity data used in the analyses.
Supplementary Data 2 | Ecological data for living marine actinopterygian and elasmobranch species.
Supplementary Data 3 | Results of generalized least square analyses on elasmobranch data.
Supplementary Data 4 | Results of generalized least square analyses on marine actinopterygian data.
Supplementary Data 5 | Results of generalized least square analyses on freshwater actinopterygian data.
Supplementary Data 6 | Results of generalized least square analyses on PDE data of elasmobranch, marine, and freshwater data using the number of fossil collections as proxy for sampling.
Supplementary Data 7 | Results of generalized least square analyses on elasmobranch data over the Oxfordian-Thanetian interval.