Digital Commons @ Michigan Tech Digital Commons @ Michigan Tech Evolution of realized Eltonian niches across Rajidae species Evolution of realized Eltonian niches across Rajidae species

. The notion that closely related species resemble each other in ecological niche space (i.e., phylogenetic dependence) has been a long-standing, contentious paradigm in evolutionary biology, the incidence of which is important for predicting the ecosystem-level effects of species loss. Despite being examined across a multitude of terrestrial taxa, many aspects of niche conservatism have yet to be explored in marine species, especially for characteristics related to resource use and trophic behavior (Eltonian niche characteristics, ENCs). We combined ENCs derived from stable isotope ratios at assemblage- and species-levels with phylogenetic comparative methods, to test the hypotheses that benthic marine ﬁ shes (1) exhibit similar assemblage-wide ENCs regardless of geographic location and (2) display phylogenetically dependent ENCs at the species level. We used a 12-species sub-set of the monophyletic group Rajidae sampled from three independent assemblages (Central California, Gulf of Alaska, and Northwest Atlantic), which span two ocean basins. Assemblage-level ENCs implied low trophic diversity and high evenness, suggesting that Rajidae assemblages may exhibit a well-de ﬁ ned trophic role, a trend consistent regardless of geographic location. At the species level, we found evidence for phylogenetic dependence of ENCs relating to trophic diversity (i.e., isotopic niche width; SEAc). Whether individuals can be considered functional equivalents across assemblages is hard to ascertain because we did not detect a signi ﬁ cant phylogenetic signal for ENCs relating to trophic function (e.g., trophic position). Thus, additional, complimentary approaches are required to further examine the phylogenetic dependence of species functionality. Our approach illustrates the potential of stable isotope-derived niche characteristics to provide insight on macroecological processes occurring across evolutionary time, which could help predict how assemblages may respond to the effects of species loss.


INTRODUCTION
The extent to which closely related species share similar ecological characteristics, such as morphology, thermal tolerance, habitat, and diet, herein defined as phylogenetic niche conservatism (PNC, Cooper et al. 2010), has long fascinated evolutionary biologists (Harvey and Pagel 1991, Holt and Barfield 2008, Losos 2008, Wiens et al. 2010. Understanding the extent to which species resemble each other in ecological space can reveal insights into the complex interplay between evolutionary and biological process, which ultimately dictate the distribution and functional roles of species in space and time (Losos et al. 2008(Losos et al. , M€ unkem€ uller et al. 2012. Methodological advancements in the past two decades have revolutionized our understanding of niche evolution in wild populations by adopting a phylogenetic approach (Pagel 1999, Blomberg et al. 2003, M€ unkem€ uller et al. 2012. Despite the proliferation of phylogenetic comparative methods, much of the published work has been limited to examining niche evolution in terrestrial fauna, such as birds (Rice et al. 2003, Maldonado et al. 2017, lizards (Harmon et al. 2003, Knouft et al. 2006, and mammals (Olalla-T arraga et al. 2017), with fewer attempts having been made in marine organisms, including fishes (e.g., Ingram et al. 2010, Egan et al. 2018). This emphasis on terrestrial fauna is likely an artifact of the logistical challenges of studying marine organisms, especially those which display complex and cryptic life histories. Considering the multitude of threats facing marine ecosystems (C ozar et al. 2014, McCauley et al. 2015, Molinos et al. 2016, understanding how ecological niches have evolved is intricately tied to our ability to identify vulnerable taxa, predict how associated distributions and functional roles are likely to change over time, and understand how protecting specific niches may help preserve marine biodiversity (Brooker et al. 2016).
The study of Eltonian niche characteristics (ENCs), defined as those encompassing resourceconsumer dynamics and their associated trophic behavior and diversity, can be used to evaluate the functional role of species within a food web (Elton 1927, Sober on 2007, Cooper et al. 2010. Eltonian niche characteristics are integral components of the fundamental niche (Hutchinson 1957) that ultimately dictates the fine-scale distributions of species in space and time. Assessments of PNC for ENCs are comparatively scarce in contrast to those investigating Grinnellian niche axes (e.g., thermal niches, geographical ranges, etc.). Existing evidence for PNC in ENCs is equivocal largely due to the sensitivity of phylogenetic comparative methodologies to the specific trait measured, the spatiotemporal variation in the trait, and its taxonomic representation (Blomberg et al. 2003, Losos et al. 2008, Revell 2010. In fact, many have argued that our ability to detect the phylogenetic signal of niche traits may be heavily biased by the specific model of evolution assumed (i.e., Brownian, niche filling, accelerating vs decelerating; Cooper et al. 2010). Thus, there is a need for examination of PNC in ENCs across a diversity of measured traits at multiple phylogenetic resolutions.
Typically, phylogenetic assessments of ENCs have been examined using stomach content data (e.g., Olalla-T arraga et al. 2017). Work by B€ ohning- Gaese and Oberrath (1999) found evidence for dietary niche conservatism in European birds. However, more recent work by Olalla-Tarraga et al. (2017) used a global dataset of mammalian diets and found limited evidence for Eltonian niche conservatism. Among-study differences of Eltonian niche estimates may be due to the use of stomach content data, which generally provides a snapshot view of diet, and may reflect shortterm variability in foraging behavior. Therefore, additional data types that allow for integration of prey consumption over time may be required to fully capture the multi-dimensional nature of resource use and trophic dynamics in space and time.
Ecogeochemical proxies, which vary as a function of species' biotic and abiotic phenotypes, offer an appealing means for examining niche dynamics in space and time (Bearhop et al. 2004, Newsome et al. 2012. Specifically, the variability of naturally occurring stable isotope ratios in the tissues of organisms are strong indicators of habitat-use (McCauley et al. 2012, Bird et al. 2018, ambient temperature (Power et al. 2003, diet Epstein 1978, 1981), and the physiological processes driving fractionation and turnover (Pinnegar andPolunin 1999, Gorokhova 2018 usefulness of stable isotope ratios lies in the fact that they encompass a hypervolume of ENCs, the variability of which can be analyzed in geometrically multi-dimensional space (e.g., Jackson et al. 2011, Swanson et al. 2015, Junker et al. 2016, Blonder et al. 2018. Few studies have applied stable isotope ratios within a phylogenetic context, but those that have present conflicting support for PNC. For example, Maldonado et al. (2017) found that closely related perching passerine birds exhibited similar degrees of individual foraging specialization (as inferred through d 15 N). Conversely, Manlick et al. (2019) used regionally scaled isotopic compositions (d 13 C and d 15 N) of American and Pacific martens (Martes americana and M. caurina) sampling from four regions across North America and found limited evidence of PNC. Limited data coupled with equivocal support for PNC suggest a need to explore the phylogenetic underpinnings of isotopically inferred niche metrics. Furthermore, efforts such as these could provide critical insight into the complementarity of species in a food-web context. Rajiformes (skates) are a speciose family of cartilaginous marine fishes (Class: Chondrichthyes) comprising~293 extant species (Stein et al. 2018) that first diverged between 157 and 208 million years ago (Aschliman et al. 2012). Compared to other chondrichthyans, such as sharks, skates exhibit relatively low evolutionary distinctiveness (Stein et al. 2018), further evidenced by a cryptic taxonomy (Bizzarro et al. 2014). Recent work suggests that closely related skate species exhibit only subtle differences in niche characteristics, such as habitat associations (Zacharias 2013) and aspects of their trophic niche (Shipley et al. 2019). Yet there have been very few examinations of Eltonian niche evolution in skates, despite the usefulness of such information for understanding distributions, inter-species interactions, and the subsequent roles of species in ecological communities. This lack of attention stems from a generic paucity of biological information pertaining to skates as a whole, which is problematic considering the growing concern for the future vitality of populations due to exploitation (Dulvy andReynolds 2002, Dulvy et al. 2014) and the multitude of emerging threats marine communities now face (Knip et al. 2010).
Here, we investigated the assemblage-and species-level similarity of ENCs displayed by three skate assemblages spanning two ocean basins. We hypothesize that (1) skates will exhibit similar ENCs at the assemblage-level and (2) species-specific ENCs will be phylogenetically dependent. This study contributes to the growing body of literature specific to understanding trophic complementarity and the functional roles of cryptic species within skates. Further, it provides a clear blueprint for combining phylogenetic comparative methods with ecogeochemical tracer techniques to test a fundamental hypothesis regarding the evolution of the ecological niche, a necessary perspective for formulating predictive solutions to global biodiversity loss.

METHODS
Carbon and nitrogen stable isotope ratios generated from white-muscle tissue of 12 skate species (n = 632) were compiled from three independent datasets: Central California, the Gulf of Alaska, and the NW Atlantic continental shelf. Muscle tissue was selected due to its relatively long isotopic turnover rate, which reflects a temporally integrated average (~1 yr, MacNeil et al. 2006, Kim et al. 2012, Hussey et al. 2012) of resource use and is thus insensitive to small-scale (seasonal) fluctuations in trophic behavior and/ or movement.
Skate white-muscle tissue samples from central California were collected during bottom trawl and longline surveys conducted by the National Marine Fisheries Service (NMFS) Southwest Fisheries Science Center between 2002 and 2005 (Bizzarro et al. 2014). Sampling was conducted inside of Monterey Bay and on the northern adjacent continental shelf and slope (36.0°-7.5°N, Fig. 1 Samples from all datasets were treated and/or corrected for lipids and urea, which are enriched in 12 C and 14 N relative to pure protein, respectively (Sweeting et al. 2006, Carlisle et al. 2017, Shipley et al. 2017. Accounting for lipids and urea in this context improved the overall comparability of Eltonian traits between the three assemblages. A detailed description of sample preparation and stable isotope analysis can be found in Appendix S1.

Assemblage-level ENCs
All statistical analyses were performed in RStudio (version 1.1.383; R Core Team 2017).
Metrics first applied by Layman et al. (2007) were used to quantify the assemblage-level ENCs of skates across the three independent assemblages. For these calculations, species sampled from multiple assemblages were treated as distinct populations because they were captured in geographically distinct regions. Metrics reflect the spatial configuration of species bivariate means in delta-space and include carbon and nitrogen ranges (CR, NR), total convex hull area (TA), mean distance to the assemblage centroid (CD), mean nearest neighbor distance (MNND) and standard deviation of MNND (SDMNND). CR and NR indicate the range of primary production sources, as well as the trophic diversity exhibited by species within an assemblage. TA represents the total area of niche space occupied by the bivariate species means, thus the combined variability of habitat (d 13 C) and prey resources (d 15 N) utilized by an assemblage. The CD is a measure of trophic diversity defined by the mean Euclidean distance between bivariate species means and the assemblage centroid (bivariate mean for the entire assemblage). The MNND defines the relative degree of species trophic relatedness (i.e., packing, Layman et al. 2007) and is calculated as the mean Euclidean distance between each species' nearest neighbor. Finally, SDMNND accounts for sample size biases in MNND estimates and represents the relative degree of species evenness (Layman et al. 2007).

Species-specific ENCs
Eight ENCs were calculated for each species, these reflect the extent and direction of isotopic dispersion in delta-space, which is ultimately Eccentricity is a value between 0 and 1 and reflects the contribution of one or both isotopes to the ellipse (Turner et al. 2010, Reid et al. 2016, where E approaches 1 variance in delta-space is driven by a single isotope, and as E approaches 0 variation is equally driven by both isotopes. Theta reflects the slope of the relationship between d 13 C and d 15 N relative to the x-axis (i.e., d 13 C). Theta can be positive or negative depending on the direction of isotopic dispersion, that is, positive or negative covariance between d 13 C and d 15 N (Reid et al. 2016). We also calculated species-level estimates of MNND, SDMNND, and CD using the R package SIBER (Jackson et al. 2011). Carbon and nitrogen ranges (CR and NR) were not calculated because estimates are extremely sensitive to variable sample sizes (Layman et al. 2007). Variation in the regional isotopic composition of baseline producers, or isoscapes, occurs due to nutrient input, dominating phytoplanktonic composition, and oceanographic conditions (McMahon et al. 2013, Magozzi et al. 2017, Lorrain et al. 2019, and precludes the direct comparison of raw d 13 C and d 15 N between skate assemblages (Manlick et al. 2019). Therefore, we calculated region-specific estimates of trophic position (TP) using the R package tRophicPosition (Quezada-Romegialli et al. 2018). This approach combines two end-member mixing models to calculate relative trophic position using a Bayesian inference, which allowed for the incorporation of region-specific, pelagic and benthic isotopic baselines. We also calculated the mixing parameter a which represents the contribution of Baseline 1 (pelagic production) to consumer TP. Literature-derived isotope data used to assign baseline taxa and their associated sampling information (date, year, etc.) can be found in the Appendix S1 (Appendix S1: Table S1). Posterior distributions for TP and a were derived from 10,000 Monte Carlo simulations, run across two chains, each with a burn-in period of 1000. We used trophic discrimination factors (i.e., the enrichment/depletion of 13 C/ 15 N predatory tissues relative to their prey, TDFs) reported for leopard (Triakis semifasciata) sharks by Kim et al. (2012) where D 13 C = 1.7 AE 0.5& and D 15 N = 3.7 AE 0.4&. Mixing models and trophic position estimates can be sensitive to variable TDFs, so we ran sensitivity models for TP and a in which D 13 C and D 15 N were reduced by 0.5&.
We examined the relative correlation of ENCs through principal component analysis (PCA, R package vegan, Oksanen et al. 2019), which also provided an opportunity to calculate relative niche envelopes for each species. This is because the principal components (i.e., scores derived for each species), based on correlation between ENCs, provided a single trait value that could be tested for phylogenetic signal (see Schnitzler et al. 2012, Pienaar et al. 2013. Others have combined species' phylogenetic distances with measures of niche differentiation (e.g., relative overlap of niche envelopes) using statistical approaches such as Mantel tests (e.g., Knouft et al. 2006). However, this approach was not appropriate because phylogenetic distance and measures of niche differentiation were derived from different distance measures (i.e., Euclidean vs. patristic). We constructed Eltonian niche envelopes for each species by first normalizing all metrics to a mean of zero and a standard deviation of 1. We then extracted the principal component scores for each species, across those axes that explained more variance than expected at random, to carry forward to phylogenetic signal analysis.

Phylogenetic signal analysis
A phylogenetic tree comprising the 12 study species was constructed using the R package Ape (Paradis and Schliep 2018) using the mitochondrial cytb, 12S, and cox1 genes (Fig. 2); full details of phylogenetic tree construction can be found in the supplementary materials (Appendix S1: Table S2).
Because two species were captured at multiple locations (Beringraja binoculata and Beringraja rhina), we performed two analyses, one using ENCs from Central California individuals, and a second using those from the Gulf of Alaska. We tested the phylogenetic signal of species-level ENCs and niche envelopes derived from PC1, PC2, and PC3 by calculating k (Pagel 1999, Freckleton et al. 2002  2003) using the R package phytools (Revell 2012). Both parameters were designed to quantitatively assess the phylogenetic structure of trait data and can be used to compare differences in traits across a phylogenetic tree to discover general patterns of relative evolutionary lability across trait types (Pagel 1999, Blomberg et al. 2003. Lambda is a scaling parameter that compares the similarity between species covariance, to the covariance expected under a random walk model (Pagel 1999), whereas K is given by the ratio: (MSE 0 /MSE) observed /(MSE 0 /MSE) expected, where MSE0 is the mean squared error of the trait data measured from the phylogenetically corrected mean, and MSE is the mean squared error of the variance covariance matrix given the phylogeny under the Brownian motion evolution (M€ unkem€ uller et al. 2012). We calculated both k and K; however, K is more sensitive to behavioral traits, which are generally considered to be more evolutionarily labile (Blomberg et al. 2003, Cooper et al. 2010. Although these parameters are calculated differently, a value of 0 for both k and K indicates that evolution of trait variation among species is not correlated with the phylogeny, and instead closely approximates the null model (i.e., a star tree that assumes equal relatedness of all taxa). Alternatively, higher values of k and K suggest greater phylogenetic dependence in a trait, a pattern consistent with the predictions of PNC (Cooper et al. 2010, Kamilar andCooper 2013). Both parameters can be used in a hypothesis testing framework to derive P values by running the observed trait data against a null model, where k is based on a likelihood ratio test and K is based on a randomization test. To visualize the variation of normalized ENCs across the phylogeny, a heatmap was generated across two scenarios using the R package phytools (Revell 2012).

Assemblage-level ENCs
Carbon ranges were consistent and did not exceed 1.1&, but nitrogen ranges were greater for Central California (1.32&) than for Gulf of Alaska and NW Atlantic assemblages (<1.0&, Table 1). The total convex hull area encompassed by each assemblage was consistently low and was <0.54& 2 for all assemblages. Assemblages  (Table 1).

Species-specific ENCs
Species-specific ENCs suggested high specieslevel variation within and among assemblages (Tables 1, 2). Standard ellipse area (SEAc) estimates were generally lowest in the Gulf of Alaska assemblage (0.31& 2 -0.99& 2 ) and highest in the NW Atlantic species (1.06& 2 -2.47& 2 ). These patterns were also reflected in MNND, SDMNND, and MCD (Table 2). Trophic position estimates varied across species and by region (Fig. 3). The highest TP estimates were observed by species captured from the Gulf of Alaska (TP SI 4.25-4.47), followed by the NW Atlantic (TP SI = 3.71-3.92), with Central California species exhibiting the lowest trophic position estimates (TP SI 2.90-3.25, Fig. 3). The Gulf of Alaska and NW Atlantic skate assemblages were predominantly supported by pelagic production pathways (a Pelagic = 0.69-0.92), with a smaller contribution of benthic-derived production. The Central California skate assemblage was supported by a strong mix of both pelagic and benthic production (a Pelagic = 0.52-0.70, Fig. 3). Across all species, trophic position and a estimates did not appear sensitive to the use of two different DTDFs, with the sensitivity model increasing trophic position estimates by 0.16-0.37 and reducing alpha values by 0.05-0.11 (Appendix S1: Table S3).

Phylogenetic signal analysis
Phylogenetic signal analysis highlighted that SEAc was phylogenetically dependent and that this trend was consistent across both scenarios for B. binoculata and B. rhina (Table 4, Fig. 5). For the Central California scenario, SEAc yielded k = 0.83 and K = 0.48 for the Gulf of Alaska scenario k = 0.99 and K = 0.93. In both scenarios, k and K were statistically significant compared to the null model (P < 0.05, Table 4). In some cases, ENCs were only found to be phylogenetically dependent in one of the two trait scenarios. Mean distance to centroid (CD) was found to carry a phylogenetic signal in one of the two scenarios. For the Central California assemblage, k = 0.7 and was approaching statistical significance when compared to the null model (P = 0.06), K = 0.31, however, was not found to be significantly different relative to the null model (P = 0.11). For the Gulf of Alaska scenario k = 0.006 and K = 0.75, both of which were found to be significantly different than the null model (P = 0.006, P = 0.003, respectively, Table 4). This was also true for a under the Central California scenario across both k (k = 1, P = 0.05) and K (K = 0.59, P = 0.002), and for PC3 for K (K = 0.4, P = 0.03; Table 4; Fig. 5). All other ENCs were not found to carry a phylogenetic signal across either of the scenarios (Table 4).

DISCUSSION
Assemblage-level ENCs suggested that skates exhibit similar ecological roles and that this trend is consistent across broad geographic scales. At the species level, niche characteristics pertaining to resource use diversity (SEAc) showed a significant phylogenetic signal suggesting that closely related species may share some similar aspects of their Eltonian niches under the expectations of Brownian motion evolution. Linking support for or against PNC with a specific evolutionary process has been strongly discouraged (Revell et al. 2008, Cooper et al. 2010); thus, we discuss the potential drivers that may contribute to the prevalence of a phylogenetic signal in some traits and the potential diversification of ENCs between less closely related species.
Assemblage-level ENCs for skates were found to be similar across geographic locations, suggesting the cumulative ecological role of skate assemblages are relatively uniform within their respective ecosystems. Assemblage-level ENCs, however, have been criticized for providing insufficient resolution to provide ecological meaningful information, because trophic inferences may be diluted by using bivariate species Notes: P values represent results from likelihood ratio (k) and randomization tests (K) to test for statistical significance against the null hypotheses (k = 0, K = null expectation, Swenson 2014). Bold indicates statistical significance at a = 0.05. Models were run across two different scenarios due to the presence of B. binoculata and B. rhina in both Central California and Gulf of Alaska assemblages. v www.esajournals.org 10 February 2021 v Volume 12(2) v Article e03368 means. It is likely that the ecosystem function of individual skate species is not synonymous and may be augmented by different levels of variation occurring at the individual-level. This is illustrated by the variable trophic position estimates observed here (TP = 2.9-4.5), which are normalized to the regionally specific isotopic baselines of each species' sampling location. The strong phylogenetic signal of ENCs describing resource use diversity (SEAc) may reflect the environmental conditions and geographic ranges over which species could have evolved trophic behaviors (Van Valen 1965, Maldonado et al. 2017. Spatial scale and environmental homogeneity are known to drive resource use diversity as measured using isotopic approaches (i.e., the resource breadth hypothesis, Rader et al. 2017, Reddin et al. 2018. For example, Reddin et al. (2018) found that the spatial scale explained up to 50% of isotopic variance in three intertidal invertebrates (Mytilus spp., Patella vulgata, and Littorina littorea). Here, we find that the broadest isotopic niches are exhibited by species from the NW Atlantic, the study system where individuals were captured across the greatest latitudinal range (and thus encompass the largest isotopic baseline variation). It is therefore plausible that spatial extent may influence resource use diversity of skates and could in part explain the strong phylogenetic signal for these ENCs, if closely related species move and forage across similarly heterogenous landscapes. Despite this hypothesis, we did not observe phylogenetic a signal for MNND, SDMNND, and CD raising questions over whether individuallevel patterns of Eltonian niche variation are in fact consistent across closely related species.
It is not possible to ascertain whether the phylogenetic dependence of ENCs related to resource use diversity translates to occupation of similar functional roles by skates in a broader food-web context. In fact, we found little support for phylogenetic signal in Bayesian-derived TP estimates, suggesting that closely related species could serve different functional roles. This is evidenced by the Gulf of Alaska assemblage, where Bayesian mixing models revealed that individuals were a TP higher than Central California and NW Atlantic species, which likely represents the high food-web complexity and heightened influence of the detrital food web in Gulf of Alaska relative to Central California (Gaichas et al. 2012, Madigan et al. 2012. Additionally, the two species found across both the Gulf of Alaska and Central California exhibited different TPs dependent upon location, suggesting that these two species could exhibit different functional roles in different ecosystems (Manlick et al. 2019). One explanation for this observation may relate to the potential for some ENCs to be inherently labile (Grundler et al. 2017, Manlick et al. 2019, thus responding to processes at a different rate to that expected under Brownian motion evolution. In fact, spatiotemporal variation of ENCs is common among studied chondrichthyan populations (MacNeil et al. 2005, Madigan et al. 2015, and the incidence of specialization across many populations (Matich et al. 2011, Munroe et al. 2014 suggests ENCs are by nature, extremely changeable. Despite limited assessments of phylogenetic dependence of ENCs across marine fishes, data from carnivores and passerine birds support the idea of Eltonian niche lability, both within and between populations (e.g., Maldonado et al. 2019, Manlick et al. 2019. One may hypothesize that the chance of Eltonian traits evolving in line with a strictly Brownian process, such as natural selection, is oversimplified, and other models may be more effective at validating phylogenetic signal (e.g., niche filling, Freckleton and Harvey 2015). Regardless, the ability of skates to potentially alter aspects of their Eltonian niche in response to fluctuating environmental conditions could be a driving force behind the long persistence of the lineage across evolutionary time and may be advantageous considering environmental change and continued anthropogenic exploitation of natural resources (Lawton et al. 2012). Certainly, significant ontogenetic changes in trophic behavior has been observed in species such as big and longnose skate (Bizzarro et al. 2007, which could potentially buffer species complexes in rapidly changing environments. Debate remains over the use of stable isotope ratios to accurately characterize Eltonian niches in ecological systems, mainly due to the biotic and abiotic processes that can strongly influence the bulk isotopic composition of organismal tissue, such as (reviewed by Newsome et al. 2007, Shipley andMatich 2020). For example, non-dietary determinants such as individual temperature (Sweeting et al. 2007, Britton andBusst 2018) and environmental stress (Reddin et al. 2016) may partially obscure an exclusively Eltonian signal, and these facets may vary across species and among environments. It is also possible that consumers may utilize different prey items that share a similar isotopic composition, thus divergent feeding behaviors are otherwise diluted by relying upon the bulk composition of a consumer (Hammerschlag 2019). We removed much of this bias using assemblage-specific TPs, which accounted for interregional variability in isotopic baselines, and drivers therein. Additionally, our two-scenario approach allowed us to account for the potential variability of ENCs (and associated drivers) in the two species found across multiple sampling regions (B. binoculata and B. rhina), where we found that traits relating to resource use diversity were phylogenetically dependent across both scenarios. This finding allowed us to remove much of the bias within the phylogenetic models that could be introduced by generating ENCs from a single assemblage/location. In some cases, we found phylogenetic signal for traits derived from one scenario, but not the other, for example, a. It must be noted that our sampling design is somewhat opportunistic in that samples were collated from multiple studies, locations, and across different years, which may account for variability in predatory-prey dynamics and isotopic baseline composition (D ecima et al. 2013). While the isotopic composition of white-muscle tissue likely dampens short-term variability, due to a long isotopic turnover rate, it is possible that detection of a phylogenetic signal in some ENCs only in a single scenario may be explained by sampling design and interannual baseline variation. Additionally, slightly different sample pre-treatment methods between communities could be a source of variability in the stable isotope ratios (Connan et al. 2019). However, the above factors are unlikely to introduce the variability necessary to bias overall ecological interpretation of Eltonian niche diversity in skates. These drawbacks, however, emphasize the importance of accounting for the potential variability introduced by sampling design, something that future studies exploring phylogenetic dependence of ENCs using stable isotope tracers should also aim to incorporate (Shipley and Matich 2020).
The ability to accurately detect a phylogenetic signal can also depend on the size of the phylogeny as well as the quality of phylogenetic information available to infer trees (Cooper et al. v www.esajournals.org 12 February 2021 v Volume 12(2) v Article e03368 2010, M€ unkem€ uller et al. 2012). Phylogenetic comparative methods utilizing small trees with a low number of species can exhibit high type 1 and type 2 errors, which are hard to reduce without increasing taxon resolution (M€ unkem€ uller et al. 2012). Our study used a relatively small sub-set of species, which may have introduced biases into phylogenetic signal analysis (Blomberg et al. 2003) and could further explain why some ENCs carried a phylogenetic signal in one scenario, but not the other. Despite these noted limitations, others have tested for phylogenetic signal of various ecological traits in smaller phylogenies comprised of species sub-sets (Corcobado et al. 2012, Easson and Thacker 2014, Maldonado et al. 2017, Grundler et al. 2017). To refine the use of stable isotope ratios in a phylogenetic framework, sensitivity of phylogenetic signal analysis to clade size should certainly be examined. This type of assessment may be possible with the use of large chondrichthyan databases that allow for generation of complete phylogenies, such as the Chondrichthyan Tree of Life (https://sharksrays.org/), when used in association with global data sharing platforms that may be used to derive ENCs such as IsoBank (Pauli et al. 2017) and the Chondrichthyan Stable Isotope Data Project (CSIDP, Bird et al. 2018).

ACKNOWLEDGMENTS
We would like to acknowledge J. Bigman, A. Carlisle, S. Brown, S. Kim, D. Ebert, and Moss Landing Marine Laboratories for providing access to stable isotope data for the Gulf of Alaska and Central California assemblages. Thanks to New York State Department for Environmental Conservation (NYDEC) for supporting logistical support. ONS, JK, JAO, MP, and MGF conceived project concept. ONS, MGF, JJB, JAO, and MP conducted fieldwork and performed laboratory analyses. ONS, JK, RMC, and MGF performed statistical analyses. ONS wrote the paper with input from all co-authors. All authors approved the submitted version of the manuscript.