Phylogenetic measures reveal eco-evolutionary drivers of biodiversity along a depth gradient Research

Energy and environmental stability are positively correlated with species richness along broad-scale spatial gradients in terrestrial ecosystems, so their relative importance in generating and preserving diversity cannot be readily disentangled. This study seeks to exploit the negative correlation between energy and stability along the oceanic depth gradient to better understand their relative contribution in shaping broadscale biodiversity patterns. We develop a conceptual framework by simulating speciation and extinction along energy and stability gradients to generate expected patterns of biodiversity for a suite of complementary phylogenetic diversity metrics. Using a time-cali-brated molecular phylogeny for New Zealand marine ray-finned fishes and a replicated community ecological sampling design, we then modelled these metrics along large-scale depth and latitude gradients. Our results indicate that energy-rich shallow waters may be an engine of diversity for percomorphs, but also suggest that recent speciation occurs in ancient fish lineages in the deep sea, hence questioning the role of energy as a key driver of speciation. Despite potentially facing high extinction early in their evolution, ancient phylogenetic lineages specialized for the deep-sea were likely preserved by environmental stability during the Cenozoic. Furthermore, intermediate depths might be a ‘museum’ (or zone of overlap) for distinct lineages that occur predominantly in either shallow or deep-sea waters. These intermediate depths (500–900 m) may form a ‘phylogenetic diversity bank’, perhaps providing a refuge during ancient (Mesozoic) anoxic events affecting the deep sea and more recent (Pliocene–Pleistocene) climatic events occurring in shallow ecosystems. Finally, the phylogenetic structures observed in fish communities at intermediate depths suggest other processes might restrict the co-occurrence of closely related species. Overall, by combining a conceptual framework with models of empirical phylogenetic diversity patterns, our study paves the way for understanding the determinants of biodiversity across the largest habitat on earth.


Introduction
Despite harbouring some of the strongest environmental gradients on earth, little is known regarding mechanisms shaping biodiversity along the ocean's depth gradient (Ramirez-Llodra et al. 2010). Along other broadscale gradients (e.g. latitude, elevation), biodiversity tends to increase with increasing energy, temperature and environmental stability (e.g. at low latitude) (Mittelbach et al. 2007, Jetz and Fine 2012, Fine 2015. Energy, temperature and environmental stability are, however, positively correlated along latitudinal and elevational gradients (Fine 2015, McClain and Schlacher 2015, Valentine and Jablonski 2015. In contrast, while energy and temperature both decrease with increasing depth, environmental stability increases (Fig. 1a) (Ramirez-Llodra et al. 2010, Fine 2015, Valentine and Jablonski 2015. Since the global anoxic events (OAEs) that caused mass extinction during the Mesozoic (Priede and Froese 2013) as recently as 94-93 Ma ago (Jenkyns 2010), environmental conditions in the deep sea have remained reasonably stable (Priede and Froese 2013). Owing to the buffering effect of the water column (seasonal variation rarely occurs deeper than 400 m [Fiedler 2010]), deeper ocean environments would have been less affected by the great seasonal variation of the Eocene/ Oligocene boundary (Ivany et al. 2000) and climatic oscillations during the Miocene and Pleistocene (e.g. glacial/ inter-glacial periods) than shallow-water environments (Pillans et al. 1998, McClymont et al. 2016. Thus, studying patterns of biodiversity along the depth gradient offers an opportunity to disentangle the often-confounded drivers of Figure 1. (a) Correlations between energy (black) and stability (gray) along latitudinal and elevation gradients and along the depth gradient. (b) Conceptual framework based on simulations (Supplementary material Appendix 1) proposed to disentangle the role of energy versus stability on evolutionary processes shaping biodiversity over depth. The central tree represents a hypothetical phylogenetic tree of a community (i), and each quadrant shows what phylogenetic structures and associated changes in phylogenetic diversity values that would be expected for communities at different positions along the energy and stability axes, given their expected effects on speciation and extinction rates. SR = species richness, PD = phylogenetic diversity, MPD = mean pairwise distance, VPD = variance of the pairwise distance, MNTD = mean nearest taxonomic distance. Right: the predicted change in phylogenetic diversity measures (PD, MPD, VPD, MNTD) versus depth when (ci, cii) energy alone, (d) stability alone, or (ei, eii) both stability and energy are driving patterns in biodiversity. energy/temperature and stability on biodiversity (Fine 2015, Valentine andJablonski 2015).
Examining the phylogenetic structures of communities, rather than just taxonomic species richness, has enhanced our understanding of broadscale biodiversity patterns and their drivers (Cavender-Bares et al. 2009, Mouquet et al. 2012. Several alpha (local) diversity measures have been proposed to capture different aspects of this phylogenetic structure (Vamosi et al. 2009, Pavoine and Bonsall 2011, Tucker et al. 2017. For example, phylogenetic diversity (PD) captures the overall amount of evolutionary history in a sampled community as a sum of the branch lengths linking species through a phylogeny (Faith 1992), whereas mean pairwise distance (MPD, Webb 2000) represents the average phylogenetic distance (or average time divergence) among species, and is influenced by ancient diversification events (Mazel et al. 2016). In contrast, mean nearest taxonomic distance (MNTD, Webb 2000) estimates the average phylogenetic distance among the closest relatives within a community, and is therefore influenced by recent speciation events (Mazel et al. 2016). Last, whereas the variance of the pairwise distance (VPD, Clarke and Warwick 2001) measures variability in phylogenetic distances among species, influenced by both recent and ancient divergences, the variance of the nearest taxonomic distance (VNTD, Tucker et al. 2017) focuses on closest relatives, so is not influenced by ancient divergence. As each measure draws on unique information from the phylogenetic relationships among species, documenting patterns for several of these measures simultaneously can provide novel insights into the mechanisms shaping biodiversity along broad-scale gradients, such as depth.
Here we present a simple conceptual framework (Fig. 1b) that predicts patterns of phylogenetic alpha diversity (hereafter phylo-α) under different scenarios of speciation and extinction along broad-scale geographical/environmental gradients (see Supplementary material Appendix 1 for simulation details). Our framework draws on macroecological hypotheses commonly proposed for both terrestrial and marine systems that suggest a positive relationship between energy/temperature and speciation, and a negative relationship between long-term environmental/geological stability and extinction (Rohde 1992, Allen et al. 2006, Davies et al. 2007, Mittelbach et al. 2007, Gillman et al. 2009, Condamine et al. 2013, Rolland et al. 2014, Siqueira et al. 2016. In our framework we explore a broad spectrum of extinction rates/scenarios for a given speciation level and vice-versa (Fig. 1b, Supplementary material Appendix 1), and acknowledge the potential for a negative relationship between energy and extinction (Evans et al. 2005), and that long-term environmental stability may also increase speciation rate through specialization (Jansson and Dynesius 2002).
We began with a hypothetical community phylogenetic tree generated with intermediate levels of speciation and extinction (Fig. 1bi) and portray how this structure, and associated phylogenetic diversity measures, would change with increasing or decreasing speciation or extinction rates. For example, high speciation and low extinction will cause trees to have higher PD and MPD values, but shorter distances among close relatives (low MNTD) and lower VPD (Fig. 1bii). This pattern typifies a 'cradle' (generating new species) as well as a 'museum' (preserving old lineages), likely to occur in regions of high energy and low environmental/ geological variability (e.g. tropical regions, Jablonski et al. 2006, Pyron and Wiens 2013, Rolland et al. 2014). An increase in the extinction probability (i.e. more frequent mass extinction events) in conjunction with a high speciation rate, as expected in energy-rich but unstable regions, would generate phylogenies with long basal branches creating distinct lineages with many short terminal branches. This phylogenetic structure would yield lower values for PD and MPD, increases in VPD, and MNTD would remain stable and low (Fig. 1biii). High extinction in association with low speciation, as expected in highly unstable and energy-poor regions, would generate depauperate phylogenies having a mixture of long and short branches, decreasing the PD and MPD, but increasing the VPD and the MNTD (Fig. 1biv). Finally, under a scenario of low speciation and extinction, as expected in stable but low-energy areas, species-poor phylogenies with an even branch-length distribution would be generated, yielding communities with relatively low PD and VPD, but higher MPD and MNTD (Fig. 1bv).
Considering this conceptual framework in the context of the depth gradient, we made specific predictions regarding expected changes in phylo-α according to underlying drivers. First, if energy/temperature boost speciation and decrease extinction (Allen et al. 2006, Pyron andWiens 2013), PD and MPD should be highest in energy-rich surface waters and decrease with depth, while VPD and MNTD should be lowest in shallow waters, then increase with depth ( Fig. 1bii, iv, ci). Alternatively, if energy is not the major determinant of speciation rates in marine systems (Rabosky et al. 2018, O'Hara et al. 2019, but rather only decreases extinction rates (Evans et al. 2005), MNTD should be constant along the depth gradient, PD and MPD should be highest in shallow waters and decrease with depth, while VPD should be lowest in shallow waters, then increase with depth ( Fig. 1bii, biii, cii). Second, if long-term stability decreases extinction, and increases the opportunity for speciation through specialization (Jansson and Dynesius 2002, Kozak and Wiens 2012, Salisbury et al. 2012, Rolland and Salamin 2016, Gajdzik et al. 2019), we expect PD and MPD to increase with depth toward more stable environments, while VPD and MNTD should decrease with depth ( Fig. 1bii, iv, d). Third, if both energy/temperature and long-term stability yield low extinction and high speciation, PD and MPD would be U-shaped versus depth, while MNTD and VPD would be unimodal (Fig. 1ei). Finally, if energy/temperature does not affect speciation (Rabosky et al. 2018, O'Hara et al. 2019) but only extinction, and long-term stability yields low extinction and high speciation, PD should be consistently low across shallow and intermediate depths before increasing at deeper depths, and MPD would be U-shaped (extinction decreases in both shallow and deep waters), while MNTD and VPD would decrease with depth (Fig. 1b,bi,bii,bv,eii).
Previous studies of biodiversity versus depth describe either a decline in species richness with depth or a peak at intermediate depths (Ramirez-Llodra et al. 2010, Costello andChaudhary 2017). Patterns in taxonomic diversity indices (e.g. AvTD) versus depth are also variable (Tolimieri and Anderson 2010, Campbell et al. 2011, Zintzen et al. 2011. Studies that have inferred evolutionary processes along large-scale gradients in the ocean using near-complete species phylogenies suggest that speciation rates may be higher in marine fishes (Rabosky et al. 2018) and ophiuroids (O'Hara et al. 2019) found in polar regions (possibly owing to instability since the development of icecaps around 41-34 Ma [Crame 2018]) and at depth, challenging the idea that energy is the primary driver of speciation rates in the ocean.
Our study provides the first analysis of complementary phylo-α measures using a time-calibrated phylogeny to examine phylogenetic structures for communities of marine ray-finned fishes sampled along broad-scale ocean gradients. Specifically, we modeled variation in each of five phylo-α indices versus depth (50-1200 m; replicated across latitude) within New Zealand waters to test our predictions regarding the relative roles of speciation and extinction ( Fig. 1c-e), potentially driven by energy/temperature and long-term environmental stability. We used multivariate cluster analysis of our phylo-α measures to identify 'phyloregions' within which fish communities demonstrate similar phylo-α signatures, hence are likely to have been influenced by similar underlying drivers.

Data acquisition
Phylo-α measures of fish communities were calculated based on presence/absence data for 149 teleost species identified from baited remote underwater video (BRUV) apparatuses deployed in a structured replicated sampling design (Zintzen et al. 2017). Fish species were recorded from 2 h of video footage obtained by each BRUV deployment. The BRUVs were deployed at seven locations spanning 21° latitude within New Zealand waters: Rangitahua (the Kermadec Islands, KER), the Three Kings Islands (TKI), Great Barrier Island (GBI), White Island (Wis), Kaikōura (KKA), the Otago Peninsula (OTA) and the Auckland Islands (AUC). For each location, 5-6 replicate videos were taken at each of seven targeted depths: 50, 100, 300, 500, 700, 900 and 1200 m.

Phylogenetic tree building
We built a Bayesian posterior distribution of time-calibrated molecular phylogenies for New Zealand marine ray-finned fishes (class Actinopterygii) following the procedure of Eme et al. (2019). We used their same 15 gene supermatrix, adding DNA sequences for 14 species never sequenced before, and also added 6 taxa that lacked sequence data. Taxa lacking sequence information were added to the tree reconstruction using conservative topological constraints based on their taxonomic status. The overall supermatrix included 803 species (797 species with DNA and 6 outgroup species not present in New Zealand) and 3679 DNA sequences. The tree reconstruction included 165 soft topological constraints to avoid spurious phylogenetic relationships due to data deficiency (incomplete taxon/gene sampling). The tree was calibrated in absolute time using 44 fossil calibration points. Four independent MCMCs were constructed using BEASTv2.4.7 (Bouckaert et al. 2014) to obtain a posterior distribution of time-calibrated trees and a maximum clade credibility tree (MCCT). We randomly sampled 100 trees from the posterior distribution to quantify phylogenetic uncertainty in downstream analyses (Rangel et al. 2015) (Supplementary material Appendix 2 for laboratory and tree-building methods, and Data availability statement).

Phylogenetic diversity metrics
We calculated five phylo-α indices to capture different aspects of the phylogenetic structure of communities (Tucker et al. 2017). Although phylogenetic diversity (PD, Faith 1992) is strongly correlated with species richness (SR), it provides insights regarding rates of speciation. For instance, low PD associated with high SR indicates a high rate of speciation (Davies et al. 2007). Mean phylogenetic distance among species (MPD, Warwick 1998, Webb 2000) is independent of SR, but is sensitive to large phylogenetic distances contributed by co-occurring species belonging to old and distinct lineages (Mazel et al. 2016), whereas mean nearest taxonomic distance (MNTD) is influenced by recent speciation only (Webb 2000, Mazel et al. 2016. Variation in phylogenetic distances among species within a community (VPD) is independent of SR and MPD (Clarke and Warwick 2001). This regularity index captures more complex phylogenetic structure, such as the presence of distinct old lineages with recent clusters of closely related species (Zintzen et al. 2011), a signature of high extinction in the past (Sanmartín and Meseguer 2016), and/or strong environmental filtering (Cavender-Bares et al. 2009, Zintzen et al. 2011). Finally, the variability of phylogenetic distances among close relatives (VNTD, Tucker et al. 2017) was used to capture the temporal breadth of recent speciation. We also calculated correlations among these metrics for our dataset (Supplementary material Appendix 2).

Simulation of community phylogenetic structures
To build our conceptual framework (Fig. 1b), we calculated SR and the phylo-α measures from hypothetical trees generated under varying speciation and extinction rates using four sets of simulations. Simulation sets differed by relaxing the extinction rate through time (including mass extinction events) and sampling different fractions of extant species. Because MNTD and VNTD showed similar trends in our simulation framework, only MNTD is presented in Fig. 1, but both metrics were included in the empirical study. For details regarding simulation procedures and results see Supplementary material Appendix 1 Fig. A1

Statistical analysis
There were 47 depth-by-latitude cells in our full sampling design (there were no BRUV deployments at 1200 m for either White Island or the Auckland Islands, due to poor weather conditions, see Zintzen et al. 2017). From the full BRUV dataset, we pooled the species sampled from BRUVs within each cell, hereafter called the 'pooled 47-cell dataset'.
We used the MCCT tree to calculate phylo-α values for each metric in each cell and performed model selection with depth and latitude as explanatory variables using Akaike's information criterion with small-sample-size correction (AICc). Model selection performed on average values obtained over 100 trees did not alter the results (r = 0.98-0.99 for all metrics). For each metric, we considered four sets of models (Supplementary material Appendix 2 Table A2): 1) a two-factor ANOVA design with depth and latitude as factors; 2) an ANCOVA design with latitude as a factor and depth as a covariate; 3) an ANCOVA design with depth as a factor and latitude as a covariate; and 4) a regression model with depth and latitude as continuous predictors including quadratic terms for both and their interactions. With each set of models, we also considered simpler sub-set models and retained the model with the lowest AICc. Model selection was done in R (R Core Team) using ordinary least squares models. To quantify phylogenetic uncertainty, we fitted the best model for each phylo-α metric for all 100 trees and computed the mean and variance of parameter estimates and their contributions towards the explained variance. To estimate the proportion of variation explained by phylogenetic uncertainty, we also performed a linear mixed-effects model with trees as a random effect (Supplementary material Appendix 2).
We used clustering methods to identify bioregions of similar phylogenetic characteristics using all five phylo-α measures based on the MCCT. First, we used the 'find.cluster' function of the adegenet R package (Jombart 2008) to perform sequential k-means analyses for k = 2-15 clusters, and chose the optimal number of clusters using the recommended 'diffNgroups' option (i.e. the biggest difference in Bayesian information criterion values between successive values of k, Jombart et al. 2010). To check the general robustness of our clusters, we also performed hierarchical agglomerative groupaverage clustering (UPGMA). Next, we repeated k-means clustering for the optimal number of clusters for all 100 trees to accommodate phylogenetic uncertainty. Analyses were repeated excluding PD to decrease the influence of SR. Finally, to identify biodiversity patterns at a finer spatial scale, all analyses were also done using values obtained for each phylo-α metric at the individual BRUV level, then averaged for each of the 47 cells.

Results
The MCCT showed general agreement with the tree published by Eme et al. (2019) (Supplementary material Appendix 2 Fig. A9).

Modeling phylogenetic diversity versus depth and latitude
The phylo-α metrics varied more strongly with depth than with latitude (Fig. 2). All metrics, except VNTD, showed a significant relationship with depth. SR and PD (which were highly correlated, Supplementary material Appendix 2 Fig.  A10, A11) decreased sharply with depth between 50 and 300 m, then showed a secondary peak at 900 m. SR and PD had significant unimodal relationships with latitude, peaking around 37°S. MPD showed an unexpected unimodal relationship with depth, peaking at 700-900 m (Fig. 2c). MNTD also had a unimodal relationship with depth, showing a maximum between 300 m and 700 m and a minimum at 1200 m ( Fig. 2d), matching the expected relationships in Fig. 1ei. VPD increased sharply with depth, especially beyond 500 m, while VNTD remained fairly constant throughout the depth and latitudinal ranges examined here (Fig. 2f, l). These patterns echoed phylogenetic structure in species' occurrences versus depth, with a dominance of Percomorpha in shallow waters being replaced by species belonging to more ancient lineages (Anguilliformes, Notacanthiformes, Alepocephaliformes, Stomiiformes, Gadiformes) at depth (Supplementary material Appendix 2 Fig. A12). In contrast, there was no obvious phylogenetic structure in species' occurrences versus latitude (Supplementary material Appendix 2 Fig. A13).
Model selection always favoured models with depth and latitude as continuous variables. The parameter estimates, their significance and the explained variance of individual terms retained in the best models remained very consistent across the distribution of 100 trees for all phylo-α metrics (Fig. 3, Table 1), indicating robustness of observed patterns to phylogenetic uncertainty. The best model for all phylogenetic indices, except VNTD, included depth (Table 1) as either a linear or quadratic term (for MPD, MNTD and VPD), and its explained variance was always higher than for latitude, except for PD (Fig. 3  Latitude ( Figure 2. Relationships of diversity indices with depth (top) and latitude (bottom), when considering the average value estimated over 100 trees on the pooled 47 cells dataset for: species richness (a, g), phylogenetic diversity, PD (b, h); mean pairwise distance, MPD (c, i); mean nearest taxon distance, MNTD (d, j); variance of the pairwise distance, VPD (e, k); and variance of the nearest taxon distance, VNTD (f, l). Black horizontal bars, black dots and boxes show the median, outliers and interquartile range respectively. Whiskers are up to 1.5 times the interquartile range. Blue dots represent the average value per depth/latitude. For visualization purposes only, the blue curves represent the fit of the best model (dashed lines represent non-significant trends) selected based on AICc among three alternative models: model with a single term, a model with single and quadratic terms, and a model with single, quadratic and cubic terms for depth (a-f ) or latitude (g-l). The model selection was performed separately for depth and latitude. For species richness, a generalized linear model with a negative binomial error model with a log link was used, for all other metrics we used an ordinary least square models with a Gaussian error, and identity link.
For PD, the models explained on average 50% of the variance and the interaction between depth and latitude was the most important term in the model (17.9% of explained variance) followed by the quadratic effect of latitude (Table 1, Fig. 3). For MPD, the models explained 51% of the variance and depth was the most important (i.e. 44%, summing depth and depth 2 ). Latitude and its quadratic term remained significant but their contribution was minor (i.e. below 7% when summing latitudinal effects, Table 1). For MNTD, the model explained 38% of the variance and the quadratic term for depth contributed the most (23%), followed by the depthby-latitude interaction (10%). For VPD, depth was the most important predictor, explaining 70% of model variation (summing depth and depth 2 ). There were sharp increases in VPD along the depth gradient across all latitudes (Fig. 2e, Supplementary material Appendix 2 Fig. A14d). For VNTD, none of the model parameters were statistically significant.
The linear mixed-effects models (with trees as a random factor) yielded identical results regarding the most important model terms (Supplementary material Appendix 2 Table A3). Phylogenetic uncertainty (i.e. the random factor) contributed little to overall variation (0-2%, Supplementary material Appendix 2 Table A3).
Similar results were obtained when analyses were done based on smaller sampling units (individual BRUV units) for all metrics, except for VNTD (Supplementary material Appendix 2 Fig. A15-A17, Table A4, A5). The latter showed a significant positive relationship with depth (Supplementary material Appendix 2 Fig. A15), explaining 36% of the variance in the best model (when summing depth and depth 2 , Supplementary material Appendix 2 Table A4).

Phylogenetic bioregions across depth and latitude
K-means clustering of the phylo-α measures for the pooled 47-cell dataset identified four clusters that were robust to phylogenetic uncertainty ( Fig. 4a-b, Supplementary material Appendix 2 Fig. A18) (VNTD was removed because it was uncorrelated with latitude or depth). The first phylogenetic bioregion (hereafter 'phyloregion') corresponded to low-latitude shallow waters (in red, Fig. 4) characterized by high PD and low MPD, MNTD and VPD, suggesting a scenario of high speciation, and intermediate-to-high extinction by reference to the conceptual framework (Fig. 1, Supplementary material Appendix 1 Fig. A5). This suggested that energy might be a driver of recent and sustained speciation in shallow and subtropical waters (see Discussion). A second phyloregion corresponded to intermediate depths (300-700 m) at low latitudes plus shallow waters (50-500 m) at high latitudes, which had high MNTD and low MPD, PD and VPD (in orange, Fig. 4). A third bioregion (in blue) corresponded to depths around 700-900 m and some shallower depths (500-300 m) at higher latitudes; it was characterized by high MPD, MNTD and intermediate VPD and PD values. These two phyloregions (orange and blue) characterized by intermediate stability and energy (i.e. mesotrophic), suggests a  scenario of low speciation and low to intermediate extinction (see the Discussion regarding other potential processes that may explain this unexpected high MPD and low VPD). The fourth phyloregion included the deepest depth strata across all latitudes. It was characterized by low MNTD and PD, but with high MPD and VPD (in black, Fig. 4). This group aligns with a scenario of intermediate to high recent speciation and high past extinction (Fig. 1, Supplementary material Appendix 1 Fig. A5), suggesting that the lack of energy might not preclude recent speciation, and that deep water might not always have been stable. These phyloregions were robust to 1) the removal of PD (Supplementary material Appendix 2 Fig. A19) and 2) the inclusion of VNTD when considering smaller-scale sampling units (Supplementary material Appendix 2 Fig. A20, A21). Results obtained using hierarchical agglomerative clustering generally mirrored the results found using the k-means approach, although there were some minor discrepancies for cells occurring at the interfaces between groups (Supplementary material Appendix 2 Fig. A18-A21).

Discussion
We have revealed new broadscale phylogenetic diversity patterns for marine ray-finned fishes. Our study indicates that both energy and environmental stability are important drivers generating phylogenetic structures in fish communities along the depth gradient, where energy and stability are inversely related. Although there are some discrepancies between the observed and expected relationships (Fig. 1) for some phylogenetic diversity metrics, overall our results suggest that energy-rich shallow waters are an engine (i.e. 'cradle') of diversity for percomorphs, and that recent speciation for ancient lineages occurring in the deep sea is not driven by energy. Despite potentially facing high extinction early in their evolution, ancient phylogenetic lineages specialized for the deep sea were likely preserved during the Cenozoic by environmental stability. Furthermore, our findings of high MPD and MNTD at intermediate depths suggest that these depths might be a 'museum' (or zone of overlap) for the distinct phylogenetic lineages that occur either predominantly in shallow or deep-sea ecosystems. Contrary to previous findings (Brown and Thaje 2014), our study proposes that these intermediate depths might have served as a 'refuge' during extreme historical environmental changes such as the Mesozoic anoxic events, and the more frequent and recent climatic events in the shallow ocean scattered through the Cenozoic, rather than being a cradle of speciation. Finally, the phylogenetic community structures observed in these mesotrophic waters suggest some other processes might restrict the co-occurrence of closely related species.

New phylogenetic diversity patterns versus depth for vertebrates
For the first time, our study documents phylogenetic diversity patterns for the most speciose group of vertebrates (Actinopterygii) across broad-scale depth and latitude gradients using a molecular phylogeny. In this study, phylogenetic diversity (PD) matched species richness, as expected (Davies andBuckley 2011, Fritz andRahbek 2012). Both generally decreased with depth (this decrease was more pronounced at low latitudes), but peaked at intermediate depths for some high-latitude locations, consistent with previous work (Leathwick et al. 2006, Ramirez-Llodra et al. 2010, Brown and Thatje 2014, Costello and Chaudhary 2017. Such an interaction between depth and latitude is likely the reason a cubic relationship (including a secondary peak at intermediate depths) was found to be the best model when depth alone was considered as an explanatory variable (Fig. 2a). The increase of VPD with depth also confirmed previous findings for fishes, supporting the co-occurrence of highly divergent fish lineages likely specialized for the deep sea that have subsequently undergone recent speciation, creating clusters of related species (Zintzen et al. 2011, Rabosky et al. 2018. In contrast, MPD and MNTD presented new and unexpected patterns with depth for vertebrates; MPD was maximal at 900 m while MNTD showed a broader unimodal pattern suggesting low levels of close relatives at intermediate depths (300-700 m), a pattern also observed for brittlestars (O'Hara et al. 2019). Finally, an increase in VNTD with depth observed only at the smaller spatial scale may be attributed to consistent co-occurrence of certain lineages (Anguilliformes, Notacanthiformes, Aulopiformes), but high species turnover of close relatives within these lineages. Such fine spatial partitioning among close relatives may be driven by small-scale niche differentiation or biotic interactions such as competition (Cavender-Bares et al. 2009). Clustering of sampled cells according to phylo-α measures revealed four main regions having similar phylogenetic structures. Phyloregions were primarily stratified by depth and secondarily by latitude (Fig. 4a). Shallow depth strata (50, 100 m) had the highest turnover of phyloregions across latitude (red, yellow, blue and black in Fig. 4a), while at increasing depths phyloregions become more homogenous across latitude (blue and black in Fig. 4a).

Eco-evolutionary hypotheses behind phylogenetic patterns
Our results generally support our third hypothesis (Fig. 1eii): both energy/temperature and long-term stability act in concert as drivers of biodiversity, and energy/temperature is most likely not only influencing extinction but also speciation, at least for the percomorphs. However, the complex relationships of the phylogenetic diversity patterns along depth and latitude that deviate from the conceptual models suggest that this hypothesis is too naïve to fully explain the observed patterns and other factors undoubtedly play a role. Discrepancies between the observed and expected patterns may reflect uncertainty (species composition, phylogenetic relationships) inherent in observational studies and/or oversimplification of our theoretical models. For instance, all species are assumed to have the same physiological limits in our conceptual framework, but due to phylogenetic inertia some lineages may be unable to adapt to the environmental conditions of the deep sea and may be filtered along spatial gradients, leading to highly complex and unanticipated patterns of phylogenetic diversity. Furthermore, our conceptual framework considered only variation along the depth gradient, whereas our results (Fig. 2) indicate that patterns of phylogenetic diversity also vary with latitude. Deviations from our simple conceptual models point the way towards novel hypotheses regarding the distribution of phylogenetic diversity at large spatial scales.

The shallow ocean is a cradle of diversity for percomorph fishes
Our results (namely, that percomorphs were dominant in shallow waters with high SR and PD, and low MPD, MNTD and VPD) are consistent with the known prevalence of percomorph fishes in coastal and warm tropical shallow waters (e.g. central Indo-Pacific, Price et al. 2015, Miller et al. 2018 and the general increase in coastal marine diversity in warm marine regions (Tittensor et al. 2010). Shallow and lowlatitude regions are thought to have high speciation rates caused by: 1) warm, energy-rich waters that may speed up the evolutionary process and maintain viable population sizes (reducing extinction); 2) increased habitat diversity providing opportunities for diversification (Alfaro et al. 2007, Price et al. 2015, Siqueira et al. 2016; and 3) sea-level variations during glacial and inter-glacial periods (e.g. eustatism during the Pleistocene) that can restrict gene flow among populations to create incipient species (Cowman 2014, Ludt and Rocha 2015, Leprieur et al. 2016a. Recent studies on fishes tend to reject the hypothesis that elevated speciation rates are the main driver of species richness in shallow and tropical waters (Price et al. 2015, Miller et al. 2018, Rabosky et al. 2018. Instead, it is suggested that the older age of tropical ecosystems (e.g. central Indo-Pacific) and their long-term climatic stability (34 My) has allowed more time for speciation and colonization, along with decreased extinction risk (Pellissier et al. 2014, Leprieur et al. 2016b, Miller et al. 2018. However, prior studies have not considered depth gradients (but see Rabosky et al. 2018, O'Hara et al. 2019, where there is a sharper decrease in energy with depth than is found with increasing latitude. In our study, we cannot rule out the simultaneous contributions of habitat diversity, instability and fragmentation as potential drivers of speciation in shallow seas. Nonetheless, the rapid decline of SR, PD and the increasing time lag between speciation events with increasing depth (except at the deepest depth we sampled), suggest that energy might also be an important driver of diversity along the depth axis of marine systems (Woolley et al. 2016).
Disentangling the relative contribution of speciation and extinction in the diversification process is difficult without an abundant fossil record (Rabosky 2010). Nevertheless, our findings support the idea that the enhanced diversity found in shallow waters (50 m) could be the result of increased recent speciation, predominantly in percomorphs (Supplementary material Appendix 2 Fig. A12). Superficially, such a conclusion may seem to contradict recent studies that found the highest rate of speciation in deep and polar waters (Rabosky et al. 2018, O'Hara et al. 2019). However, we observed the lowest values of MNTD at 1200 m indicating that recent speciation is highest in the deep sea, followed by shallow waters at 50 m (Fig. 2d). Recent studies (Rabosky et al. 2018, O'Hara et al. 2019) used a coarser resolution of depth strata (i.e. above and below 200 m) that may have failed to detect raised speciation rates at shallow depths. Nonetheless, the mean DR statistic and relative PD (corrected for species richness) for the shallowest strata in tropical and temperate waters reported by O'Hara et al. (2019) indicated high recent speciation rates that were comparable to deep polar waters, suggesting speciation rates at shallow depths might only be enhanced in conjunction with warm photic surface water.

Intermediate ocean depths are a museum of diverse phylogenetic lineages
Our findings suggest that intermediate depths of the lower mesophotic zone (500-900 m) are a museum rather than a cradle of diversity, as proposed elsewhere (Brown and Thatje 2014, Gaither et al. 2016, Costello and Chaudhary 2017. In our study, communities were species-poor but represented a wide diversity of phylogenetic lineages characteristic of both shallow (Percomorpha) and deep (Gadiformes, Aulopiformes, Notacanthiformes) regions leading to high MPD values. Such high MPD values can result from either a high or low speciation rate along with a low to moderate extinction rate (Supplementary material Appendix 1 Fig. A5, A7) making it difficult to define the explicit contribution of speciation. Nonetheless, the high MNTD values found in parallel with high MPD values at intermediate depths (Fig. 2d) suggest that low speciation and extinction typify these depths, characteristic of a museum. Intermediate depths might provide a transitional zone where the lower and upper physiological boundaries for shallow and deep-sea specialists, respectively, overlap (Carney 2005, Brown andThatje 2014). Moreover, these fish communities lacked closely related species -a pattern previously described for brittlestars (O'Hara et al. 2019) -indicating that other factors may limit species diversity at these depths. Intermediate depths are characterized by rapid declines in temperature and particulate organic matter (Levin 2003), but also the presence of oxygen minimum zones (Levin 2003;although unlikely in New Zealand, Chiswell et al. 2015). This may present a challenging environment for fishes, with limited resources and increased competition among close relatives (Priede and Froese 2013). However, over a long period of time, the mesophotic zone might act as an important refuge for diverse lineages; both shallow and deep-sea specialists may survive at intermediate depths when faced with extreme events, such as climatic variation in the shallows (e.g. Pleistocene climatic oscillations; Moffitt et al. 2015) or rarer and ancient (Mesozoic) oceanic anoxic events (OAE) in the deep sea (Priede and Froese 2013, Song et al. 2014, Rogers 2015, enabling subsequent post-event recolonization (Song et al. 2014). As a potential 'museum', intermediate depths form a 'phylogenetic diversity bank' that could be crucial in the future to mitigate anthropogenic threats affecting shallow and/or deep-sea marine ecosystems (Rogers 2015, Danovaro et al. 2017).

The deep sea as a filter for ancient fish lineages
The deepest fish communities sampled here resembled a selective museum of old phylogenetic lineages that have undergone recent speciation (also suggested by Zintzen et al. 2011Zintzen et al. , 2012. Given that marine diversity most likely originated in shallow coastal waters (Sallan et al. 2018), several clades of fishes must have invaded the deep sea early in their evolutionary history (Haedrich 1996, Priede and Froese 2013, Sallan et al. 2018) and subsequently radiated (Barros-García et al. 2018). Potentially strong phylogenetic inertia for the adaptive traits required to cope with food scarcity, extreme pressure and low temperatures at depth may explain the low prevalence of recently derived, shallow-origin fish lineages in the deep sea (Priede 2017). Differential filtering of phylogenetic lineages may explain the increase of VPD and deviations from the relationships expected in our conceptual framework (Fig. 1), where all species were considered equal. The phylogenetic structure of deep-sea communities furthermore suggests that ancient actinopterygian lineages might have suffered extinction in the past (i.e. long internal branches). It is unclear whether old deep-sea lineages survived the Mesozoic massextinction events within in situ deep-sea refugia (Guinot et al. 2013, Thuy et al. 2014 or if they recolonized the deep sea from intermediate depths (Priede andFroese 2013, Priede 2017). Our results showed that, on average, closely related species at 1200 m diverged around 45 Ma ago, which postdates the last global OAE (OAE2, 96 Ma ago, Takashima et al. 2006) and one of the most recent regional OAE (56 Ma ago, Dickson et al. 2014), after which the deep-sea environment has remained relatively stable (Priede and Froese 2013).
Phylogenetic structures of fishes at depth suggested recent, and potentially high, rates of speciation. This pattern supports other recent findings of high speciation rates in deep-sea environments or cold polar regions for fishes, ophiuroids and octocorals (Chrysogorgiidae) (Pante et al. 2012, Bribiesca-Contreras et al. 2017, Rabosky et al. 2018, O'Hara et al. 2019, and in energy-poor environments of terrestrial systems (Weir and Schluter 2007, Schluter and Pennell 2017, Quintero and Jetz 2018. The increase of recent speciation rates in depauperate, cold and polar regions strongly suggests that energy might not be the key driver of speciation in marine systems (Rabosky et al. 2018, O'Hara et al. 2019). The instability of polar seas during the Cenozoic, including several cooling phases that formed ice-caps (41-34 Ma and 14 Ma, Crame 2018), may have driven extinctions, leaving vacant niches that fostered recent speciation (Crame 2018, O'Hara et al. 2019. In temperate deep-sea environments, perhaps long-term stability, coupled with localised physical seafloor disturbances which create spatial heterogeneity (Levin et al. 2001) increased the potential for specialization and speciation. These recently formed hypotheses demand further examination; the depth gradient may yet prove to be highly useful for decoupling the relative roles of environmental stability versus energy.

Conclusion
Our development of a simple conceptual framework, and associated description of multiple phylogenetic diversity patterns across broad depth and latitudinal gradients for marine ray-finned fishes (Actinopterygii) offers new insights regarding potential links between energy and long-term stability with the evolutionary mechanisms driving biodiversity. Our findings support the hypothesis that high energy is not a strict requirement for high speciation rates (Rabosky et al. 2018). Furthermore, long-term environmental stability might be crucial to preserve phylogenetically distinct lineages of fishes and thus buffer their extinction (Mittelbach et al. 2007). Last, we found that mean phylogenetic distances among species peaked at intermediate depths, despite their being species poor. This was primarily due to the co-occurrence of shallow and deep specialists, suggesting that intermediate depths might have provided a refuge against frequent climatic instability in shallow waters and rarer deep-sea extinction events (Song et al. 2014).
From the complexity of the phylogenetic diversity patterns we observed, it is apparent that several other abiotic and biotic factors play a role in shaping phylogenetic structures of fish communities, not only through their effects on speciation and extinction processes, but also on dispersal. Scenarios and hypotheses presented here neglect the crucial role of dispersal (Miller et al. 2018), which might be investigated through in-depth documentation of patterns in phylogenetic beta-diversity (Graham and Fine 2008). A better understanding of eco-evolutionary mechanisms shaping distributions of fishes will require integrated studies of taxonomic, phylogenetic and functional diversity, as well as emerging mechanistic models such as the spatial diversification model of lineages through time (SPLIT, Leprieur et al. 2016b, Descombes et al. 2018, Pontarp et al. 2019). Such methods model speciation, extinction and dispersal in spatially-explicit and temporally dynamic environments (e.g. tectonic movements, Leprieur et al. 2016b, Descombes et al. 2018. However, disentangling the contribution of energy and environmental stability along continuous gradients will also require consideration of the evolvability and heritability of species' niche breadths through time (Brayard et al. 2005), allowing for phylogenetic niche conservatism, peripatric/parapatric speciation and selective extinction. In addition, to test for the potential indirect contributions of energy, such mechanistic models could be expanded to include diversity-dependent speciation models driven by ecological opportunities (Etienne et al. 2019), indexed on the amount of available energy. We consider that our conceptual framework and quantitative tests provide a first step toward the development of these more elaborate approaches. Integrating empirical and theoretical models in this way, we can gain a greater understanding of biodiversity determinants for the largest habitat on earth, which will be crucial to achieve better management of ocean ecosystems under increasing future anthropogenic threats (Rogers 2015, Danovaro et al. 2017.