Phylogenetic and functional evidence suggests that deep-ocean ecosystems are highly sensitive to environmental change and direct human disturbance

An understanding of the balance of interspecific competition and the physical environment in structuring organismal communities is crucial because those communities structured primarily by their physical environment typically exhibit greater sensitivity to environmental change than those structured predominantly by competitive interactions. Here, using detailed phylogenetic and functional information, we investigate this question in macrofaunal assemblages from Northwest Atlantic Ocean continental slopes, a high seas region projected to experience substantial environmental change through the current century. We demonstrate assemblages to be both phylogenetically and functionally under-dispersed, and thus conclude that the physical environment, not competition, may dominate in structuring deep-ocean communities. Further, we find temperature and bottom trawling intensity to be among the environmental factors significantly related to assemblage diversity. These results hint that deep-ocean communities are highly sensitive to their physical environment and vulnerable to environmental perturbation, including by direct disturbance through fishing, and indirectly through the changes brought about by climate change.


Introduction
Competition between species has long been recognized as an important factor determining the ecological diversity and structure of organismal communities [1][2][3][4]. Intense interspecific competition for scarce resources can result in the exclusion of certain taxa [2,4], shaping species' realized niches and distributions, and influencing ecosystem functioning [5,6]. Numerous studies have investigated the role of competition in structuring terrestrial (such as plant [7]) and shallow-water (such as coral reef [8]) assemblages, but the importance of interspecific competition in structuring the expansive and functionally important communities of the deep ocean has been a matter of debate because the discovery of high alpha diversity in deep-water sediments [9][10][11]. Some researchers have emphasized an important role of biological interactions in structuring deep-seafloor communities, theorizing a dynamic balance between competitive forces and predation [12,13]. Others have argued that extensive niche differentiation, coupled with typically low organismal densities and the availability of space, mean that competitive interactions are unlikely to be significant in structuring modern-day deep-ocean communities [9,10,14,15].
Empirical evidence in support of either of these viewpoints is limited. Studies that have examined the morphological or trophic characteristics of deep-sea assemblages present some evidence for the displacement of ecologically similar taxa by their competitive dominants [16][17][18]. Conversely, studies that have investigated the taxonomic or phylogenetic structure of assemblages provide some evidence that variation in physical environmental conditions may be of greater influence than competitive interactions in structuring deep-ocean communities [19,20]. However, investigations to date have been limited in their analytical scope by the challenges associated with sampling and/or conducting experiments in the deep ocean.
Knowing whether interspecific competition or the physical environment dominates in structuring natural communities is important because it further enhances our understanding of the sensitivity of ecosystems to environmental change. Communities that are structured predominantly by interspecific competition are typically more stable under environmental stress than those whose structure is governed by the physical environment [21]. Understanding the sensitivity of deep-sea ecosystems to environmental change is of pressing concern because they are predicted to experience increasing direct and indirect anthropogenic pressure over this century [22][23][24]. For example, fishing fleets are operating at ever-increasing depths [25], there is growing commercial interest in the mining of seabed minerals [26], and greenhouse gas emissions are increasing oceanic temperature, reducing pH and dissolved oxygen concentrations, and altering food supply to the deep ocean [23,24,27].
In this study, we perform 'community phylogenetic' analyses to investigate the importance of interspecific competition versus the physical environment in shaping the composition of deep-seafloor assemblages in the Northwest Atlantic Ocean-a region predicted to experience particularly rapid environmental change over this century [28]. Application of a community phylogenetic approach avoids many of the problems associated with conducting experiments in deep-ocean environments, enabling an investigation of previously unprecedented scale. Under this approach, the dispersion of taxa within samples across a phylogeny or functional trait dendrogram is compared with that which would be expected by chance [29]. If taxa within samples are found, on average, to be less similar to one another than would be expected by random draw from the available taxa pool, assemblages are described as 'over-dispersed'; this is typically considered evidence of a dominance of competitive exclusion in shaping assemblage structure, because phylogenetically/functionally similar taxa are assumed to be ecologically similar [3,29,30]. Conversely, if taxa within samples are found, on average, to be more similar to one another than would be expected by chance, assemblages are described as 'under-dispersed'; this is typically considered evidence of a dominance of the physical environment in determining assemblage structure, because phylogenetically/ functionally similar taxa are assumed to share the particular traits that are necessary for survival under the prevailing environmental conditions [3,29,30] (although see [31]).
Our results provide evidence that deep-seafloor communities may typically be both phylogenetically and functionally under-dispersed. We therefore also investigate and discuss a number of physical environmental parameters which may be of importance in structuring the enigmatic but widespread communities of deep-ocean sediments.

Material and methods (a) Sampling of deep-seafloor assemblages
We analysed 312 sediment samples, forming the largest macrofaunal sample set yet collected from the deep ocean. Samples were collected with a box corer (area 0.25 m 2 ) from the continental slopes of the Northwest Atlantic Ocean (depth range: 582-2294 m) (figure 1) between May-August 2009 and June-August 2010, and form part of the international 'NEREIDA' programme (https://www.nafo.int/About-us/International-Cooperation), a project instigated by the Northwest Atlantic Fisheries Organisation (NAFO) in order to investigate the impacts of high seas fisheries on Vulnerable Marine Ecosystems. Sediment subsamples were taken for geochemical and particle size analyses and remaining sediment was washed over a 1 mm mesh sieve for faunal analyses; 20 245 specimens of peracarid crustacean were identified to the genus level (177 genera within 74 families, in total). Peracarid crustaceans were chosen for analysis because of their low dispersal potential, extremely high taxonomic diversity, superabundance in marine sediments and ecological importance as prey, predators and ecosystem engineers [32][33][34].  To investigate the phylogenetic structure of the sampled peracarid assemblages, we constructed a 'supertree' [35] (figure 2a). We used Google Scholar to identify 59 studies containing suitable evolutionary source trees. From each study, only unique source trees were retained for analysis to ensure that there was no duplication of source tree topology (an emergent characteristic of evolutionary trees as phylogenetic hypotheses, and the information directly used during supertree construction) that would otherwise unfairly weight the analysis as a result of pseudo-replication [35]. One hundred and twenty-seven evolutionary trees were retained for analysis (electronic supplementary material, table S2), and monophyletic taxonomic groups were labelled following World Register of Marine Species (WoRMS) systematic nomenclature.
Supertrees were constructed using MULTILEVELSUPERTREE (MLS) v. 1.0 [36] run on the Oxford University Advanced Research Computing supercomputer ARCUS (Phase B) (http:// www.arc.ox.ac.uk/content/home). Because of prohibitive run times, the program was run individually for the peracarid orders Amphipoda, Isopoda, Tanaidacea and Cumacea. For each run, a taxonomy tree was used to guide the program.
To provide reference branch length information for the supertree, we constructed two further phylogenetic trees based on 18S SSU rDNA, 16S rDNA, cytochrome c oxidase 1 (COI) and histone H3 gene sequences downloaded from GenBank (electronic supplementary material, appendix S1; the topologies of these trees are available upon request). Genes were aligned individually using MAFFT v. 7.273 [40] running on the MAFFT online server (http://mafft.cbrc.jp/alignment/server/). Alignments were scrutinized using TRIMAl v. 1.2 [41]. The final alignment was concatenated using SEQUENCEMATRIX v. 1.8 [42] and consisted of 285 taxa and 2586 bp. PARTITIONFINDER v. 1.1.1 [43] was used to select the most appropriate model of evolution and partitioning scheme. ML and Bayesian topologies were estimated using RAXML v. 8.2.8 and MRBAYES v. 3.2.6, respectively, on the CIPRES Science Gateway v. 3.3 online server [44].
During a fifth round of supertree construction, we used MLS 1.0 to combine the output trees from the four previous MLS runs and the ML and Bayesian analyses with all source trees focusing on order-level relationships within Peracarida and Malacostraca to produce a final supertree with 1487 terminal taxa. This supertree topology was then trimmed to include only those taxa present in the GenBank sequence concatenated alignment. We estimated maximum-likelihood branch lengths for this topology using RAXML v. 8.2.8. Common nodes between the supertree and ML branch length tree were labelled using PHYLOCOM v. 4.2 [45] and the labelled ML branch length tree was used as an input for the program R8S v. 1.8 [46] in order to obtain a dated phylogeny. Twenty-two nodes were constrained with age estimates based on fossil data (electronic supplementary material, table S1). Non-parametric rate smoothing (NPRS) with Powell optimization was selected as the analysis method. The BLADJ function of PHYLOCOM [45] was then used to obtain a fully dated supertree (figure 2a; see electronic supplementary material for a 'Newick' format representation of the supertree to enable detailed examination of its topology).

(c) Functional dendrogram construction
To characterize the functional structure of the sampled peracarid assemblages, we constructed a dendrogram (figure 2b) describing the functional similarity of sampled families based on their scoring for a selection of traits (electronic supplementary material, table S3). Trait groupings and traits were chosen based on ecological relevance and data availability. We used a fuzzy coding [47]    rspb.royalsocietypublishing.org Proc. R. Soc. B 285: 20180923 [48] to select the most appropriate distance metric and clustering method as Euclidean distance and unweighted pair group method using arithmetic averages (UPGMA) clustering. Analyses were conducted in R v. 3.0.2 [49].

(d) Testing for phylogenetic and functional assemblage structure
We investigated phylogenetic assemblage structure using the phylostruct function of the R package 'picante 1.6 -2' [50] based on the constructed supertree (figure 2a) and complete genuslevel peracarid assemblage matrix. To investigate assemblage structure at smaller spatial scales, we used ESRI ARCGIS v. 10.1 to produce seven data subsets, each consisting of 30 box cores chosen at random from within a set radius (50, 100, 200, 300, 400, 500 and 600 km) of the central-most sampling point. To quantify the phylogenetic dispersion (diversity) of peracarid genera within each sample, we calculated the metric 'phylogenetic species variability' (PSV) [51]. We employed permutation tests (100 000 permutations) to determine whether the average PSV of all samples, and sample subsets, was significantly different from that expected under two null hypotheses-'Null 1' and 'Null 2' [51]. Under Null 1, phylogenetic structure was removed from both taxon prevalence and sample composition by the randomization of taxon presence within samples. Under Null 2, phylogenetic structure was removed from sample composition, but not from taxon prevalence, by the randomization of taxon occurrence between samples.
Since functional dendrograms are analogous to phylogenies in form, we employed the same methods as outlined above to quantify the functional dispersion (referred to herein as 'functional species variability' (FSV)) of peracarids contained within each sample based on the constructed functional dendrogram (figure 2b) and compared this to the expected outcome under the two null hypotheses stated above.
(e) Testing for phylogenetic signal As the interpretation of community phylogenetic patterns relies on knowledge of the evolution of taxon traits [29], we tested for phylogenetic signal across the peracarid trait matrix by applying a Mantel test [52]. Based on the constructed supertree (figure 2a), we calculated the square root of patristic distance as a measure of phylogenetic distance between taxa using the R package 'ape 4.1' [53]. Euclidean distance was calculated as a measure of trait similarity between taxa based on the peracarid trait matrix. The Mantel test was performed using the R package 'vegan 2.0-9' [54] (100 000 permutations). To assess the strength of phylogenetic signal in individual traits, based on the supertree of Peracarida (figure 2a), we calculated Pagel's l [55] for each trait in the peracarid functional trait table using the R package 'phylosignal 1.2' [56] (100 000 permutations).

(f ) Characterization of the deep-sea physical environment
To investigate relationships between the PSV/FSV of sampled peracarid assemblages and the prevailing environmental conditions, we examined the following environmental parameters: bathymetry (depth, slope, aspect, seafloor rugosity, bathymetric position index (BPI)); fishing intensity (vessel monitoring system (VMS) signal density and total trawl length per km 2 ); geological context; seafloor sediment particle size ( percentage clay/silt/sand); carbon availability (percentage inorganic, organic and total carbon, surface chlorophyll a and particulate organic carbon (POC) concentrations, modelled transport of POC to depth); physical oceanographic variables (temperature, salinity and current speed); and month and year of sample collection.
Water depth at each sampling location was extracted using ARCGIS 10.1 based on multibeam bathymetric surveys (5625 m 2 cell size). Slope, eastness and northness, roughness (225 Â 225 m analysis window) and standard deviation of multibeam bathymetry values (225 Â 225 m analysis window) were calculated using the Spatial Analyst extension of ARCGIS 10.1. Benthic Terrain Modeler [57] was used to calculate BPI over a range of radii as well as seafloor rugosity (375 Â 375 m and 1875 Â 1875 m analysis windows).
We quantified bottom trawling intensity using VMS signal locations. Individual trawl paths were identified based on boat identity, speed, location, date and time using ARCGIS 10.1, and the Line Density Tool (Spatial Analyst extension) was used to measure the total length of trawls per km 2 within a set radius (1, 3 or 5 km) from each box core.
The sediments of the study area were classified into 12 discreet geological categories based on their acoustic characteristics, depth and slope [58]. We extracted the relevant geological category for each sampling location using ARCGIS 10.1.
Sediment per cent clay/silt/sand was calculated for each core subsample based the following particle size categories consistent with the 'Phi' (F) scale. We calculated particle size diversity following the methodology of Etter & Grassle [59] and Leduc et al. [60].
Sediment total carbon and organic carbon content were determined using a Leco TruSpec CHN analyser. Inorganic carbon was determined by the difference between the total carbon and organic carbon measurements for each sample [61].
We obtained surface chlorophyll a and POC concentrations from the Giovanni ocean colour radiometry online data system (https://giovanni.gsfc.nasa.gov/giovanni/). MODIS AQUA 4 km resolution data were downloaded for the years 2008 -2010. These data were interpolated to 2500 Â 2500 pixels (525 m resolution) in QGIS 2.2. We estimated POC delivery to the seafloor from surface POC concentrations following region-specific equations [62].
Seafloor temperature, salinity, and meridional and zonal current speed values were extracted from a modelled monthly average data layer for the study area [63] and averaged both across the year prior to sample collection and across the year of sample collection. These values were interpolated to 3000 Â 3000 pixels (578 m resolution) in QGIS 2.2. Absolute current speed was calculated using Pythagoras's theorem; 10-year minimum/maximum values for temperature and current speed, respectively, and 10-year average values for both variables were calculated to capture longer term variability.

(g) Statistical analyses
We removed highly correlated environmental variables following consideration of variance inflation factors (VIFs). Variables were removed from the analysis in a stepwise manner (those with highest VIF first). The resulting dataset contained 19 variables (electronic supplementary material, table S4). VIF calculations were undertaken in R 3.0.2 [49] using the package 'HH 3.1-32' [64].
We constructed generalized additive models (GAMs) using the R package 'mgcv 1. criterion [65]. Explanatory terms included in each GAM were refined by backwards stepwise selection considering variable p-values and model AIC until a minimum AIC value was reached (see the electronic supplementary materials for additional detail relating to the methodology employed by this study).

(a) Phylogenetic and functional structure of deepseafloor assemblages
We found the average PSV value of the deep-sea assemblages analysed to be significantly smaller than that which would be expected under both null hypotheses (mean PSV observed ¼ 0.8728; mean PSV Null 1 ¼ 0.8817, mean PSV Null 2 ¼ 0.8830; probability mean PSV observed taken from Null 1 distribution , 0.001; figure 3a; probability mean PSV observed taken from the Null 2 distribution , 0.0001; figure 3b). Similar results were obtained for all data subsets analysed. Based on these results, we conclude that the deep-sea assemblages analysed are phylogenetically 'under-dispersed' (i.e. that the peracarid taxa within a sample are on average more closely related to one another than would be expected by chance).
Further, we found the average FSV value of the assemblages analysed also to be significantly smaller than that which would be expected under both null hypotheses (mean FSV observed ¼ 0.8220; mean FSV Null 1 ¼ 0.8354, mean FSV Null 2 ¼ 0.8259; probability mean FSV observed taken from Null 1 distribution , 0.0001; figure 3c; probability mean FSV observed taken from the Null 2 distribution , 0.0001; figure 3d). Similar results were obtained for all data subsets analysed. We therefore conclude that the deep-sea assemblages analysed are functionally 'under-dispersed' (i.e. that the peracarid taxa within a sample share, on average, more functional traits with each other than would be expected by chance).

(b) Phylogenetic signal
Significant phylogenetic signal was identified in the peracarid trait matrix (Mantel test: r ¼ 0.3583, p , 0.001; Pagel's l: 34 of 38 traits exhibit significant phylogenetic signal; electronic supplementary material, table S5). We thus conclude that, on average, phylogenetically similar taxa tend also to be functionally similar. We found assemblage PSV to relate negatively with average seafloor temperature for the year of sample collection ( p ¼ 0.0349; figure 4a), and maximal current speed values for 10 years prior to sample collection ( p ¼ 0.0214; figure 4b). However, we found PSV to vary unimodally with average surface chlorophyll a concentration for the year of sample collection ( p ¼ 0.0011; figure 4c), and in a complex but weakly positive manner with sediment organic carbon content ( p ¼ 0.0001; figure 4d). We also found assemblage PSV to be significantly related to the month ( p , 0.0001; figure 4e) and year ( p , 0.0001; figure 4f ) of sample collection. We found assemblage FSV to be negatively related to bottom trawling intensity ( p ¼ 0.0452; figure 4g), but positively   figure 4i). Sediment total carbon content was found to relate in a complex but weakly positive manner with FSV ( p ¼ 0.0185; figure 4j). We also found assemblage FSV to be significantly related to the month of sample collection ( p , 0.0001; figure 4k).

Discussion
The complementary phylogenetic and functional analyses performed here provide evidence for a compositional underdispersion of the focal deep-sea assemblages at the spatial scales investigated. This under-dispersion may reflect the selection of favourable phenotypic traits that are shared between similar taxa [29]. Typically, the selecting agent in question is the physical environment, and, as a result, this process is known as 'environmental filtering' [29]. Although not global in extent, and focusing only on continental slope depths, the results of our study provide evidence that the physical environment may be more important than interspecific competition in shaping the composition of deep-water communities, emphasizing a potentially high sensitivity of deep-sea ecosystems to environmental perturbation [21].
Our findings challenge those studies that hypothesize an importance of competition and character displacement as a significant ecological structuring agent in the deep sea [12,13], and conflict with investigations that have examined the morphological or trophic characteristics of deep-sea assemblages [16][17][18]. Instead, they substantiate the hypothesis that competitive interactions between species in the deep ocean are weak and unlikely to be significant in structuring communities at the spatial scales investigated [14,15]. Further, they support the results of previous analyses of lesser spatial scope that have investigated the taxonomic and phylogenetic structure of deep-sea assemblages [19,20].
Comparison of our results with those of other studies employing the PSV metric suggests that the phylogenetic signal observed here is comparably strong or stronger than that reported for many non-marine assemblages, including temperate lake fish assemblages [30,51], tropical plant assemblages [66], archaea assemblages [67] and tropical bird assemblages [68].
Although significant phylogenetic signal was apparent in the functional trait matrix, our results provide some evidence for the convergence of functional traits between relatively distantly related crustacean taxa. For example, in figure 3, mean FSV Null 2 is closer to mean FSV observed than mean FSV Null 1 is ( figure 3c,d ), suggesting that a portion of the observed underdispersion reflects elevated functional similarity of the most prevalent taxa (Null 2 removes phylogenetic/functional structure only from assemblage composition, maintaining any structure in relative taxon prevalence, while Null 1 removes phylogenetic/functional structure from both assemblage composition and taxon prevalence). However, while the more prevalent taxa are more functionally similar to each other than would be expected by chance, they are not correspondingly phylogenetically similar; mean PSV Null 2 is not closer to mean PSV observed than mean PSV Null 1 is ( figure 3a,b), suggesting that their functional similarity is convergent to an extent. Examples of apparent functional convergence can be identified in figure 2b-the amphipod family Lysianassidae, and the isopod family Cirolanidae, for example. Our results suggest that traits related to fecundity and armament exhibit greatest propensity for convergent evolution among the peracarid taxa analysed (electronic supplementary material, table S5).
Our investigation into the possible environmental drivers of assemblage variability demonstrates that both natural and anthropogenic factors may influence the structure of the deep-sea assemblages (figure 4). The negative relationship between bottom trawling intensity and assemblage FSV (figure 4g) suggests that physical disturbance by bottom trawling reduces soft-sediment functional diversity, with the resulting assemblages exhibiting a reduced subset of the functional traits that would otherwise be present in an undisturbed assemblage. While generally concordant with the small number of studies that have investigated trawling impacts on deep-sea macrofauna and meiofauna [69][70][71][72][73], this finding adds a new facet to our understanding of the impacts of bottom trawling in the deep ocean.
Our analyses demonstrate a negative relationship between seafloor temperature and the phylogenetic diversity of the sampled peracarid assemblages (figure 4a), indicating that the physiological tolerances of peracarid taxa to temperature change are preserved within evolutionary lineages. That this relationship is apparent across a temperature range of only approximately 1.28C (figure 4a) suggests that even the superficially small increases in deep-ocean temperature that are predicted to occur over this century as a result of climate change [74,75], particularly in the high seas of the North Atlantic [28], will significantly reduce the phylogenetic diversity of the communities found there, potentially impacting deep-ocean ecosystem functioning [76]. An altered phylogenetic profile of deep-sea ecosystems may eventually lead to a change in the cycling, storage and sequestration pathways of nutrients and chemicals, such as carbon.
Under current climate change scenarios, global patterns of the export of surface production to the deep ocean are expected to change in a complex manner [23,24,27,77]. Food supply to the deep ocean may dwindle in some regions, such as the North and South Atlantic Oceans [24], while being enhanced in others, such as the Arctic and Southern Oceans [24]. Our analyses suggest that changes in food availability in the deep ocean may affect both the phylogenetic and functional variability of communities, but in a complex manner (figure 4c,d,j ). This, in turn, may affect the availability and variety of food for demersal and pelagic organisms that feed on sediment-dwelling prey. This multifaceted relationship is complex and still poorly understood.
Overall, the results of our analyses suggest that deep-water soft-sediment ecosystems, which constitute the majority of global seafloor area, may be particularly sensitive to environmental change. Such ecosystems are central to a number of important ecosystem services, including carbon sequestration [78], and are predicted to come under increasing direct and indirect anthropogenic pressures [22][23][24]27]. Even superficially small changes in natural and anthropogenic disturbance regimes, temperature, food availability and bathymetry may significantly alter the phylogenetic and functional variability of deep-seafloor communities (figure 4), and this may alter ecosystem functioning and the provision of ecosystem services by the deep ocean [76]. Our findings are therefore relevant to the understanding of anthropogenic pressures on deep-sea ecosystems, including, for example, the prediction of possible mining impacts on deep-sea fauna. We advocate that the precautionary principle be exercised in all circumstances where anthropogenic actions may disrupt the natural ecology of deep-water ecosystems. rspb.royalsocietypublishing.org Proc. R. Soc. B 285: 20180923 Ethics. This research was performed in accordance with all applicable international, national and/or institutional laws, guidelines and ethical standards. Competing interests. The authors declare no competing interests. Funding. This research was supported a Natural Environment Research Council (NERC) Collaborative Awards in Science and Engineering (CASE) studentship (NE/K006886/1), by Merton College, University of Oxford, UK and through the European Union's Horizon 2020 research and innovation programme under grant agreement no. 678760 (ATLAS). This output reflects only the authors' views and the European Union cannot be held responsible for any use that may be made of the information contained therein. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. A.D.R. also acknowledges the Oxford Martin School in the funding of the Oxford Martin School Sustainable Oceans Project, which has also contributed to the current paper.