Trait similarity in reef fish faunas across the world’s oceans

To quantify the distribution of reef fish traits, we first created a multidimensional trait space of 2,786 reef fish species (using principal coordinates analysis, PCoA) according to five categorical traits. Trait variation among species was primarily explained by water-column position (benthic, demersal, site-attached pelagic, or mobile pelagic), gregariousness (solitary, paring, or schooling), and active period (diurnal or nocturnal). Along the first axis, there was a clear gradient from benthic to pelagic species, and solitary and nocturnal species were distinguished from schooling and diurnal species (Fig. 1 A and C). Species also varied largely in diet, with invertivores being distinct from planktivores and piscivores along the first axis. Along the second axis, species primarily varied in body size, with a clear gradient from small to large species.

Across five tropical and three temperate marine realms, we found a common pattern in the distribution of species in trait space despite a 10-fold gradient in species richness (Fig. 1D). In all realms, species clustered in the same area of trait space and occupied similar ranges of the overall space (Fig. 1D). The proportions of different trait categories were also strikingly similar across realms (SI Appendix, Fig. S1). Thus, reef fish trait composition appears remarkably consistent worldwide when considered at the scale of marine realms, despite major differences in species number and identity. Across realms, species richness ranged from 6 to 54% of the observed global pool, whereas the number of functional entities (i.e., unique trait combinations) ranged from 24 to 77%, indicating that species-poor realms maintain disproportionately high trait diversity (SI Appendix, Table S1). Mouillot et al. (19) previously found that the breadth of reef fish traits was consistent across the tropics, with a ninefold gradient in species richness corresponding to only a threefold gradient in the number of functional entities. Here, we extend those patterns to temperate oceans, showing astonishing similarity between reef fish faunas worldwide despite drastic differences in species richness, environmental regimes, and habitat types (i.e., coral versus rocky substrates). Together, these findings show that key ecological functions are retained, even in species-poor realms, through high differentiation and limited redundancy that capture the overall breadth of reef fish traits (19).

We next identified functional entities shared by all marine realms, both temperate and tropical (Fig. 1B). A total of 21 entities (of 356 in the global species pool) define a common “backbone” (20) in trait space filled by reef fishes from all ocean basins. These 21 core entities encompass a broad range of trait values, including a variety of body sizes, behaviors, and diets (SI Appendix, Fig. S2). In particular, the most common backbone traits were invertivorous diet, demersal water-column position, diurnal active period, solitary behavior, and small to medium body size, most frequently represented by families such as wrasses, sea basses, groupers, damselfishes, and gobies. We also identified specific backbones for temperate and tropical realms, which revealed consistent ecological structuring within temperate (42 entities) and tropical (45 entities) oceans (Fig. 1B and SI Appendix, Fig. S2). However, null models revealed that the numbers of shared functional entities within the global and tropical backbones were lower than expected by random (SI Appendix, Fig. S3). This suggests that similar trait distributions across realms were not due to identical functional entities but rather similar trait proportions. Hence, similar overall trait composition can arise from a variety of functional entities, depending on their relative frequencies, particularly when functional entities share multiple traits in common. Furthermore, this result indicates that intense filtering by environmental and ecological constraints has limited the number of viable functional entities that can persist in all ocean basins.

To evaluate similarity in species assemblage trait composition at the regional scale, we examined trait space patterns across 89 ecoregions globally. We generated probability densities of PCoA 1 and PCoA 2 scores (i.e., species distributions on the first and second axes of trait space) for regional species pools and calculated the centroid position (i.e., mean value along each axis, SI Appendix, Fig. S4) for each ecoregion, as well as species richness and the number of functional entities. Across ecoregions, we found that PCoA 1 centroids had a clear unimodal relationship with latitude, being similar in temperate regions in both hemispheres and converging toward the equator (Fig. 2), highlighting similar trait compositions at similar latitudes worldwide. Linear models of PCoA 1 centroids using the absolute value of latitude revealed that including hemisphere as a predictor added no explanatory power (latitude: R² = 0.327; latitude + hemisphere: R² = 0.332), further demonstrating latitudinal similarity regardless of hemisphere. Reef fish assemblages in temperate ecoregions were distinguished from those in tropical ecoregions along the first axis of trait space (Figs. 1 and 2), with temperate ecoregions more characterized by benthic, solitary, nocturnal, and invertivore species and tropical ecoregions more characterized by demersal, pairing, diurnal, and grazer species (SI Appendix, Fig. S5 and Table S2). PCoA 2 centroids also showed a unimodal relationship with latitude, indicating that body sizes were larger at higher latitudes and that temperate regions were more characterized by piscivores (Fig. 2C and SI Appendix, Fig. S5).

We further found that both species richness and the number of functional entities were lowest at temperate latitudes and increased toward the equator; however, species richness increased 2.5 times faster than the number of functional entities (2.8 species versus 1.1 functional entities per degree latitude; SI Appendix, Fig. S6). Thus, as species richness increased, an average of 2.5 species were packed into existing functional entities before a new functional entity emerged. Together, these results indicate that as species richness increases in the tropics, species are packed into existing functional entities, increasing redundancy rather than extending the breadth of trait space (19, 21).

To determine if trait similarity resulted from shared species composition (taxonomy) or evolutionary history (phylogeny), we examined dissimilarity patterns between realms and between ecoregions. We defined trait dissimilarity as the Euclidean distance between species-pool centroids in the first four dimensions of trait space and defined taxonomic and phylogenetic dissimilarity as the turnover component of pairwise Jaccard dissimilarity (22). While both taxonomic and phylogenetic dissimilarity are naturally bounded between 0 and 1, trait dissimilarity (distance between ecoregion centroids) has no intrinsic upper limit. We therefore used a null model approach to define the maximum expected trait dissimilarity between realms under random assembly (see Materials and Methods). For consistency, this approach was also used for both taxonomic and phylogenetic dissimilarity. We found that marine realms were highly similar in trait composition but highly dissimilar in taxonomic and phylogenetic composition, suggesting that common trait distributions are not a result of shared species or evolutionary histories (Fig. 3). We also found that when ecoregions became farther apart in latitude, they became more dissimilar in taxonomic and phylogenetic composition. However, trait composition was not related to latitudinal distance since both geographically near and far ecoregions had similar trait compositions. Thus, ecoregions separated by up to 100 degrees of latitude displayed equivalent trait compositions despite hosting radically different species pools from separate evolutionary lineages (Fig. 4 and SI Appendix, Fig. S7).

Debate over the relative influences of contemporary environmental conditions and shared evolutionary history has focused extensively on species composition rather than trait composition (13, 16). Therefore, to determine if trait similarity was more strongly driven by environmental conditions than by taxonomic or phylogenetic composition, we evaluated the relationship between trait dissimilarity and taxonomic, phylogenetic, and environmental dissimilarity using multiple regression of distance matrices (MRM) (23). We defined environmental dissimilarity as the Euclidean distance between ecoregions according to mean sea-surface temperature (SST), SST range, nitrate, phosphate, and net primary productivity (see Materials and Methods). We found that environmental dissimilarity was the primary driver of trait dissimilarity (effect size = 0.60 ± 0.01 and R² = 0.56), while taxonomic (effect size = 0.09 ± 0.02 and R² = 0.20) and phylogenetic dissimilarity (effect size = 0.11 ± 0.02 and R² = 0.24) added little explanatory power (all variables together R² = 0.58; Fig. 4).

Trait similarity among geographically distant assemblages inhabiting similar environmental conditions could arise through both environmental filtering—in which regional trait composition is filtered from the realm pool—and evolutionary convergence—in which species with distinct trait combinations independently evolve in response to similar environmental pressures (24, 25). Both processes have likely contributed to contemporary patterns in trait composition yet disentangling their relative effects is challenging and contentious (25). Therefore, to assess the likelihood of convergent evolution in shaping spatial trait similarity among assemblages, we tested whether the observed relationship between trait dissimilarity and taxonomic, phylogenetic, and environmental dissimilarity was different from expected under neutral evolution by simulating trait dissimilarity according to a Brownian model of evolution. Following Mazel et al. (25), we simulated species’ positions on the first axis of trait space, recalculated trait dissimilarity between ecoregions, recalculated the relationship between trait dissimilarity and taxonomic, phylogenetic, and environmental dissimilarity, and compared null expectations to the observed relationship. We also calculated the magnitude of difference in observed versus expected trait dissimilarity between ecoregions [i.e., standardized effect size (SES), termed “phylogenetically standardized trait distance” in Mazel et al. (25)]. We found a positive, albeit weak, relationship (slope = 0.036; SI Appendix, Fig. S8) between SES and environmental dissimilarity, with 9% of ecoregion pairs sharing more similar trait compositions than expected under neutral evolution. More importantly, we found that the observed relationship between trait dissimilarity and environmental dissimilarity was stronger (R² = 0.45) than expected under neutral evolution (mean of R² null distribution = 0.21, 95% = 0.43; SI Appendix, Fig. S9), whereas the relationship between trait dissimilarity and both taxonomic and phylogenetic dissimilarity did not differ from neutral expectation (observed taxonomic dissimilarity: R² = 0.06, mean of R² null distribution = 0.08, 95% = 0.13; observed phylogenetic dissimilarity: R² = 0.08, mean of R² null distribution = 0.12, 95% = 0.21; SI Appendix, Fig. S9).

Fish data were collected through the RLS, a standardized monitoring program that surveys fish communities on coral and rocky reefs worldwide (www.reeflifesurvey.com). RLS was chosen for this study because 1) all surveys follow a standardized protocol, and patterns between different biogeographic regions are readily comparable; 2) RLS is conducted uniquely on shallow reef habitats; and 3) RLS has extensive spatial coverage, including both temperate and tropical regions. RLS is comprised of teams of professional scientists and trained recreational divers who complete underwater visual censes of fish communities using 50 × 10 m belt transects. On each transect, divers record the species identity, abundance, and estimated size of all fishes encountered within 5 m of either side of the transect. For full details refer to Edgar and Stuart-Smith (17). For this study, fish data were compiled from 10,913 transects (conducted between 2006 and 2018), encompassing 2,786 taxa from 89 marine ecoregions spanning 10 marine realms (two realms were not included in realm-level analyses, see Traits, Trait Space, and Common Trait Backbones). Of the 2,786 observed taxa, 2,582 were recorded at species level, with the remaining recorded primarily at genus level (162 at genus level and 43 at family level); however, for simplicity we refer to all taxa as species throughout the paper. Ecoregions and realms were defined according to the Marine Ecoregions of the World categorization (18), which was also used to classify ecoregions as temperate or tropical. For this study, we considered species occurrences and did not integrate abundance or biomass. Species occurrences were chosen to examine trait composition and redundancy within different realms and biogeographic regions and to identify how these patterns are shaped by environmental conditions and evolutionary history. Species occurrences reveal changes in trait composition due to differences in the existence or redundancy (number of species) of functional entities (i.e., unique trait combinations), whereas species abundances can lead to drastically different trait compositions even when species pools have identical trait spaces. For instance, substantial variation exists in species evenness between temperate and tropical realms, which has strong implications for trait distribution among individuals (39). All indices and metrics throughout the manuscript were therefore calculated using presence/absence data; however, sensitivity tests were performed using species abundances (see Sensitivity of Results to Occurrence versus Abundance Data).

We compiled five categorical traits for each of the 2,786 species in the data set, which were chosen to describe life history, trophic ecology, habitat preference, and behavior. While we recognize that these traits only represent a small portion of the true multidimensional species niche, they are recognized proxies of fish functions (5) and have been used to examine reef fish biogeographic patterns (19, 38, 39) and impacts on ecosystem functioning (40). Diet characterized species’ main food items and was coded using seven trophic categories: grazer, microphage (feeding on cyanobacteria and microorganisms), planktivore, omnivore, corallivore, invertivore, and piscivore. Microphage was chosen in lieu of “scraping herbivore” or “excavator” based on recent reviews by Clements et al. (41) and Bellwood et al. (42). Water-column position was coded using four ordered categories: benthic (sedentary on the bottom), demersal (swimming near the bottom), pelagic site attached (swimming off the bottom within a reef), and pelagic mobile (swimming off the bottom between reefs). Active period was coded as either diurnal (active during the day) or nocturnal (active during the night). Gregariousness was coded using three ordered categories: solitary, pairing (or sometimes forming small schools), and schooling. Body size was coded using six ordered categories: 0 to 7 cm, 7.1 to 15 cm, 15.1 to 30 cm, 30.1 to 50 cm, 50.1 to 80 cm, and >80 cm. Trait data came from Stuart-Smith et al. (39), the majority of which were derived from FishBase and expert knowledge. The multidimensional trait space was created by applying PCoA to a Gower similarity matrix of the species × traits table and was primarily examined by plotting the first two Principal Coordinate axes (43), which cumulatively explained 39% of the total variance. Trait space was additionally examined using the third and fourth Principal Coordinates axes which together explained 25% of the total variance (SI Appendix, Figs. S10 and S11). Thus, the first four axes of trait space accounted for 64% of overall variation, and because a test of mean absolute deviation (44) identified four axes as the optimal number for describing trait space, additional axes were not examined.

While we analyzed trait biogeography by considering single trait values for each species, it should be noted that trait–environment relationships operate at the scale of individuals, and substantial trait variation may exist within species (45). Traits may also be highly plastic and can evolve in response to environmental or ecological pressures (46). In turn, this intraspecific variability can influence ecosystem functioning (47). Here, we were unable to integrate intraspecific trait variability due to insufficient data; however, by using broad categorical traits (e.g., diet, water-column position, and size bin) rather than highly specific continuous traits (e.g., length at maturity and fecundity), we capture the ecological profile of adult individuals.

To evaluate reef fish trait composition globally, we began by examining the distribution of species within trait space for three temperate and five tropical marine realms. Temperate South America and Temperate Southern Africa were not included in realm-level analyses due to low and patchy sampling effort, as few sites and ecoregions have been sampled within these realms. For instance, only 25 species have been observed in Temperate Southern Africa in comparison to 154 species in Temperate North Atlantic, the most species-poor realm included. Thus, we were not able to confidently assess trait-space patterns for Temperate South America and Temperate Southern Africa. However, the individual ecoregions sampled within these realms were included in ecoregion-level analyses, which were separate from realm-level analyses. We calculated 50, 75, and 95% density contours for each realm based on the position and density of species in trait space using the R package emdbook (48) function HPDregionplot. We then generated probability density plots of PCoA 1 and PCoA 2 scores for all species in each realm to further examine the distribution of species along each of the trait-space axes. To examine which traits were most responsible for trait-space structure and to examine how species were distributed according to their trait combinations, we plotted convex hulls of each trait category within the trait space.

Functional entities refer to unique combinations of trait categories, that is, the number of functional entities for a given species pool is the total number of combinations of the different trait categories present (19). Species in the same functional entity also occupy identical positions in trait space. Here, the 2,786 species and five traits considered led to 1,008 possible trait combinations, of which 356 were realized in the global species pool. Common trait “backbones” were calculated by identifying functional entities that were shared by assemblages in 1) all realms, 2) temperate realms only, and 3) tropical realms only. Temperate Australasia is considered a temperate realm but borders the Great Barrier Reef along the east coast of Australia and has been extensively sampled by the RLS program, leading to high species richness and similarities to Indo-Pacific realms. We therefore examined how temperate and tropical backbones changed if we reclassified Temperate Australasia as a tropical realm or if we removed it entirely. In both cases, there was little change, as reclassifying Temperate Australasia as tropical reduced the tropical backbone from 45 to 43 functional entities and removing it from temperate realms expanded the temperate backbone from 42 to 44 entities.

Several null models were used throughout this study to assess the sensitivity of results and to identify whether results could have arisen because of chance (random assembly) or from confounding patterns like latitudinal gradients in species richness. These null models all use the approach of randomizing species composition across assemblages (realms or ecoregions) while maintaining observed species richness and species occurrence frequency, using the “independent swap” algorithm (49) in the R package picante (50). We chose to randomize species composition among assemblages rather than randomize species-trait values to avoid creating species with ecologically unfeasible trait combinations.

To determine whether the number of shared functional entities comprising global, temperate, and tropical backbones was different from expected by chance, we randomly shuffled species composition between realms while maintaining initial species richness and species occurrence frequency and calculated the number of shared functional entities. We repeated this process 999 times to build null distributions of shared functional entities and compared the actual observed numbers of functional entities to the null distributions. Both the global and tropical backbones contained fewer shared functional entities (global = 21 and tropical = 45) than 95% of the null distribution (global null 5% = 26 and tropical null 5% = 50), while the temperate backbone contained fewer entities than 90% of the null distribution but not fewer than 95% (temperate = 42 and temperate null 5% = 40; SI Appendix, Fig. S3).

PCoA centroids refer to the central positions of species pools within trait space and reflect the mean trait values of the observed species. PCoA centroids were therefore calculated for each ecoregion as the mean positions of observed species pools on each of the PCoA axes (SI Appendix, Fig. S4). To calculate PCoA centroids and functional entities for each ecoregion while accounting for uneven sampling effort across regions, we randomly selected six transects per ecoregion to generate species pools, calculated metrics for the resulting species pools, repeated this process 99 times with replacement, and calculated mean values for the 99 repetitions. For PCoA centroids, in each of the 99 repetitions, the ecoregion centroid value was calculated as the mean centroid of the six transects. We chose six transects because it was the most common number of transects performed within any ecoregion, that is, the mode. The mode was chosen in lieu of mean or median because the distribution of sampling effort (number of transects per ecoregion) was highly skewed, with few ecoregions containing disproportionate numbers of transects (SI Appendix, Table S3). Transects were randomly sampled with replacement because some ecoregions had fewer than six transects; however, we tested the robustness of these choices (see below).

To help visualize latitudinal similarity in ecoregion centroids, we used second-order polynomial linear regressions. However, the purpose of these regression models was primarily visualization—providing a trend line to highlight the unimodal latitudinal pattern in ecoregion centroids (and additionally for a null model, see Null Model for the Influence of Species Richness). Yet, we do note that the relationship between PCoA 1 centroids and latitude was substantially better fit by a second order polynomial regression (R² = 0.46 and Akaike Information Criterion [AIC] = −254.8) than by a simple linear regression (R² = 0.007 and AIC = −203.4). To provide quantitative evidence for latitudinal similarity in trait composition, we also used simple linear regressions of PCoA 1 centroids in functional of the absolute value of latitude and examined whether including hemisphere as a predictor improved explanatory power.

To test the robustness of the number of transects and procedure used to calculate ecoregion centroids, we recalculated PCoA 1 centroids and reexamined latitudinal similarity patterns by randomly selecting 1, 5, 10, 15, and 20 transects per ecoregion instead of 6. This was done first with replacement and no deletion, second with replacement but after deleting all ecoregions with less than the selected number of transects, and third without replacement after deleting all ecoregions with less than the selected number of transects. Latitudinal trait similarity was robust to all choices except deleting ecoregions with fewer than 20 transects, because this removed over 40% (38 of 89) of ecoregions (SI Appendix, Fig. S12).

As visible in Fig. 2A, distances between ecoregion centroids were much smaller than the extent of the overall trait space, implying that variation across ecoregions was limited. However, the extent of the trait space is defined by the maximum distance between the most dissimilar individual functional entities, which is an inappropriate context for evaluating variation among regional species assemblages. To determine the expected maximum distance between ecoregion centroids, we randomly shuffled species composition among ecoregions while maintaining initial species richness and species occurrence frequency and calculated the maximum distance between ecoregions in the first four dimensions of trait space. We repeated this process 999 times to build a null distribution of expected maximum distances between ecoregion centroids. The maximum null value was 0.60, substantially lower than the maximum distance between any two functional entities (1.24). The actual observed value of maximum distance between ecoregion centroids (0.550) was greater than 95% of the null distribution (95 null = 0.549), demonstrating that ecoregions displayed substantial variation in trait composition that was greater than expected under random assembly (SI Appendix, Fig. S13).

To examine the contribution of each trait to latitudinal patterns in PCoA 1 centroids, and to examine the distinction between temperate and tropical ecoregions, we calculated average trait values per ecoregion as the proportions of different trait categories present in each species pool using the R package FD (51) function functcomp. We then calculated Pearson correlations between PCoA centroid values and trait values for each ecoregion, and we used multiple factor analysis (MFA) to examine associations between traits and ecoregions. MFA is a modified PCoA that can account for categorical data and therefore give equal weighting to traits with different numbers of categories, for example, diet, which has seven categories, is not overweighted. Associations between traits and temperate versus tropical ecoregions were assessed by the contributions of trait values to the first and second axes of the MFA analysis (SI Appendix, Fig. S5). MFA was conducted using the R package FactoMineR (52) function MFA.

In this study, trait dissimilarity was calculated as the Euclidean distance between species-pool centroids in the first four dimensions of trait space, which is not bounded by an upper limit of 1. To define the expected upper limit of trait dissimilarity between marine realms, we randomly shuffled the species composition of each realm while maintaining initial species richness and species occurrence frequency and recalculated the Euclidean distance between realms. We repeated this process 999 times and calculated the maximum trait dissimilarity between realms across the 999 repetitions. The observed maximum value (0.26) was then used to define the upper limit of trait dissimilarity between realms in Fig. 3A. This procedure was chosen to provide a conservative estimate of maximum dissimilarity between realms, rather than defining maximum dissimilarity between two species in the trait space (the theoretical maximum distance between assemblages that would each have only one of these two extreme species). For consistency, this procedure was also applied to taxonomic and phylogenetic dissimilarity, which resulted in a maximum taxonomic dissimilarity of 1 and a maximum phylogenetic dissimilarity of 0.91.

To evaluate taxonomic and phylogenetic dissimilarity, we use the additive partitioning of pairwise Jaccard dissimilarity proposed by Baselga (22) for taxonomy and adapted by Leprieur et al. (53) for phylogeny, which decomposes beta diversity into nestedness and turnover components. Because we were interested in knowing whether similarity in trait composition resulted from similarity in species composition or evolutionary lineages, we focused on taxonomic and phylogenetic turnover—dissimilarity that occurs when species or lineages present in one region are absent in another region but are replaced by other species or lineages not found in the first region independently of the difference in species richness. As with PCoA centroids, we randomly selected six transects per ecoregion to generate species pools, calculated taxonomic and phylogenetic turnover for the resulting species pools, repeated this process 99 times, and calculated mean values for the 99 repetitions. Phylogenetic trees were derived from the R package fishtree (54), and taxonomic and phylogenetic turnover were calculated using the R package betapart (55). For this analysis, two ecoregions with only five species present in the phylogeny (Malvinas/Falklands and Channels/Fjords of Southern Chile) were removed. Environmental dissimilarity between ecoregions was calculated according to SST mean and range, nitrate, phosphate, and net primary productivity. Following extensive discussion among coauthors, these environmental variables were chosen based on existing knowledge of reef fish functional responses to environmental gradients (56). However, we tested the sensitivity of results to the choice of environmental variables, including adding pH and salinity (see Sensitivity of Results to the Choice of Environmental Variables). Data for mean SST and SST range were obtained from the National Oceanic and Atmospheric Administration Coral Reef Watch (57). Data for nitrate and phosphate were obtained from the Bio-ORACLE environmental data set (58). Net primary productivity data were derived from the standard Vertically Generalized Production Model (59). We calculated mean values of environmental variables per ecoregion and standardized them (using the R function scale) before calculating Euclidean distances between all pairs of ecoregions. As for realms, we calculated trait dissimilarity between ecoregions as the Euclidean distance between ecoregion centroids in the first four dimensions of trait space. We then used MRM [R package ecodist (60) function MRM (23)] to quantify the relative influences of taxonomic, phylogenetic, and environmental dissimilarity on trait dissimilarity. Because trait dissimilarity and environmental dissimilarity were not bounded between 0 and 1, all distance matrices were normalized between 0 and 1 prior to analysis. Because the MRM function does not allow direct model assessment, we assessed model quality and assumptions (e.g., residual normality) by performing a standard multiple linear regression with distance matrices as vectors, which provides equivalent results to the MRM function (SI Appendix, Fig. S14).

To compare observed results to those expected under neutral evolution, we simulated trait dissimilarity using Brownian motion. A Brownian motion model was first fit to the observed phylogeny and the corresponding sigma value was used for the simulations. Following Mazel et al. (25), we simulated species’ positions on the first axis of the trait space (PCoA 1) using the R package geiger (61). We repeated this process 999 times and, in each repetition, we recalculated trait dissimilarity (distance between ecoregion centroids) and reexamined the relationship between trait dissimilarity and taxonomic, phylogenetic, and environmental dissimilarity to develop null distributions of R² values between variables. We then compared observed R² values, also calculated using trait dissimilarity on the first PCoA axis only, to the null distributions to determine whether the observed relationship was stronger than neutral expectations. We also calculated the magnitude of difference between observed and expected trait dissimilarity [i.e., SES, termed “phylogenetically standardized trait distance” in Mazel et al. (25)] as observed dissimilarity—mean expected dissimilarity/SD of expected dissimilarities—and examined the linear relationship between SES and environmental dissimilarity. Ecoregion pairs were considered more similar or less similar than expected if observed dissimilarity was less than or greater than 97.5% of simulated dissimilarity values (two tails used to allow for both “convergent” and “divergent” pairs). We simulated species’ positions on the first axis only because PCoA axes are perfectly uncorrelated by construction and simulating multiple PCoA axes under Brownian motion would induce strong correlations, leading to bias in the neutral data (see ref. 25).

We tested whether latitudinal similarity in PCoA 1 centroids was robust to the choice and number of traits by deleting each trait, one at a time, recalculating PCoA 1 centroids for each ecoregion, and reexamining latitudinal patterns. Latitudinal similarity was again robust to all choices (and in some cases became more pronounced), indicating that the results were not dictated by any single trait but rather the overall composition of traits (SI Appendix, Fig. S15).

We tested whether latitudinal similarity in PCoA 1 centroids as well as the relationship between trait dissimilarity and environmental dissimilarity were robust to removing extreme centroid values—highly positive PCoA 1 centroids at high latitudes (i.e., ecoregions with the highest leverage). We sequentially removed each of the top 15 ecoregions with the highest PCoA 1 centroid values and reexamined latitudinal patterns for the remaining PCoA 1 centroids, as well the relationship between trait dissimilarity and environmental dissimilarity. However, only 13 ecoregions were removed for assessing the relationship between trait dissimilarity and environmental dissimilarity, because Malvinas/Falklands and Channels/Fjords of Southern Chile were not included in the original MRM analysis (see Influences of Taxonomic, Phylogenetic, and Environmental Dissimilarity). Latitudinal similarity in PCoA 1 centroids (i.e., unimodal relationship between centroids and latitude) was again robust to removing ecoregions with the most extreme values; however, the unimodal pattern became less pronounced as ecoregions were removed (SI Appendix, Fig. S16). The strong linear relationship between trait dissimilarity and environmental dissimilarity was highly robust and remained clear even with 13 ecoregions removed. (SI Appendix, Fig. S17).

We tested whether the results of the MRM analysis were robust to the choice of environmental variables used to calculate environmental dissimilarity. We removed each environmental variable, one at a time, recalculated environmental dissimilarity, and reran the MRM analysis. We also tested the effect of adding pH and salinity, both individually and together. Salinity data were obtained from the MULTIOBS_GLO_PHY_REP_015_002 database (62) available through Copernicus Marine Environment Monitoring Service (https://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=MULTIOBS_GLO_PHY_REP_015_002), and pH data were derived from the Norwegian Earth System Model forced ocean simulation (63). Again, the results were robust to all combinations of environmental variables, with mean effect sizes (across all models) of 0.60 (±0.03) for environmental dissimilarity, 0.13 (±0.04) for phylogenetic dissimilarity, and 0.09 (±0.04) for taxonomic dissimilarity (SI Appendix, Table S4).

Because not all taxa included in this study were identified to species level, and because different species in the same family may have similar traits, we tested whether results of the MRM analysis were robust to calculating taxonomic dissimilarity at the family level rather than species level. Thus, we recalculated taxonomic dissimilarity between ecoregions and reran the MRM analysis modeling trait dissimilarity in function of family-level taxonomic dissimilarity, phylogenetic dissimilarity, and environmental dissimilarity. Again, results were robust with effect sizes of 0.54 for environmental dissimilarity, 0.11 for phylogenetic dissimilarity, and 0.09 for taxonomic dissimilarity.

We tested whether latitudinal similarity in PCoA centroids could have resulted from latitudinal patterns in species richness using two null models where we randomly shuffled species composition among ecoregions while maintaining initial species richness and species occurrence frequency. We therefore shuffled species composition, recalculated PCoA centroid values for each ecoregion, and calculated the strength (R² value) of the relationship between randomized PCoA 1 centroid values and latitude using 1) second-order polynomial regression and 2) simple linear regression with the absolute value of latitude. We repeated this process 999 times to build null distributions of R² values and compared observed R² values to the null distributions. R² values from both observed models were greater than 95% of the values in the null distributions, indicating that latitudinal trait similarity was not a result of latitudinal patterns in species richness alone (SI Appendix, Fig. S18).

In this study, we used species occurrences rather than species abundances (see Fish Survey Methods). To determine whether results were robust to this choice, we recalculated PCoA 1 centroids for all ecoregions as abundance-weighted centroids and tested the correlation between the original centroid values and the abundance-weighted centroid values. We then reran the MRM analysis in which trait dissimilarity was calculated as the Euclidean distance between abundance-weighted ecoregion centroids in the first four dimensions of trait space, and taxonomic turnover was calculated with species abundances using Bray–Curtis dissimilarity. The Pearson’s correlation value between the original (occurrence-based) PCoA 1 centroid values and the abundance-weighted PCoA 1 centroid values was 0.97, demonstrating that trait similarity was robust to both species’ occurrences and abundances. Effect sizes from the MRM analysis, and R² values from individual regressions, using species abundances were nearly identical to those of the original occurrence-based analysis (Original analysis: environmental dissimilarity: effect size = 0.60 and R² = 0.56; taxonomic dissimilarity: effect size = 0.09 and R² = 0.20; phylogenetic dissimilarity: effect size = 0.11 and R² = 0.24. Abundance-based analysis: environmental dissimilarity: effect size = 0.51 and R² = 0.56; taxonomic dissimilarity: effect size = 0.003 and R² = 0.12; phylogenetic dissimilarity: effect size = 0.15 and R² = 0.28), again demonstrating that results were robust to both occurrence and abundance data.

We examined whether results were robust to how trait dissimilarity was quantified by rerunning the MRM analysis using functional beta diversity from Chao et al. (64) instead of distances between ecoregion centroids in trait space, where functional beta diversity is calculated as a generalized version of pairwise Rao's quadratic entropy between assemblages. In agreement with the original MRM analysis, and individual regressions, environmental dissimilarity had the highest effect size (0.31) and individual R² value (0.28), while taxonomic and phylogenetic dissimilarity had weaker effect sizes (taxonomic: −0.012; phylogenetic: −0.002) and R² values (taxonomic: 0.04; phylogenetic: 0.06) and added little explanatory power (R² full model: 0.281).