Landscape-scale spatial heterogeneity in phytodetrital cover and megafauna biomass in the abyss links to modest topographic variation

Sinking particulate organic matter (POM, phytodetritus) is the principal limiting resource for deep-sea life. However, little is known about spatial variation in POM supply to the abyssal seafloor, which is frequently assumed to be homogenous. In reality, the abyss has a highly complex landscape with millions of hills and mountains. Here, we show a significant increase in seabed POM % cover (by ~1.05 times), and a large significant increase in megafauna biomass (by ~2.5 times), on abyssal hill terrain in comparison to the surrounding plain. These differences are substantially greater than predicted by current models linking water depth to POM supply or benthic biomass. Our observed variations in POM % cover (phytodetritus), megafauna biomass, sediment total organic carbon and total nitrogen, sedimentology, and benthic boundary layer turbidity, all appear to be consistent with topographically enhanced current speeds driving these enhancements. The effects are detectable with bathymetric elevations of only 10 s of metres above the surrounding plain. These results imply considerable unquantified heterogeneity in global ecology.

The ocean basin-scale distribution and abundance of life on the abyssal seafloor is set extensively by the supply of sinking particulate organic matter (POM) from overlying surface waters, a limiting food resource for most deep-sea life [1][2][3][4] . At regional to local scales, it is suggested that bathymetry, habitat terrain type, and lateral transport of POM play important roles in determining the distribution of biomass and biological assemblage composition [5][6][7][8][9] . About 85% of the world's seafloor lies at abyssal depths, and contains many millions of hill features 10,11 . If these bathymetric features, rising only 10 s to 100 s of metres above the abyssal plain, support differing quantities of biomass, this intermediate heterogeneity would be important to global biogeochemistry, ecology, and biogeography.
As a result of remineralisation in the water column, particulate organic carbon (POC) flux decreases with increasing water depth 12 . The Martin curve 13 allows for the estimation of variation in POC flux with water depth (Z), i.e. flux ∝ Z −0.7 (Z in metres) 14 . Similarly, the corresponding decline in the expected standing stock of seafloor biomass supported by that flux can be linked to water depth via other studies, e.g. megafauna biomass ∝ 10 −0.4Z (Z in kilometres; Supplementary Text) 3 . The vertical fluxes of POC, particle volume, and particle mass over the Porcupine Abyssal Plain (PAP) are all highly correlated 3,15 . These relationships provide a null framework against which to judge the potential impact of abyssal hill terrain on the supply of POM to the seafloor and consequently biomass. In this contribution we assess POM supply as the areal cover (%) of the seafloor by POM aggregates, frequently referred to as phytodetritus, which can be visually distinguished from the seafloor sediment surface 16,17 . Biomass Distribution. Megafauna biomass exhibited highly significant variation between depth bands, with the two deepest bands significantly different from all shallower bands, and from each other ( Table 1,  Supplementary Table 1). These latter contrasts were of substantial magnitude (d = 1.4-1.9, "very large" practical significance). There was a factor of 2.5 increase in megafauna biomass from the plain to the shallowest depth band (Fig. 2B). Similarly, a separate assessment of the small hill and surrounding plain in Area D (Fig. 1A) indicated significant variation with depth, with the top of the small hill having significantly higher biomass than deeper bands (Supplementary Table 1). Comparison of the high resolution survey areas (Fig. 1A, Areas A-C) indicated highly significant variation between hill and plain areas, but no significant difference between the two plain areas (Supplementary Table 1). Primary terrain classes exhibited highly significant differences in biomass, with all three classes significantly different from one another, Hill the highest and Plain the lowest (Table 1, Fig. 3A,C). Similarly, biomass varied among secondary terrain classes, with the 'Plain' class significantly different from all others with the exception of the 'Slope D' class (Table 1, Fig. 3B,D).
Turbidity and Seabed Sediments. Benthic boundary layer (BBL, water column < 10 m above seafloor) turbidity (suspended particulate load) exhibited a modest but marked increase over the elevated terrain, with values some 1.04-1.09 times higher than the plain (Fig. 2C). Water column turbidity > 10 m above the seafloor was broadly consistent over the depth range examined and very similar to BBL turbidity over the abyssal plain (Fig. 2C). Across the 21 sites sampled with a Megacorer (Fig. 1A), sediment mud content (Mud %), total nitrogen content (TN %), and total organic carbon content (TOC %) exhibited significant positive correlations with water depth (Table 2; Fig. 2D-F). Similarly, TN % and TOC % exhibited significant positive correlations with Mud % (p < 0.05, p < 0.01, respectively). Comparisons of all three parameters between abyssal plain (sites > 4840 m) and elevated terrain (sites < 4840 m) locations indicated significant differences (p < 0.02) in all cases (Table 2). In contrast, the TOC-to-TN ratio (C/N) exhibited no significant correlation with water depth, and no significant difference between abyssal plain and elevated terrain locations (Table 2). Environment-POM-Biomass. Links between POM cover, biomass, and water depth were assessed by reference to predictions based on the influence of water depth alone ( Fig. 2A,B). In each case, the observed increase in POM cover and biomass at shallower depths was greater than the null framework predicted, with POM cover increasing by a factor of 1.05 and biomass by a factor of 2.5 in comparison to the predicted values of 1.01 and 1.07 respectively. Inter-relationships between POM cover, biomass, and environmental variables were assessed by Spearman's rank correlation (Table 3). POM cover and biomass exhibited significant correlations with each other and with most other variables tested (Table 3). We also examined which relationships remained significant after controlling the influence of other variables through partial correlation. The partial correlation between POM cover and biomass was notably non-significant, and biomass only exhibited a significant but moderate negative partial correlation with water depth (Table 3). POM cover continued to display an appreciable negative partial correlation with both water depth and distance from hill crest, and a lesser positive relationship with seabed slope angle.

Discussion
We detected substantial spatial variation in both POM cover and megafauna biomass in relation to abyssal hill terrain on the PAP. The results suggested that water depth, terrain type, and distance from hill crest were important in determining the distribution of POM cover at the seabed. The interaction of bottom currents with seabed terrain is most likely the mechanism through which POM cover was controlled. Various measures of local terrain elevation (water depth, BPI, seabed slope angle, etc.) provide reasonable proxy variables for the complex interactions between bottom currents and seafloor features 9,19,23,24 .
The influence of bottom current speed on POM cover may be non-linear. For example, above a threshold speed, POM is resuspended and potentially removed, while below threshold, POM accumulation at the seabed may be positively correlated with current speed. Seafloor time-lapse photography and near-bottom current meter data from the Porcupine Seabight 23 (4,025 m water depth) showed resuspension of POM at c. 7 cm s −1 . During Hill; Area C, Southern Plain; Area D, small hill and its surrounding plain). Seabed sampling sites are also indicated (*). (B) Spatial variation in seabed POM cover is illustrated by colour-coded symbols for those map cells having 50 or more data values. The primary hill-bounding bathymetric contour (4840 m), and the defined primary terrain types are also illustrated. Map projection is UTM WGS 1984 zone 28N (ArcGIS v10.2, Environmental Systems Research Institute); bathymetric contours and primary terrain classes have been edited to remove small features (noise) over the abyssal plain (note, un-edited data were used in all analyses, and are illustrated in Supplementary Fig. 1).
the period of our study, current speed was monitored at 30-minute intervals 100 m above the abyssal plain c. 6 km to the east of our seabed survey area 21 ( Supplementary Fig. 4), the mean speed recorded was 3.8 cm s −1 , with 8% of observations exceeding 7 cm s −1 , time-integrated flow was to the SE and highest mean current speeds directed to the SSE. Clearly, speeds may have exceeded this threshold more frequently over the elevated terrain features studied. Similarly, current speeds are likely to have been generally reduced in the BBL over the abyssal plain by comparison to the water column 100 m above the plain. Modified bottom water current speed, direction, and turbulence, have been observed over elevated areas elsewhere, including small abyssal hills 5,25,26 . Our observation of significantly enhanced BBL turbidity over elevated terrain is also suggestive of increased speed/turbulence, and consequently resuspension of POM. Systematic tidal variation in current speed and direction ( Supplementary  Fig. 5) could introduce a temporal bias to our observations; however, this seems unlikely given the extended period (9-days) of our observations. For example, in our comparisons of the high-resolution study areas (Fig. 1A, Areas A-C), each represents 6-hours of continuous observation, such that the full range of tidal current speed is likely to have been encompassed in each area.
Downslope, gravity-related transport of POM at the seafloor is important on the continental slope, and in major geomorphological features such as canyons and trenches 6,7,27 . It is conceivable that some of the differences and trends in POM cover we observed were influenced by downslope transport processes. Nevertheless, POM cover was consistently higher on elevated terrain than on the plain in all of our comparisons, suggesting that downslope transport was of limited significance in our case.
Our finding that Hill terrain had 2.5 times the megafauna biomass of Plain terrain was consistent with previous work at the PAP, where a series of three hill summits and one hillside location were found to have a mean megafauna biomass of 3.1 times that of the mean of four plain locations spread up to ~40 km apart 18 . Fishes notably exhibited no observed distribution structure over the hills examined in the present study, as determined from the same AUV missions 20 . However, over 90% of the fish recorded were presumed scavengers including carrion feeders (e.g. Macrouridae), and consequently unlikely to respond to POM cover 20 , i.e. the carrion on which they feed has a lower potential for near-seabed current influence/redistribution. Our megafauna biomass data exhibited substantial negative correlations with water depth and distance from hill crest, though partial correlations suggested that seabed elevation was the strongest single correlate. The observed magnitude of the effect (factor of 2.5, Hill relative to Plain) greatly exceeded that expected from depth-only attenuation of POM at the study site (factor 1.01) 14 , the observed seabed POM cover (factor 1.05), and that expected for megafauna biomass (factor 1.07; Supplementary Text) 2 .
Interpreting these differences in apparent effect scale, 1.01-1.05 in food supply versus 2.5 in standing stock biomass, requires a number of assumptions. By reference to the Metabolic Theory of Ecology 28 , its potential application to seabed biological communities generally 29 , and deep-sea fauna specifically 30 , we can suggest five important points: (a) the megafauna (as quantified here) are likely to represent a substantial fraction of total seafloor biomass; (b) megafauna biomass is likely to scale proportionately with total seafloor biomass; (c) that the PAP megafauna taxa are likely to have lifespans in excess of a decade (typical individual body mass 5 g, able 1. Mean particulate organic matter (POM) cover and megafauna biomass in water depth and terrain classes. Bathymetric position index (BPI) and seabed slope angle define terrain classes as indicated (LT, less than; GT greater than; full classification is provided with Supplementary Fig. 1). (n, number of POM images or number of biomass tiles; fwwt, fresh wet weight; 95% CI, confidence interval of mean value). environmental temperature 2.5 °C), and consequently; that (d) megafauna biomass is likely to scale directly with integrated multi-year food supply (e.g. POM, phytodetritus); and that (e) biomass links to resource consumption (e.g. respiration) 3 . The apparent scale difference between POM cover (1.05), as a proxy for potential food supply, and megafauna biomass (2.5) enhancement factors from Plain to Hill, therefore suggests that a substantial additional organic matter supply supports the observed Hill biomass-that additional supply is most likely laterally transported POM from the adjacent water column, driven by bathymetrically enhanced bottom current speeds.
We assume that instantaneously observed seabed POM cover represents a dynamic balance between supply (vertical and lateral) rate and removal (decomposition and consumption) rate. For example, the apparent consumption of POM (phytodetritus) by some megafauna taxa can lead to the removal of most visible POM within weeks 31,32 . Similarly, resuspended POM may be further consumed within the water column of the BBL, a processes that has previously been observed to occur for 2-3 months following the peak flux period over the PAP 33 . The mechanisms that underpin such bathymetrically-linked relationships are challenging to resolve, in part, because visible accumulations of settling POM on the seafloor are ephemeral often only lasting several weeks, and even then potentially subject to tidal variation 31,[33][34][35] . Given that the vertical flux rate is only expected to be a factor of 1.01 higher on the Hill relative to the Plain, and that removal rate may scale directly with megafauna biomass, then a near-corresponding factor increase is implied in laterally transported organic matter supply (i.e. factor 2.4). The potential significance of laterally transported material is supported by our observation of a c. five-fold increase in megafauna suspension feeder biomass from Plain to Hill. Our observation of increased suspended particulate load (turbidity) in the BBL over the hill is consistent with this difference in suspension feeder biomass. Similar, potential links between bathymetrically enhanced current speeds and seafloor biomass have been recorded in other deep-sea environments 5,18,[36][37][38] .
The likely significance of bathymetrically enhanced bottom current speeds is also suggested in our observations of significant positive correlations between mud, total nitrogen, and total organic carbon content of seabed sediments and water depth-i.e. reduced on Hill relative to Plain, in opposition to the POM and biomass trends. Our results are consistent with previous studies of other hill and plain locations at the PAP 18,19 . The reduced sediment mud content on hills reflects winnowing of the sediments by enhanced current speeds (i.e. loss of fine, silt and clay, particles). This effect is very readily observed on PAP abyssal hills by the ubiquitous presence of ice-rafted dropstones at the sediment surface, that are absent from plain sites 18 ; post-glacial sedimentation has buried the abyssal plain dropstones, but not those on hills where the sediment column has been winnowed by topographically enhanced currents.
The change in sediment mud content (elevated terrain median 74%, plain median 85%; 0-1 cm sediment horizon) that we have observed is likely accompanied by a change in the dominant sediment mineralogy, as has been previously documented for other hill and plain sites at the PAP 39 . This shift in sediment type and mineralogy may be significant when considering the sediment inventory of organic carbon. TOC is frequently associated  Supplementary Fig. 1). with clay particles or minerals 40 , and organic carbon-clay systems are thought to be important in the preservation of organic matter in marine sediments 41 . Our data do exhibit a statistically significant positive relationship between sediment mud content and TOC content (0-1 cm sediment horizon; Spearman's rank correlation r s = 0.583, n = 21, p = 0.006). Turnewitsch et al. 19 have previously considered organic matter supply to an adjacent much larger hill on the PAP, noting that reduced sedimentary organic carbon content may be the result of both (a) reduced deposition of phytodetritus, and (b) reduced organic matter preservation, both processed being driven by topographically enhanced bottom water current speeds. Our data are certainly consistent with the latter process (b). However, over the modest terrain elevations we have studied in detail, we have not detected any reduction in the deposition of phytodetritus (as assessed by instantaneous seabed POM cover). Despite limiting our study to very modest terrain elevations (< 80 m) in deep water (4850 m), we have detected substantive change in seafloor food supply and megafauna biomass, that appears to be driven by bathymetrically enhanced bottom current speeds 42 . The abyssal hill terrain studied supported a biomass 2.5 times greater than the surrounding abyssal plain, suggesting a corresponding enhancement of biological rate processes and ecosystem functions. Given the prevalence of abyssal hill terrain on our planet, and that water depth-only null predictions yield enhancement factors of only 1.0-1.1, we believe that our observations are particularly significant. We anticipate that these results will be mirrored in other hill areas. We expect that the relative magnitude of the increase in both food supply and biomass may increase, to some extent, with greater seabed elevations, acting through locally enhanced bottom water current speeds.

Confidence
This new perspective suggests that elevated terrain may experience magnified responses to change in vertical POC flux from surface waters, as may result from climate change 43 . There may also be practical implications for establishing appropriate environmental baselines, and in making impact assessments, over varied bathymetry in abyssal regions subject to seafloor mining 44 . Our results would also suggest the need for caution in global meta-analyses of the relationships between POM cover, or sedimentary TOC, and benthic standing stocks where local seafloor terrain is not considered. We also note that this new appreciation of hills driving seafloor processes does not alter previous interpretations of substantial temporal change in the seafloor community on the plain at PAP 45 .

Conclusions
Our results demonstrate that changes in bathymetry as small as a few tens of metres can substantially impact the local supply of organic matter to the seafloor and consequently megafauna biomass. There appears to be strong evidence to suggest that this effect is driven by bathymetrically enhanced bottom water currents. The enhanced megafauna biomass on hills is likely mirrored by increased ecological function and carbon cycle roles (e.g. carbon standing stock and total seafloor respiration). Abyssal hill terrain may represent the most widespread landform on the planet 46 , with an estimated 25 million hills of 100-1000 m elevation present globally, and even higher numbers of smaller features suspected 10 . The influence of abyssal hill terrain on deep-sea ecology and carbon cycling is likely to be ubiquitous. Understanding this landscape-scale heterogeneity appears to be vital to the biogeography, ecology and biogeochemistry of the world's seafloor.

Methods
During RRS Discovery cruise 377, 5-27 Jul 2012, a series of Autonomous Underwater Vehicle (AUV) seabed photographic surveys were carried out covering two spatial scales 47 . One comprised a grid having 1 km line spacing over and around an abyssal hill, with a single line extending a further 10 km across the abyssal plain to the south of the hill. The second scale comprised of three grids with 100 m line spacing, one on the abyssal plain immediately to the north of the abyssal hill, one on the hill, and one on the plain 10 km to the south of the hill (Fig. 1A).
Colour images of the seabed were obtained using a vertically-mounted Point Gray Research Inc. Grasshopper 2 camera, with a 2/3″ sensor (2448 × 2048 pixels), attached to the AUV Autosub6000. The camera was fitted with a 12 mm lens, yielding in water acceptance angles of 26.7° and 22.6°. Images were taken every 0.9 seconds from a target altitude of 3.2 m. Full details of the field methodology are given elsewhere 22 .  To estimate megafauna biomass 48 , 64,690 images were mosaicked in groups of 10 to produce 6,469 image tiles each representing ~14 m 2 of seafloor, these tiles were then randomly assigned to multiple investigators 49 , who annotated the tiles to generate taxon-specific numerical density and body size data (methods detailed elsewhere) 22 . Fresh wet weight megafauna biomass was estimated via morphotype-specific equations relating image-measured body dimension to individual body mass 18,50 . Each tile was assigned a terrain category and corresponding environmental variables based on geolocation. Biomass data were log transformed prior to parametric statistical analyses to account for right skew in the data 51 . Total megafauna biomass was partitioned among nominal feeding types based on previous studies of stable isotopes 52 and individual behaviour 53 carried out in the near vicinity.
Seabed cover by particulate organic matter (POM) was estimated from individual (un-mosaicked) images. The images were cropped to partially correct for non-uniform illumination, the remaining image (2248 × 1548 pixels) representing ~1.6 m 2 of seafloor. Percentage seabed cover by POM was estimated using a custom MATLAB (The MathWorks Inc.) routine based on colour image segmentation and resultant object detection measurement (Supplementary Fig. 2; Supplementary Text). Two classes of POM were identified ('light' and 'dark'), this simplified and improved image segmentation, and these two classes were summed to produce 'total' POM cover (%). POM data were logit transformed prior to parametric statistical analyses to account for the proportional nature of the data 54 .
Bathymetric data from the survey area were collected during RRS James Cook cruise 071, 29 Apr-12 May 2012, using a hull-mounted Kongsberg EM120 12 kHz multibeam system 55 , these data were processed and gridded using Caris HIPS and SIPS software (Teledyne CARIS, Inc.) to give a final bathymetric model resolution of 100 m. The bathymetric model data were further processed to yield seabed slope angle, profile curvature, and rugosity via native functions in ArcGIS v10.2 (Environmental Systems Research Institute). Bathymetric position index (BPI), a second order derivative of the model surface, relating local elevation to the wider landscape, was calculated using the package 'Benthic terrain modeller' 56 . BPI was calculated in an annulus neighbourhood at a scale factor of 1200 m (i.e. 12 cell outer radius) with an inner radius of four cells (400 m), and combined with seabed slope angle to distinguish three primary terrain classes: 'Hill' , 'Slope' , and 'Plain' ( Table 1, Fig. 1B). Finer divisions of BPI and slope angle were employed to produce a secondary classification comprising twelve terrain types ( Supplementary Fig. 1).
Variations in megafauna biomass and POM seabed cover were assessed with water depth and terrain variables. Images and tiles were binned into natural 12.5 m depth intervals (i.e. integer multiples of 12.5); to avoid low sample size in the first and last bins, data were amalgamated with adjacent bins (Table 1). Each image and tile was assigned a primary and secondary terrain class based on geolocation. After transformation (see above), both POM cover and megafauna biomass exhibited significant inequalities in variance between environmental categories (p < 0.001; Levene's test 57 ). To acknowledge this heteroscedasticity, and the unbalanced samples sizes, we employed Welch's modification of the one-way ANOVA 58 , and used the Games-Howell method for subsequent pair-wise comparisons 59 . We also considered effect size using Cohen's d-statistic 60 . The potential impact of spatial autocorrelation on our analyses and interpretations was addressed by examination of variograms that suggested the phenomenon occurred but was unlikely to impact our interpretations (Supplementary Text; Supplementary  Fig. 3). We compared our field estimates of biomass and POM cover with water depth only-driven null predictions (Supplementary Text).
For assessment with terrain variables and classes, seabed POM cover and megafauna biomass were averaged in spatial grid cells corresponding with those of our bathymetric model. Only cells with ≥ 50 images (POM) and ≥ 5 tiles (biomass) were included in these analyses. Potential correlations between POM, biomass, and individual terrain parameters were assessed via Spearman's rank correlation 61 . Partial Spearman's rank correlation coefficients were also examined (implemented in R 62 using package 'ppcor' , function 'pcor' 63 ) to provide additional insight into the practical significance of the observed correlations.
The Autosub6000 vehicle was fitted with a Seapoint Turbidity Meter (Seapoint Sensors, Inc.), detecting light (880 nm) scattered (15-150° scattering angle) by suspended particles in close proximity to the sensor (< 5 cm). Sensor data, calibrated to Formazin Turbidity Units (FTU), were recorded at 2-second intervals throughout the course of our surveys. For assessment these data were binned into natural 12.5 m depth intervals (as for analysis of POM cover and megafauna biomass) and assigned to two groups: (a) water column, where vehicle altitude above seafloor was > 10 m, and (b) benthic boundary layer (BBL), where vehicle altitude was ≤ 10 m.
Seabed sediments were sampled at 21 locations ( Fig. 1A; Supplementary Table 2) using a Bowers & Connelly Megacorer 64 during RRS Discovery cruise 377 47 . Sediment particle size distributions were determined by laser diffraction (Malvern Mastersizer) 18 and results reported as mud content (% wt/wt, particles < 63 μ m) for the 0-10 mm sediment depth horizon. Sediment total nitrogen (TN) and total organic carbon (TOC) content (% wt/wt) were analysed in duplicate for the 0-10 mm sediment depth horizon using the vapour-phase acid de-carbonation method 65 . These mud, TN and TOC data were assessed for interrelationships and relationships with water depth via Spearman's rank correlation 61 , and potential relationships illustrated as simple least squares linear fits (Fig. 2D-F). Comparisons between abyssal plain (sites > 4840 m) and elevated terrain (sites < 4840 m) locations were carried out using Mood's median test 61 , with confidence intervals calculated using the Wilcoxon signed rank method 66 , both as implemented in Minitab 17 (Minitab, Inc.).