High functional diversity in deep‐sea fish communities and increasing intraspecific trait variation with increasing latitude

Abstract Variation in both inter‐ and intraspecific traits affects community dynamics, yet we know little regarding the relative importance of external environmental filters versus internal biotic interactions that shape the functional space of communities along broad‐scale environmental gradients, such as latitude, elevation, or depth. We examined changes in several key aspects of functional alpha diversity for marine fishes along depth and latitude gradients by quantifying intra‐ and interspecific richness, dispersion, and regularity in functional trait space. We derived eight functional traits related to food acquisition and locomotion and calculated seven complementary indices of functional diversity for 144 species of marine ray‐finned fishes along large‐scale depth (50–1200 m) and latitudinal gradients (29°–51° S) in New Zealand waters. Traits were derived from morphological measurements taken directly from footage obtained using Baited Remote Underwater Stereo‐Video systems and museum specimens. We partitioned functional variation into intra‐ and interspecific components for the first time using a PERMANOVA approach. We also implemented two tree‐based diversity metrics in a functional distance‐based context for the first time: namely, the variance in pairwise functional distance and the variance in nearest neighbor distance. Functional alpha diversity increased with increasing depth and decreased with increasing latitude. More specifically, the dispersion and mean nearest neighbor distances among species in trait space and intraspecific trait variability all increased with depth, whereas functional hypervolume (richness) was stable across depth. In contrast, functional hypervolume, dispersion, and regularity indices all decreased with increasing latitude; however, intraspecific trait variation increased with latitude, suggesting that intraspecific trait variability becomes increasingly important at higher latitudes. These results suggest that competition within and among species are key processes shaping functional multidimensional space for fishes in the deep sea. Increasing morphological dissimilarity with increasing depth may facilitate niche partitioning to promote coexistence, whereas abiotic filtering may be the dominant process structuring communities with increasing latitude.


| INTRODUC TI ON
Studying biodiversity across large-scale environmental gradients plays a key role in aiding scientists to understand potential mechanisms shaping species' distributions. Analyses of taxonomic diversity (e.g., species richness) are useful, but a more integrative understanding and insights regarding potential mechanisms driving biodiversity and ecosystem services can be obtained through analyses of phylogenetic diversity and functional diversity (Díaz et al., 2016;Swenson, 2011). Previous studies have documented a general decrease in species richness with increasing latitude, elevation, and depth (Costello & Chaudhary, 2017;Gaston, 2000;Hillebrand, 2004). However, functional diversity displays a variety of different patterns along gradients Stuart-Smith et al., 2013;Villéger et al., 2010).
Furthermore, trends in functional diversity will depend on the particular traits that are measured and the indices that are calculated from these.
Functional diversity is inherently multivariate, where each trait is a variable, and species (or individuals) occupy a particular position in multivariate trait space. There are many ways to measure functional diversity (Laliberté & Legendre, 2010;Mouchet et al., 2010;Villéger et al., 2008), and each metric may capture a different aspect of functional diversity. For example, one might measure the functional hypervolume covered by a set of species ("space") or the packing (i.e., the proximity of neighboring species) within this volume Swenson & Weiser, 2014). Patterns in these different aspects of functional diversity along environmental gradients vary, depending on the focal taxa and traits included, as well as the resolution and spatial scale of the study. For example, both functional richness and species richness were found to decline with increasing latitude for plants in the New World . In contrast, high species richness for birds at low elevations was characterized by a greater density (packing) of species in functional space rather than an increase in the volume (functional richness) of the space occupied (Pigot et al., 2016). Simultaneous analysis of multiple functional diversity measures, along with species richness, yields more comprehensive biodiversity information and insight into the potential processes driving community assembly.
External processes, such as abiotic filtering, and internal processes operating within the community, such as density-dependent interactions (e.g., competition), work in tandem to shape functional diversity (Kraft, Adler, et al., 2015;Swenson & Weiser, 2014;Violle et al., 2012). Abiotic filtering is assumed to select the number or type of functional strategies adapted to the environmental conditions encountered by a community (i.e., a species pool), determining its volume (or functional richness), whereas internal community processes such as competition are generally expected to affect the packing (or density) of species or individuals within the functional volume occupied by the community (Swenson & Weiser, 2014). For example, increasing competition for limited resources in environmentally harsh conditions is expected to decrease the packing of taxa in a functional space, whereby similar species are excluded from the community (Swenson & Weiser, 2014). Similarly, both external and internal drivers can also affect intraspecific functional variation (Violle et al., 2012).
Variation in traits among individuals within a species may be particularly important in species-poor regions, for taxa with narrow geographic ranges, or that live in stressful environmental conditions (Hoffmann & Merilä, 1999;Siefert et al., 2015). For many species, the deep sea represents an example of both a species-poor environment, and one with stressful environmental conditions (Priede, 2017), providing a unique opportunity to investigate patterns of inter-and intraspecific functional diversity along an increasing stress gradient. Disentangling individual-level variability versus species-level variation in functional trait space, including along environmental gradients, can provide novel insights into how communities are structured and maintained (Bolnick et al., 2011;Violle et al., 2012Violle et al., , 2014. The depth gradient is one of the steepest environmental gradients on earth, yet is one of the least studied in terms of functional diversity. Changes with increasing oceanic depth are dramatic, including decreasing light, temperature, dissolved oxygen, food, and increasing pressure (Priede, 2017). These changes strongly influence the spatial distribution of species, their functions, and morphologies (Mindel et al., 2016;Myers et al., 2019;Zintzen et al., 2017). One of the most striking patterns is a decline in species richness with increasing depth. Assuming that a transition from shallow to deep waters or from subtropical to sub-Antarctic latitudes represents a progression from a benign environment to an abiotically harsher, more extreme environment, we generated three contrasting conceptual models (a-c) regarding expected functional changes in communities along these gradients. Each (a) Abiotic filtering will be stronger in harsh environments, such as along increasing depth or latitude gradients (Swenson, 2011). This will lead to a decrease in functional hypervolume, but species packing should remain constant ( Figure 1a). (b) Biotic interactions, such as competition for limited resources, will intensify with increasing depth or latitude. These interactions will decrease the packing of species (or individuals) in functional space, but the overall volume K E Y W O R D S biodiversity, biotic interactions, deep-sea fishes, depth gradient, environmental filtering, functional trait, morphology, niche partitioning of the functional space will remain unchanged ( Figure 1b). (c) Both abiotic filtering and biotic filtering will jointly affect functional diversity, decreasing both the volume and the packing of species (or individuals) in deeper or higher latitude regions ( Figure 1c).
Here, we tested these predictions by quantifying patterns of change in functional diversity for 144 species of marine ray-finned fishes along large-scale depth (50-1,200 m) and latitudinal gradients (29°-51° S) in New Zealand waters. We obtained morphological measurements of individual fishes in situ using Baited Remote Underwater stereo-Video systems (stereo-BRUVs), allowing quantification of both interspecific trait variation and intraspecific trait variation. We focused on 8 functional traits related to food acquisition and locomotion and calculated 7 complementary functional diversity indices to examine their relationships with depth and latitude using univariate and multivariate approaches. By examining trends for a suite of different functional diversity metrics versus depth and latitude, we were able to characterize, more specifically, the type of change occurring in the functional space, and thus begin to differentiate among competing underlying models that might explain functional variation along these large-scale gradients.

| Fish community data
We observed 144 species of marine ray-finned fishes (Class Actinopterygii) on Baited Remote Underwater stereo-Video (stereo-BRUV) footage. Stereo-BRUVs were deployed in a stratified random sampling design, with n = 5-7 replicate units at each of seven targeted depths: 50, 100, 300, 500, 700, 900, and 1200 m, within each of seven locations spanning 21° of latitude in New Zealand waters. F I G U R E 1 Schematic diagram of three conceptual models of changes in functional space (and functional diversity metrics) when moving from a relatively benign (gray), to a more extreme environment (black). Points represent species (or individuals), and the squares represent the bounds of functional space. (a) Abiotic filters cause a decrease in functional volume, but species packing remains unchanged. (b) Biotic filters, such as increased competition for limited resources, do not affect volume, but decrease the packing of species (or individuals) in functional space. (c) Harsh abiotic conditions, decreased resources, and increased competition combine to yield a decrease in volume and also in the packing of species. FHV, functional hypervolume; MPFD, mean pairwise functional distance; MNND, mean nearest neighbor distance; VPFD, variance in pairwise functional distance; VNND, variance in nearest neighbor distance

| Measurements from fish and derivation of traits
Stereo-BRUVs are typically used to estimate relative abundance and biomass of fishes , but can also be used to make fine-scale point-to-point measurements of distances between specified morphological features (e.g., Harvey et al., 2003). To ensure the accuracy and precision of measurements, stereo-BRUVs were calibrated before and after deployments (Boutros et al., 2015).
Measurements were made using EventMeasure software (www. seagis.com.au) with the following rules: Individuals must be within 7 m of stereo-BRUVs, have a maximum root-mean-square error of 20 mm (a quality measure indicating accuracy), and have a maximum precision-to-length ratio of 10% (Harvey & Shortis, 1998).
We identified 15 morphological measurements on the basic body plan for fishes to calculate functional traits pertaining to locomotion and food acquisition (Supplement 1, Figure S1 From the raw morphological measurements, we derived eight functional traits (Table 1). Locomotion was captured by traits documenting eye position, pectoral fin position, caudal peduncle throttling, and elongation, while food acquisition traits included eye size, oral gape position, and jaw length/head length. Body size (considered to be a universal trait, Bellwood et al. (2019)) was measured as total length.

| Data imputation
Complete morphological measurement data were not available for all of the species observed on the stereo-BRUV footage. To fill in missing trait information, we followed a two-step procedure. First, to guarantee an estimate of every trait for every species seen on the stereo-BRUV footage, we obtained raw morphological measurements directly from n = 291 well-preserved museum specimens (belonging to 144 species) from the National Fish Collection at the Museum of New Zealand Te Papa Tongarewa, Wellington (see Myers et al. (2019) and Myers et al. (in press) for details, including species list and voucher registrations). Next, for individuals missing 3 or fewer morphological measurements (out of 15), and for which there were 5 or more occurrences of that species with complete trait information across the dataset (n = 877; 586 rows of data from the stereo-BRUV footage, and 291 from the museum dataset), we imputed missing values. This was done using a random-forest machine-learning model built from all conspecific individuals across the full study design that had the full set of 15 morphological measurements. The random-forest approach handles complex nonlinear relationships, is computationally fast, and estimates imputation error (Penone et al., 2014;Stekhoven & Bühlmann, 2011). We used the "missForest" R package (Stekhoven & Bühlmann, 2011) and set the maxiter argument to 20 iterations and the ntree to 300 decision trees. We performed 20 runs and used the averaged output to obtain a single imputed value. This data imputation increased our stereo-BRUV trait dataset by 136 to yield a total of 722 individuals that had a complete set of 15 morphological measurements for subsequent analysis.

| Choice of individuals to use to represent individual species' traits using a hierarchical approach
We analyzed data on the basis of the species present (observed in video footage) within each of the 47 depth-by-location cells. There were 144 species recorded and 509 species-by-cell occurrences (many species naturally occurred in more than one cell). Our original dataset was comprised of a complete set of 15 raw morphological measurements for a total of 722 individuals (136 of these required some random-forest imputation, and missing traits were remeasured for four individuals) obtained directly from Stereo-BRUV footage.
From this original dataset, two datasets were generated: one for all species-level functional metrics (i.e., FHV, MPFD, MNND, VPFD, and VNND; see Section 2.5 below) and one for all metrics focusing on intraspecific trait variability (i.e., MPFD.I and Prop.I; see Section 2.6 below). To generate the first dataset, we created 100 tables of 509 unique species-by-cell occurrences (rows) for the 8 traits (columns) by randomly drawing 1 individual from the list of all complete individuals for each species. To maintain any spatial structures in trait variability as well as possible, we drew an individual for each species within each cell from conspecific individuals that were (in order of preference): (a) within that depth-by-location cell, (b) at the same depth, or (c) from anywhere within the Stereo-BRUV study design (n = 722) or (d) from a museum specimen (n = 291). All species-level functional metrics were calculated for each replicate table, and we calculated the mean across all 100 tables for every metric for subsequent analyses.
For metrics focusing on intraspecific trait variability, we used data solely from the in situ stereo-BRUV footage (i.e., the dataset comprising 722 individuals). Due to the inability to measure every species observed on the stereo-BRUVs, this dataset represents a reduced number of species (62 out of 144) and cells (43 out of 47). Within this dataset, we were able to measure intraspecific trait variability for 42 species (2 or more individuals per species per cell).
MPFD is the functional analogue to average taxonomic distinctness (Clarke & Warwick, 1998) and mean phylogenetic pairwise distance (Swenson, 2014). It is highly correlated with functional dispersion (Laliberté & Legendre, 2010) and Rao's quadratic entropy (Botta-Dukát, 2005), but is independent of species richness and only weakly influenced by outliers. MNND has been used previously in both phylogenetic (Webb et al., 2002) and functional contexts (Pigot et al., 2016;Swenson & Weiser, 2014). It measures the average minimum distance in the functional strategies of co-occurring species and has previously been used to estimate functional originality (Leitao et al., 2016;Mouillot, Graham et al., 2013).
VPFD quantifies the regularity of the distances among species in the functional space. This measure was originally proposed by Clarke and Warwick (2001) in a taxonomic setting. Here, we calculate it in functional space (as suggested by Somerfield et al., 2008), and it is independent of species richness and MPFD. Similarly, VNND also quantifies regularity, but focuses on functional similarity between nearest neighbors. It was proposed by Swenson (2014) in a phylogenetic context, but we calculated it here in functional space.
We also performed principal component analysis (PCA) on the normalized traits in order to calculate functional hypervolume (FHV; Blonder et al., 2014Blonder et al., , 2018. FHV was calculated using the first 4 principal component axes (which accounted for 70.2%-74.4% of the total variation in the 8D functional trait space across the 100 species-bytrait tables). We did not retain all 8 dimensions due to difficulties associated with calculating FHV when few species were present. FHV has been used as a proxy to estimate niche space, including high-dimensional, irregular spaces (Cooke et al., 2019;Lamanna et al., 2014). The Gaussian kernel density estimation method was chosen to form a "loose wrap" around the data (Blonder et al., 2018).
We used a quantile threshold of 0.05, retaining 95% of the total probability density, and estimated the bandwidth vector using the Silverman method (Blonder et al., 2018).

| Quantifying intraspecific trait variability
Next, we calculated individual-level trait variation and determined its contribution toward total functional trait variation  Figure S1 for an illustration) where 2 I is the estimated individual-level trait variance, 2 S is the estimated species-level trait variance, MS Res is the residual mean square, MS Species is the "Species" factor's mean square from the PERMANOVA partitioning, and n 0 is a divisor that allows for unequal numbers of individuals (n i per species I = 1,…S) within each cell, calculated as follows (Sokal and Rohlf (1981), p. 214): The proportion of total trait variation (within each cell) attributable to individual-level variation (i.e., Prop.I) is then calculated as. This is the first time, to our knowledge, that individual-level functional variation and species-level functional variation have been partitioned and calculated in this way for full-dimensional multivariate trait space. Importantly, the estimators 2 I and 2 S are independent of species richness. Also, while small or unbalanced sample sizes (n i , the number of individuals per species) will necessarily affect the precision of the estimates (i.e., small numbers of individuals or species will yield more variable estimates), they will not, however, affect the accuracy of these measures, which are unbiased. There was very Some previous approaches to disentangle intra-from interspecific trait variation for community-level data (e.g., Lepš et al., 2011) have focused simply on differences between two weighted aver- The method described by de Bello et al. (2011) is aimed at quantifying functional turnover among habitats primarily for singular traits, rather than quantifying functional diversity within habitats and across multiple traits, which is our primary focus here. These authors do usefully mention the PERMANOVA approach and note its equivalences, where appropriate, with partitioning according to quadratic entropy diversity measures. They do not articulate or distinguish, however, the crucial additional calculations required to transform the raw additive sums of squares (arising directly from either a PERMANOVA or ANOVA partitioning) into estimates of independent variance components. The relative sizes of sums of squares from a partitioning will depend on their degrees of freedom, whereas variance components are unbiased estimators of variance derived from expectations of mean squares. Such components, calculated here in a functional context, can be compared across different factors (e.g., in hierarchical/nested designs; see Anderson et al., 2005) and are independent of sample size (numbers of species or individuals).

| Univariate models
For modeling, we considered depth and latitude as continuous predictor variables and included quadratic and cubic terms to allow for nonlinearities in response variables along these gradients. We chose to use a parametric model rather than a generalized additive model (GAM), because there were only 7 depth strata and 7 latitudinal bands. We consider that GAMs are better suited for data where there are more continuous x-values along the gradient (independent variable) of interest. We normalized depth and latitude to ensure polynomial terms remained orthogonal. We analyzed all metrics described above as response variables; VPFD and VNND were transformed to y � = log 10 () prior to analysis to improve normality. In order to explore the strength and generality of the potential relationship between the functional diversity metrics and either depth or latitude, we used linear mixed-effects models, first treating depth as a fixed factor and latitude as a random factor, then treating latitude as fixed and depth as random (see Supplement 2, Tables S1-S2 for further details). The mixed model allows a test and quantification of the effects (if any) of a chosen fixed gradient of interest, over and above any potential variation in those effects (i.e., possible interaction) across a second (random) factor (e.g., as in Quintero & Jetz, 2018). In addition, we have considered depth and latitude as fixed factors and tested their interaction. We used a linear model to test this interaction when only the linear term for both depth and latitude was retained by the selection of the best linear mixed model. If at least one polynomial term was retained for depth or latitude, we used a GAM to assess the significance of a tensor product smooth term to test the nonlinearities in the interaction between depth and latitude. Despite GAM being better suited for a more continuous independent variable, when testing the nonlinearities in the interaction between depth and latitude, we have now included 47 depthby-latitude combinations which is more in line with the use of GAM, rather than investigating the depth or latitude pattern alone considering only seven modalities. We constrained the smoothing terms for a maximum of 3 degrees of freedom to avoid overfitting. GAMs were performed with the mgcv R package (Wood, 2004).

| Multivariate analyses
To visualize changes in multiple functional metrics along depth and latitude gradients simultaneously, we did a metric multidimensional scaling (mMDS) ordination on (a) depth centroids and (b) latitude centroids (see Figure S4 for an mMDS ordination of all 47 depth-bylocation cells). We superimposed bubbles corresponding to species richness and vectors to show partial correlations of functional metrics with mMDS axes. Three-dimensional shade plots were obtained to visualize potential interactions between depth and latitude for each of FHV, MPFD, MNND, MPFD.I, and Prop.I ( Figure S2). All multivariate analyses were done on the basis of Euclidean distances for p = 7 functional metric variables using PRIMER v7 (Clarke & Gorley, 2015) with the PERMANOVA+add-on (Anderson et al., 2008).

| RE SULTS
Functional hypervolume (FHV) was stable across the depth gradient, which was in line with our second prediction derived from our conceptual model (Figure 1b). The relationship between FHV and depth was not statistically significant (p =.36); however, FHV tended to generally increase with increasing depth, with the greatest volume found at the deepest depth (1200 m) (Figure 2a) Table S1, see also Table S2). The interaction between depth and latitude was significant and positive for FHV, MPFD, VPFD, and negative for Prop.I (Table S3). Despite tensor products being nonsignificant for MNND and VNND a certain degree of interaction could be detected visually (see Figure S5 and Table S4).
There were clear gradual changes in the functional diversity of fishes along the depth gradient from shallow to deep environments (i.e., from left to right along MDS axis 1 in Figure 4a). This was primarily characterized by increases in MNND and MPFD.I at deeper depths ( Figure 4a). There were also clear sequential latitudinal changes in functional diversity metrics from north to south (i.e., from left to right along MDS axis 1 in Figure 4b Figure S4).

| D ISCUSS I ON
This study examined several functional diversity metrics simultaneously, and also partitioned variation in functional space into intra- In addition, we found that functional hypervolume, functional dispersion, and functional regularity decreased with increasing latitude. These results are consistent with our first conceptual model (Figure 1a), in which environmental filtering is the dominant process. Interestingly, although species-level metrics decreased with increasing latitude, individual-level metrics increased (Figure 3), suggesting that there was a gradual change in the source of functional diversity, with intraspecific trait variability becoming increasingly important at higher latitudes.
The two gradients of depth and latitude showed contrasting patterns with respect to functional hypervolume (Figures 2-3), which remained stable with increasing depth (the positive trend was not significant), but decreased with increasing latitude. The stability of the functional hypervolume metric across the depth gradient is a striking result because it provides a rare example of a mis-match between species richness and functional richness; generally, when the number of species in a community increases, it is more likely to lead to a greater functional volume of that community (Villéger et al., 2008). The weak positive trend with depth, albeit nonsignificant, was an unexpected result and is certainly inconsistent with the idea that functions will be filtered more strongly in harsh environments (Swenson, 2011). Instead, abiotic conditions such as limited trophic resources and habitat availability, decreasing temperature, and increasing pressure may represent key selection pressures on individuals living in the deep sea (Ramirez-Llodra et al., 2010). Species living in extreme conditions may be subject to greater disruptive selection and/or character displacement, potentially contributing to distinct morphologies, trait combinations, or functional strategies (Weiher & Keddy, 1995) that enable a greater variety of unique biotic adaptations for resource acquisition (Leitao et al., 2016). This may allow greater partitioning of limited resources, consistent with the limiting similarity hypothesis (MacArthur & Levins, 1967).
Morphological dissimilarities among deep-sea species reflect low niche overlap, which can promote coexistence in a low-resource environment (Kumar et al., 2017). Decreasing functional hypervolume with increasing latitude, however, follows a more classical stress-gradient hypothesis, whereby traits are filtered more strongly in harsh than in benign environments (Swenson, 2011;Weiher & Keddy, 1995). This pattern has been documented for plants versus latitude , for birds versus altitude (Pigot et al., 2016), and for macroinvertebrate assemblages versus depth (Ashford et al., 2018).
The packing, or density, of species within functional space decreased with increasing depth, as the mean distance between species and mean nearest neighbor distances both increased ( Figure 2a). Previous work has shown that variance in pairwise F I G U R E 2 Relationships between functional diversity metrics (mean ± 1SE across 7 locations) and depth: (a) specieslevel metrics (S = species richness) and (b) intraspecific trait variability metrics. Lines correspond to fitted values for the best model when depth was considered a fixed factor and latitude a random factor (see Table S1 in Supplementary Material). Solid lines show statistically significant trends (p < .05); dashed lines show trends that did not reach statistical significance at the 0.05-level. The conditional R 2 (Condi.R 2 ) and marginal R 2 (Margi.R 2 ) are overlaid for each metric and represent variation explained by both the fixed and random effects jointly and variation explained only by the fixed effects, respectively (Nakagawa et al., 2017). Blue dots represent the average value per depth. Black horizontal bars, black dots, and boxes show the median, outliers, and interquartile range, respectively. Whiskers represent 1.5 times the interquartile range phylogenetic distance (Eme et al., 2020) (Swenson & Weiser, 2014). In traitbased ecology, studies advocating measurement of intraspecific trait variability are becoming increasingly common (Des Roches et al., 2018;Siefert et al., 2015;Violle et al., 2012). This is the first time, to our knowledge, that multivariate variation across multiple functional traits has been partitioned into intra-and interspecific components using a hierarchical PERMANOVA approach (but see However, the proportion of total functional trait variation attributable to individual-level variability increased with increasing latitude (Prop.I; Figure 3b), supporting the idea that intraspecific trait variation becomes important in species-poor communities with narrow environmental breadth (Siefert et al., 2015). For example, the highlatitude subantarctic Auckland Islands (AUC) had the highest value for Prop.I and was characterized by low species diversity, and a narrow temperature ranges from shallow to deeper environments (9.3-5.5°C). High levels of intraspecific phenotypic variability may contribute to the higher rates of speciation in fishes in the deepsea and at high latitudes (Rabosky et al., 2018), boosting the evolutionary potential of populations (Jump et al., 2009). Intraspecific variation begets speciation, which, in turn, begets interspecific variation (Darwin, 2009;Pfennig & Pfennig, 2010).

F I G U R E 3
Relationships between functional diversity metrics (means ± 1SE across 7 depth strata) and latitude in degrees south (KER = Rangitāhua, Kermadec Islands, TKI = Three Kings Islands, GBI = Great Barrier Island, WI = Whakaari, White Island, KKA = Kaikōura, OTA = Otago, and AUC = Auckland Islands): (a) species-level metrics; (S = species richness) and (b) intraspecific trait variability metrics. Lines show fitted values for the best model for each metric when latitude was considered a fixed factor, and depth a random factor (see Table S2 in Supplementary Material). Solid lines show statistically significant trends (p < .05); dashed lines show trends that did not reach statistical significance at the 0.05-level. The conditional R 2 (Condi.R 2 ) and marginal R 2 (Margi.R 2 ) are overlaid for each metric and represent variation explained by both fixed and random effects jointly and variation explained only by the fixed effects, respectively (Nakagawa et al., 2017). Blue dots represent the average value per latitude. Black horizontal bars, black dots, and boxes show the median, outliers, and interquartile range, respectively. Whiskers represent 1.5 times the interquartile range Despite intraspecific variation being at the core of evolutionary biology, community ecology has focused mostly on interspecific variations due to the lack of data on intraspecific variation for large communities. However, greater connections between evolutionary and ecological studies (Johnson & Stinchcombe, 2007;Mouquet et al., 2012), and the emergence of larger datasets aid the revival of studying intraspecific variation within community ecology (Siefert et al., 2015). Our intraspecific results were based on a limited number of individuals per species per site (4.32 on average) that should, ideally, be improved. However, measuring more individuals per species per site is not always possible depending on the species' behavior (e.g., territorial or long-ranging species). These findings, nevertheless, illustrate the importance of considering intraspecific variation in trait-based studies at the community level.
Overall, our results, supporting inter-and intraspecific competition as a potential driver of niche partitioning, question the primary role of abiotic filtering in these harsh environments (Priede, 2017).
Competition may have a greater role in shaping communities of the deep sea than previously thought. Our study has provided novel insights into how functional diversity changes along environmental gradients at the local (alpha) scale. We consider that future work should examine turnover and nestedness components of functional beta diversity to yield further potential insights into how ecological processes may structure communities. We also consider that functional traits from undersampled taxa and environments need additional study (Borgy et al., 2017). In addition, multivariate analyses of morphological traits should be extended to include behavioral traits, life-history strategies, trophic positions, and/or physiological traits, for a more holistic measure of biologically relevant trait space (Bellwood et al., 2019;Violle et al., 2007).

| CON CLUS IONS
These results suggest that interspecific and intraspecific competitions act as key processes shaping the functional diversity of fishes in the deep sea. Increasing morphological dissimilarity with increasing depth may help to facilitate niche partitioning and promote coexistence, whereas external abiotic filtering may be the dominant factor structuring communities with increasing latitude.
In an era characterized by rapid and unprecedented change to deep-sea environments, with increasing anthropogenic pressures from fishing, deep-sea mining, and global climate change (Levin & Le Bris, 2015;Levin et al., 2016;Watson & Morato, 2013), understanding how functional diversity changes along large spatial gradients may help to predict potential responses of ecological communities to disturbances. In summary, this study quantified trait variation in marine fishes across broad-scale depth and latitudinal gradients, shedding new light on the potential roles of abiotic filtering, biotic interactions, and niche partitioning, to further our understanding of the mechanisms underlying large-scale patterns in biodiversity.

ACK N OWLED G EM ENTS
We thank James Seager from SeaGIS for developing the head morphometrics feature of EventMeasure specifically for this project. F I G U R E 4 Metric multidimensional scaling (mMDS) ordination of normalized functional diversity metrics on the basis of Euclidean distances among (a) depth centroids (50 m to 1,200 m) and (b) location centroids. Overlaid arrows follow a shallow to deep, and north to south trajectory, respectively. Bubble sizes are proportional to mean species richness (also provided as a value inside each bubble). Vectors (right) show multiple partial correlations for each of p = 7 functional metrics with the mMDS axes

CO N FLI C T O F I NTE R E S T
The authors declare the absence of conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The dataset is available at the Dryad Digital Repository with DOI accession number: https://doi.org/10.5061/dryad.xgxd2 54gt