Assessing Benthic Responses to Fishing Disturbance Over Broad Spatial Scales That Incorporate High Environmental Variation

Lundquist, Carolyn J.; Bowden, David; Cartner, Katie; Stephenson, Fabrice; Tuck, Ian; Hewitt, Judi E.

doi:10.3389/fmars.2018.00405

ORIGINAL RESEARCH article

Front. Mar. Sci., 29 November 2018

Sec. Marine Conservation and Sustainability

Volume 5 - 2018 | https://doi.org/10.3389/fmars.2018.00405

Assessing Benthic Responses to Fishing Disturbance Over Broad Spatial Scales That Incorporate High Environmental Variation

$\r\nCarolyn J. Lundquist,*$ Carolyn J. Lundquist^1,2*

David Bowden³

Katie Cartner¹

Fabrice Stephenson¹

Ian Tuck^2,4

Judi E. Hewitt^1,2

¹National Institute of Water and Atmospheric Research Limited, Hamilton, New Zealand
²Joint Graduate School in Coastal and Marine Science, University of Auckland, Auckland, New Zealand
³National Institute of Water and Atmospheric Research Limited, Wellington, New Zealand
⁴National Institute of Water and Atmospheric Research Limited, Auckland, New Zealand

Marine benthic habitats are modified by a number of human-related disturbances. When these disturbances occur at large scales over areas of high environmental variability, it is difficult to assess impacts using metrics such as species richness or individual species distributions because of varying species-specific responses to environmental drivers (e.g., exposure, sediment, temperature). Impact assessment can also be problematic when assessed at broad spatial scales because of regional heterogeneity of species pools. Even when effects on individual species can be detected, it is difficult to upscale from individual species to ecosystem scale effects. Here, we use a functional group approach to assess broad scale patterns in ecological processes with respect to fishing and environmental drivers. We used data from field surveys of benthic communities from two large, widely separated areas in New Zealand’s EEZ (Chatham Rise and Challenger Plateau). We assigned 828 taxonomic units (most identified to species) into functional groups related to important ecosystem processes and likely sensitivity to, and recovery from, fishing disturbance to the seafloor. These included: opportunistic early colonists; substrate stabilisers (e.g., tube mat formers); substrate destabilisers; shell hash-creating species; emergent epifauna; burrowers; and predators and scavengers. Effects of fishing disturbance on benthic functional composition were observed, even at this broad spatial scale. Responses varied between functional groups, with some being tolerant of fishing impacts and others showing rapid declines with minimal fishing effort. The use of a functional group approach facilitates assessment of impacts across regions and species, allowing for improved generalisations of impacts to inform management and decision making.

Introduction

Assessment of impacts of human and natural disturbances on seafloor communities is challenging (Kaiser et al., 2006), particularly in soft sediments where diffuse and chronic disturbances cover broad areas. Benthic macrofaunal communities in marine soft sediment ecosystems are highly diverse, with high spatial variability, and abundant biogenic structures produced by both epifaunal and infaunal species on both continental shelves and coastal soft sediment ecosystems (Snelgrove, 1999; Gray, 2002; Ellingsen et al., 2007). The diversity of species-specific responses to both disturbances and the underlying natural environmental variability means that overarching metrics such as species richness may be inconclusive when averaged across diverse species pools (Bremner, 2008). Sampling over broad spatial scales can further confound interpretation of responses particularly when sampling effort per unit area is limited, such as frequently occurs in regional surveys of marine systems (Hewitt et al., 2010).

Functional group approaches are a potential tool to elucidate the roles of particular species in maintaining ecosystem structure and function across regionally heterogeneous species pools, and inform generalisations about the scales of disturbance at which diverse communities can persist (Bonsdorff and Pearson, 1999; Bremner, 2008). To allow general predictions about seafloor community dynamics that are composed of hundreds of individual species, species can be grouped by biological traits, such as provision of biogenic structure, body size, feeding mode, or mobility. Common to most functional group approaches is identifying a relevant combination of biological traits that reflect key ecological processes and performance of ecosystem function (Bonsdorff and Pearson, 1999; Bremner et al., 2006; Mouillot et al., 2013).

The analysis of biological traits has proved to be a useful approach, highlighting community responses which are difficult to detect using an individual or key species approach, due to high natural spatial and temporal variability in individual species abundance (e.g., Bremner et al., 2003; Tillin et al., 2006; de Juan et al., 2007, 2009). Here, we define biological traits at a species level as biological characteristics that can be measured in all species at a location. To define functional groups based on their roles in the soft sediment ecosystems and their potential effects on ecosystem functioning, we use ‘response’ traits that characterise a species’ ability to survive with respect to different environmental conditions and ‘effects’ traits that reflect each species’ effects on the environment, and indirectly on other species through ecological interactions such as substrate modification (Beauchard et al., 2017). The functional group approach enables comparison of responses to environmental impacts among regions that may have large differences in species composition (Bremner, 2008; Hewitt et al., 2008). Functional group approaches provide advantages from approaches based solely on traits in that they avoid over-inflation of individual traits, and are not limited to solely examining impact responses that may have no relation to ecological function (Beauchard et al., 2017). Changes in the functional composition of communities with disturbance are hypothesised to be consequences of the adaptations of organisms to their environment and their response to stress (de Juan et al., 2007). As such, disturbance is likely to cause shifts in the relative proportion of functional groups, as sets of species with particular traits are replaced by other sets with traits better suited to the disturbance regime (Mouillot et al., 2013). As biodiversity loss is increasingly associated with loss of ecosystem services, functional group approaches also provide a context for predicting changes in community structure and the services associated with communities within a context of increasing disturbance intensity and frequency (Bremner, 2008; Cardinale et al., 2012).

Physical disturbances to the seafloor by fishing trawls and dredges can modify habitats and reduce biodiversity and productivity through their effects on a range of species, including those that provide biogenic structure (e.g., sponges, tubeworms, anemones), or modify the sediment (Hall, 1994; Dayton et al., 1995; Thrush and Dayton, 2002; Hiddink et al., 2006; Kaiser et al., 2006). Changes in community composition are not limited to decreases in abundances; rather, in some cases, significant increases in abundance and/or biomass of scavenging or other opportunistic species can occur, leading to further impacts on community structure through species interactions (Lambert et al., 2017). Biological traits (e.g., morphology, life history, dispersal characteristics) often correlate both with the sensitivity of different species to disturbance, and with their ability to colonise disturbed habitats. Biological traits associated with sensitivity to bottom disturbance by trawling and dredging include being sedentary or having low mobility while protruding from the sediment, and having fragile body forms such as shells or branched morphology (Thrush et al., 1998; Thrush and Dayton, 2002; Hewitt et al., 2011a). In contrast, mobile species, those with robust body morphology, and small opportunistic species have in some experiments shown positive responses to disturbance (e.g., Ramsay et al., 1998; Collie et al., 2000; Lindholm et al., 2001). Sensitivity to, and rate of recovery from, disturbance also depend on life history traits such as generation times and potential for rapid recolonisation of disturbed areas via larval dispersal (Auster and Langton, 1999; Bremner et al., 2003; Tillin et al., 2006; Lundquist et al., 2010; Thrush et al., 2013; Lambert et al., 2014).

Here, we use an extensive dataset of benthic faunal distributions covering > 200,000 km² on Chatham Rise and Challenger Plateau, New Zealand, to investigate the suitability of a functional group approach to determining impacts of bottom fishing disturbance on benthic communities. We investigate patterns of abundance of individuals within functional groups in relation to fishing effort at broad spatial scales, as well as with physical correlates of natural disturbance and other environmental drivers that may be influencing abundance. Eight functional groups were defined through expert elicitation, representing key aspects of the way organisms in seafloor communities modify their environment and interact with each other, and how they would respond to and recover from disturbance to the seafloor (Lundquist et al., 2013). These groups include: (1) opportunistic early colonists with limited substrate disturbance; (2) opportunistic early colonists with considerable substrate disturbance; (3) substrate stabilisers (e.g., tube mat formers); (4) substrate destabilisers; (5) shell hash-creating species; (6) emergent epifauna; (7) burrowers; and (8) predators and scavengers. We hypothesise that opportunistic early colonists and predators and scavengers (groups 1, 2, and 8) would have neutral or positive responses to disturbance, due to increased supply of prey and detrital material (Ramsay et al., 1998; Tillin et al., 2006; Hewitt et al., 2016); that shell hash creators and deep burrowing fauna (groups 5 and 7) would show some declines with fishing effort, but be resilient due to having hard shells or being able to burrow deeply (Hewitt et al., 2016); and that substrate stabilisers and destabilisers (groups 3 and 4), both of which are primarily surface or sub-surface dwellers and likely subject to direct physical impacts of trawls, would show declines in abundance with increasing fishing effort (Lohrer et al., 2004). Finally, we hypothesised that emergent epifauna (group 6) would exhibit the largest negative response to fishing disturbance due to their location on the surface, and sensitivity to physical impacts due to their erect and often fragile morphology (Thrush et al., 2001).

Materials and Methods

Chatham Rise and Challenger Plateau Benthic Community Datasets

Chatham Rise is east of New Zealand in the Pacific Ocean, while Challenger Plateau is in the Tasman Sea, to the west of New Zealand. Both regions were sampled as part of Ocean Survey 20/20 voyages that surveyed extensive areas (>200,000 km²) of New Zealand’s EEZ (Bowden, 2011) (Figure 1). To characterise benthic invertebrate community assemblages, we used samples collected using two different gear types on voyages occurring in 2007 (Chatham Rise: TAN0705 April 2007; Challenger Plateau: TAN0707 June 2007): video camera transects using NIWA’s Deep Towed Imaging System (DTIS) (Hill, 2009) and epibenthic sleds (SEL) (Bowden, 2011; Clark and Stewart, 2016). Data from both gear types were used in the functional group analysis to enable inclusion of both smaller and infaunal taxa that were sampled better by SEL, and larger, more widely distributed species that were sampled better by DTIS.

FIGURE 1

FIGURE 1. Sites sampled for benthic invertebrate assemblages across Chatham Rise (East of the South Island) and Challenger Plateau (West of the North Island) during Ocean Survey 20/20 voyages TAN0705 and TAN0707. Filled dots were sampled with both DTIS video and SEL epibenthic sleds; open circles were sampled with DTIS video only.

Video transect data (n = 147 sites; 107 on Chatham Rise and 40 on Challenger Plateau) included counts of mega- and macro-epifauna (sized > 50 mm), biogenic features, and substrate types extracted from along the entire length of each video transect. DTIS was deployed for 1 h on the seabed per transect (approximate tow dimensions of 1.5 m by 1 km in depths ranging 98 – 1425 m), and abundances were standardised to number of individuals 1500 m^-2 of seabed. Taxonomic resolution from the video transects was generally available to species level (Hewitt et al., 2011b). SEL datasets [n = 123; Chatham Rise (n = 76 samples covering 64 sites) and Challenger Plateau (n = 47 samples covering 41 sites)] enumerated infauna and epifauna (sized > 25 mm), identified to the lowest practicable taxonomic unit, usually to species (Hewitt et al., 2011b). Epibenthic sleds (∼1 m wide) were deployed for 15 min on the seabed per tow (approximately 1 km long in depths ranging 64 – 1950 m). Abundances were left as counts per deployment because the area covered by the epibenthic sled gear is not generally considered as reliably quantitative due to potential for gear saturation before the end of the tow.

Chatham Rise and Challenger Plateau Environmental Datasets

Environmental variables were available for a range of metrics from the New Zealand Marine Environments Classification database (Table 1), at a 1 km² grid scale within the sampling regions (Hadfield et al., 2002; Snelder et al., 2006). Modelled datasets included tidal current velocity, seabed slope, seabed rugosity, and primary productivity. Metrics directly sampled during the OS 20/20 voyages included depth, sediment grain size (from multi-corer samples), and substrate type observations from video (Table 1).

TABLE 1

TABLE 1. Environmental variables available for Chatham Rise and Challenger Plateau.

Chatham Rise and Challenger Plateau Fishing Effort

Extensive areas of Chatham Rise and parts of Challenger Plateau are commercially fished. Data on the intensity and frequency of bottom trawling by commercial fishing vessels were sourced from the New Zealand Ministry for Primary Industries for the 16 years from 1989–1990 to 2004–2005 at depths to 1600 m, the maximum recorded trawl depth in the region during this period (Baird et al., 2011). Fishing effort in the Ocean Survey 20/20 survey region was estimated from individual tow polygons representing the estimated swept area of each trawl event. The cumulative total area swept by trawl events was calculated and standardised as the cumulative proportion of area trawled within a 25 km² resolution grid. It is important to note that no sites were located in the areas of highest fishing intensity (the maximum cumulative swept area being about 70% on Chatham Rise), because the OS 20/20 surveys were not designed to sample across gradients of fishing effort (Baird et al., 2011; Hewitt et al., 2011a).

Fishing effort was generally lower on Challenger Plateau than on Chatham Rise (Baird et al., 2011), but OS 20/20 sampling sites covered a similar range of fishing effort for both regions, with the majority of sites within lower fishing effort regions and few replicate sites in high fishing effort regions (Table 2).

TABLE 2

TABLE 2. Fishing effort classes for Ocean Survey 20/20 Chatham Rise and Challenger Plateau datasets used to compare functional group abundance at site replicates.

Derivation of Trait-Based Conceptual Groups

Functional groups were selected through expert elicitation (as described in Lundquist et al., 2013) as the key components of seafloor communities to explore the response to benthic fishing disturbance, and effects of changing the functional group composition on facilitation or hindering of recovery to disturbance through their interactions with other species and with the seafloor. Traits associated with each functional group were then selected that define key ecological responses to stressors and effects of each group on the seafloor environment and on ecological functioning. A fuzzy logic trait approach was used to allocate taxa into functional groups based on biological trait categories and modalities that were perceived to define responses to disturbance and effects on potential recovery from disturbance (Table 3). Taxa were allocated to functional groups based on multiple biological traits that were selected to include both response traits that characterise the ability of each functional group to survive in response to different environmental conditions or disturbances, and effects traits that reflect how each functional group modifies the environment and creates either opportunities or barriers to recolonisation and recovery after disturbance (Table 3). Allocations of taxa (n = 828) to traits were through an expert-derived database of biological traits based on morphology and available natural history information (Bremner et al., 2006; Ellingsen et al., 2007; Hewitt et al., 2008). Traits included body size, trophic modes, mobility, habitat structure, morphological form, bioturbation ability, and sediment depth at which the animal is typically found (de Juan et al., 2007; Hewitt et al., 2008) (Table 3).

TABLE 3

TABLE 3. General description of fauna classified under eight conceptual functional groups based on key functional traitmodalities.

Fuzzy coding (e.g., Chevenet et al., 1994) was used to describe the relative affinity of a taxon for a given trait modality, ranging between 0 and 1 and summing to 1 across the trait category. For example, the feeding trait category included the trait modalities of deposit feeders, suspension feeders and predator/scavengers. A taxon that is solely a predator/scavenger would receive a score of 1 for this trait modality, whereas one that uses an equal mix of deposit and suspension feeding would receive a score of 0.5 for each trait modality. The criterion for conceptual group membership was that the integrated trait measure for a taxon must exceed a predefined threshold based on the sum of multiple traits describing that conceptual group. Thresholds typically required values of 1 for each trait modality defining a particular functional group (Table 3), except for functional group classifications that included feeding traits where values of 0.5 for feeding traits were determined sufficient to define functional group membership. For example, functional group 4 (substrate destabilisers) required a minimum threshold of 4 composed of values of 1 for traits of ‘freely motile’ for mobility, ‘not’ small sized (i.e., medium or large body size), living within surface sediments, and likely to result in surface mixing of sediment; deposit feeding and ‘not’ short-lived were also common traits of this group. Microalgae, macroalgae, salps, and foraminifera were excluded from the analysis, representing both poor functional group allocation of these groups, and inconsistent sampling as these taxa were often discarded and not enumerated.

Recognising the benefits and limitations of both DTIS and SEL datasets at enumerating all functional groups (see Results), we combined data from both sampling methods to derive functional group composition at each sampling location. Relative abundances of functional groups at each site were calculated by summing abundance across all species allocated to each functional group at a site, calculated independently for DTIS and SEL samples. For sites at which both DTIS and SEL were used (n = 95), we used the maximum value of abundance for a given functional group, recognising that some functional groups would be better estimated by DTIS video, while others (e.g., small infauna) would be better estimated by epibenthic sled samples. To determine whether we could include sites that were only sampled with one methodology (thus increasing replication across different disturbance and environmental strata), we compared the functional composition at these sites using a non-parametric analysis of similarities using PRIMER software version 6, detecting no significant difference between sampling approaches (ANOSIM Global R = 0.045, p = 0.18). Thus, all further analyses were based on data from sites with either SEL, or DTIS, or both. As Chatham Rise and Challenger Plateau datasets showed similarity in terms of relative abundance across functional groups, results represent analyses based on combined data from both regions.

Statistical Analysis

The effect of trawling on the abundances of individual functional groups was investigated. As a linear response was not expected for all functional groups, and sampling locations were distributed unequally across the fishing effort gradient observed in the regions, a categorical variable of fishing intensity was developed, using five classes based on visual analysis of data clusters of trawl effort. Fishing classes represented a gradient from no fishing (‘absent’, 0%) and very low (0.01 – 1% cumulative fishing effort within a cell over 16 years of available trawl data), to high fishing (rates of 25–57% of the area in a cell trawled within a 16 year period; Table 2). Differences in the abundance of functional groups between fishing classes were determined using generalised linear models with a Poisson error structure, a log link function and a quasi-likelihood scale. To investigate correlations between fishing disturbance classes (Table 2) and the variability in abundances within functional groups, we used Pearson’s R on raw and log transformed data and included 2 degree polynomial transformations.

Potential drivers of functional communities including both fishing metrics and environmental drivers were examined using canonical ordination (DISTLM; Primer software version 6 plus PERMANOVA; Anderson and Robinson, 2003) based on Bray-Curtis similarities of square root transformed abundances of the eight functional groups. Continuous variables included as potential explanatory factors were: VGPM; roughness; tidal current; maximum depth; slope; proportion of the area swept by trawling over 16 years; maximum count and effort of trawling over the 16 years; and total count and effort for the proceeding 5, 10, and 16 years. Sediment data (% weight of particles sized > 500 μm, 250–500 μm, 125–250 μm, 63–125 μm, and <63 μm) were only available for a limited subset (63% of sites) so two ordinations were run: one on the full dataset without sediment information as explanatory variables and one on the data subset for which detailed sediment information was available. Variables were normalised and log transformations were included for variables where bivariate scatter plots suggested semi-log relationships between abundance of specific functional groups and the variable. Backward selection using AIC as a stopping criterion was used to select the variables most important in explaining variability in the functional trait composition. Categorical variables were developed for six environmental variables to allow visualisation of patterns of abundance of functional groups with each environmental variable: depth (5 classes); seabed roughness (5 classes); sediment (4 classes); seabed slope (5 classes); tidal current (6 classes); VGPM as indicator of primary productivity (4 classes).

Results

General Patterns in Community Structure

In general, SEL datasets sampled a larger total number of taxa than DTIS, with more diversity of taxa within functional groups. Relative proportional abundances of each functional group were typically similar for those sites that were sampled using both DTIS and SEL, though taxonomic differentiation was typically higher in SEL with taxonomic entities typically identified to species, whereas DTIS methods allowed for identification to functional group, but not always to species (i.e., orange globulous sponge #1). In functional group 6 (emergent epifauna), for example, more than 100 individual taxa were recorded from DTIS samples, but more than 200 were recorded from SEL samples. Functional groups differed in relative abundance with 156 and 298 individual taxa (typically species) enumerated for group 4 and 6, respectively. In particular, small-bodied taxa (primarily groups 1 and 2) were poorly enumerated in both SEL and DTIS data, resulting in few taxa being classified to these two smaller, opportunistic functional groups (Table 3). As neither sampling method available at these broad scales did an adequate job of representing the smaller infauna groups, no further analysis of functional groups 1 and 2 was performed. Few taxa were enumerated in functional group 3 (tube mat formers) but these taxa were considered to be sampled effectively by both DTIS (tube mats typically clearly visible) and SEL. Therefore, limited taxonomic diversity in group 3 was assumed to reflect the restrictive functional trait characteristics defining this functional group, and group 3 was retained within the analysis.

Impacts of Fishing Effort on Functional Group Abundance

Response to fishing disturbance differed between the eight functional groups, generally supporting our hypotheses (Figure 2). The generalised linear models reported overall significant differences between fishing effort classes for functional groups 4, 5, 6, and 8 (p < 0.001 for all groups), but not for group 3 (p = 0.152) or group 7 (p = 0.134). For functional group 4 there was no difference between the three classes of least disturbance (Absent, Very low, Low), but the abundances in these classes were significantly higher than the abundances in the Medium effort class, which in turn was significantly higher than the High effort class (Figure 2). For functional groups 5 and 6, abundances in the lowest fishing effort class (Absent) were significantly higher than all other fishing effort classes (Figure 2). For functional group 8, abundances in the Very low – Medium fishing effort class were significantly higher than the lowest and the highest fishing effort class (Absent and High respectively) (Figure 2).

FIGURE 2

FIGURE 2. Abundance of each functional group for the Ocean Survey 20/20 offshore dataset for different fishing effort classes. Abundance values for groups 4 and 6 are plotted on the secondary y axis. Error bars represent one standard error.

Variability of abundances within each functional group was also affected by fishing disturbance (Table 4). The coefficient of variation for functional group abundances within fishing effort classes was lowest for the High effort class for functional groups 3 – 7 and lowest for the Absent fishing effort class for functional group 8.

TABLE 4

TABLE 4. Coefficient of variation (CV) between fishing effort classes and functional group abundance.

Impacts of Fishing Effort and Environmental Drivers on Functional Group Abundance

Not surprisingly, abundances of the six functional groups showed complex and non-linear patterns with respect to individual environmental variables (Figures 3A–F). For example, functional group 4 showed a discontinuous depth distribution with near absence within the primary trawl depths for the hoki fishery (400–800 m) (Figure 3A), whereas this group exhibited monotonic declines in abundance with increasing trawl intensity (Figure 2). In contrast, functional group 5 (predominantly predatory and scavenging whelks and volutes) exhibited higher average abundances at these heavily trawled depths on the continental slope (Figure 3A).

FIGURE 3

FIGURE 3. Abundance of each functional group for the Ocean Survey 20/20 offshore dataset for environmental variables: (A) depth; (B) seabed roughness; (C) sediment; (D) seabed slope; (E) tidal current; (F) VGPM as indicator of primary productivity. Abundance values for groups 4 and 6 are plotted on the secondary y axis. Error bars represent one standard error.

Only minor amounts of overall variability in functional composition were explained by the environmental variables and the continuous fishing effort data used in the ordinations, with 3 and 7.5% of the variability explained for the full dataset and for the sediment data subset, respectively. However, in both cases, fishing effort data was significantly related to variation in functional community composition. For the full dataset, the only variables selected by the DISTLM to explain the variability in functional composition were proportion of the area swept by trawling over 16 years and the maximum count of fishing trawls in a cell over the same time period. For the sediment data subset, six variables were selected: maximum depth; proportion of the area swept by trawling over 16 years; total area swept by trawling over the same time period; sand; mud; and tidal current. Forcing the model to only use the fishing effort related variables (proportion of area swept and total area swept) resulted in 3% explained, the same amount as for the full dataset. A forward selection procedure, restricted to the six variables, selected sand first (1.2% explained), followed by mud (3.3%), total area swept (0.9%), proportion swept (1.4%), depth (0.07%) and finally tidal current (0.02%). These results suggest overlapping effects between fishing effort and the other variables and between tidal current, depth and the other variables. However, the overlap between fishing effort and the other variables was small as the difference between the amount explained when only fishing effort variables were used (3%) was little different to that when the other variables were allowed to be selected first [0.9% (total area swept) + 1.4% (proportion swept) = 2.3% cf. 3%].

Discussion

Our results suggest that we can use functional groups to generalise predictions of the impacts of seafloor trawling disturbance across benthic communities within diverse seafloor communities across broad geographic regions. Even with potentially confounding factors of high environmental variability (depth, sediment, exposure) and differing regional species pools, the functional group approach resulted in fishing effort being selected as an important predictor of benthic functional group composition in New Zealand’s EEZ. The low percent explanation is not uncommon in multivariate ordination techniques, due to (i) the inclusion of a large gradient across both fishing and environmental variability, and (ii) the limited number of samples, including both intrusive (benthic sled) and non-intrusive (video) sampling methods (ter Braak and Verdonschot, 1995; Hewitt et al., 2002; Forcino et al., 2015). However, this does not preclude determination of key drivers of functional group abundance within current limitations of large scale studies in natural marine systems. Importantly the responses of the functional groups were generally non-linear, both to the environmental variables and fishing effort, and reflected hypothesised differences between different functional groups with respect to their sensitivity to seafloor disturbance.

Our analyses suggest shifts in functional composition with increased fishing effort, with largest declines in functional groups representing emergent epifauna and substrate destabilisers (e.g., surface burrowers like Echinocardium spp.), and shifts in other groups as shown in Figure 2. Emergent epifauna (functional group 6), and to a lesser extent, shell-hash creators (functional group 5), showed the least tolerance, with large and immediate declines to lower abundance even at low rates of fishing effort, implying high sensitivity to fishing effort as evidenced throughout the literature for this functional group (Turner et al., 1999; Jennings et al., 2001; Thrush et al., 2001). Sensitivity to disturbance was also exhibited by substrate destabilisers (functional group 4), with high variability in functional group abundance for the lowest three classes of fishing effort, and significant declines in abundance occurring only at the highest two fishing effort classes.

These different functional responses and reductions in functional groups at different levels of fishing effort suggest increased functional homogenisation with disturbance, as less tolerant functional groups are excluded with increasing rates of disturbance, and functional diversity is reduced (Bremner et al., 2003; Thrush et al., 2006; Hewitt et al., 2010). Decreased variability in functional group abundance was observed with increasing fishing disturbance for all but functional group 8, a group that we predicted to either be unaffected or to respond positively to fishing disturbance (Lambert et al., 2017), again suggesting a homogenising effect of trawling on functional group representation.

Declines in abundance and composition of particular functional groups support broad generalisations about the effects of fishing on ecosystem functions. Emergent epifauna showed the largest declines with increasing fishing effort, and declines in this group would be predicted to result in declines in biogenic habitat structure, the provision of settlement habitat and refuge sites from predation, modification of biogeochemical processes and exchanges, modification of flows, and sediment stabilisation (Auster and Langton, 1999; Turner et al., 1999; Thrush et al., 2001; Diaz et al., 2003; Hewitt et al., 2005; Scharf et al., 2006). Declines in abundance of burrowing fauna would be predicted to result in reductions in ecosystem services performed by this group which include increasing the oxygen flow into the sediment and extending the total zone of denitrification, stimulating nutrient cycling (Beaumont et al., 2007; Norling et al., 2007). Long-lived, large-bodied groups support ecosystem productivity and food production services (Jennings et al., 2001; Beaumont et al., 2007).

Others have shown significant environmental drivers that are associated with either sensitivity or resilience to seafloor disturbance (Bremner et al., 2006). Common examples include coarse sediments and high tidal currents, which are assumed to result in high natural disturbance rates, and thus species that are tolerant of disturbance regimes (Pitcher et al., 2017). Our results instead suggest overlapping effects between fishing effort and environmental variables such as tidal current, depth and sediment. This effect is partially due to the broad spatial scale covered by our analysis, i.e., combining samples across a wide range of depths, exposures, and sediments to examine general patterns in functional group declines with fishing effort. However, it is also likely to be in part due to the somewhat homogeneous nature of the region analysed, with over half of our samples being characterised as low current flow environments with mud or muddy-sand substrates.

High resolution spatial data on fishing activity have been available since 1989 for the study areas but commercial bottom-contact trawl fishing goes back to the early 1970s (McKenzie, 2017). Thus, Chatham Rise, in particular, is a chronically disturbed system and effects stemming from fishing effort not accounted for in our analyses may weaken the strength of detected effects. For instance, large, long-lived sessile benthic fauna will have been largely removed by direct impacts in the early days of trawl use, and widespread trawl-mediated sedimentation beyond the seabed contact footprint of the fisheries will have effects on distributions of suspension-feeding fauna (Clark et al., 2016). The cumulative fishing footprint since 1989 includes 24% of the trawlable depths (i.e., <1600 m) in the New Zealand Exclusive Economic Zone, with approximately 3% of the fishable depths trawled each year (Baird and Wood, 2018). However, the opportunistic sampling that we used for this analysis was not arranged across a gradient of fishing intensity, and the environmental gradients sampled unfortunately did not include the highest fishing intensities. Regardless, sampling of deeper marine soft sediment habitats are expensive and logistically challenging, and are rarely matched to both fishing and environmental gradients at scales of replication sufficient for determining significant impacts. It is encouraging that our methods were able to elucidate fishing as a significant driver of benthic community composition even across the large spatial scale and environmental gradients that were sampled, and using a combination of both intrusive and video techniques, suggesting that remote sensing can be used to successfully interpret responses of benthic communities to environmental and human-induced seafloor disturbances.

The functional group approach is useful and provides a first, high-level, analysis showing effects of bottom contact trawl fisheries. Some challenges in our broad analysis include the minimal replication across the full range of fishing disturbance in the benthic datasets, and low replication of samples within coarser sediment regions. Despite these factors, our results still indicate significant impacts of fishing effort, and the lack of coarse sediments being selected as an important predictor suggests this will not be a major factor affecting our ability to detect fishing impacts. Others have suggested (e.g., Lambert et al., 2014; Pitcher et al., 2017) that areas of coarse sediments or high tidal currents, often in shallow areas subject to frequent natural disturbance, appear to recover more quickly because they generally are associated with lower abundance of macrofauna, and are composed of species with functional traits that result in lower sensitivities to physical disturbances. In contrast, another New Zealand study found that when the coarse sediment is shell hash, as is the case in most of the sites sampled here, generally higher abundances of macrofauna, and long lived surface dwellers in particular, are found (Thrush et al., 1998).

The sampling methodology did not accurately represent all functional groups, particularly smaller opportunistic taxa. We combined video and infaunal (epibenthic sled) sampling to allow us to obtain broader scale information on more widely distributed, larger taxa such as epifauna, while concurrently collecting adequate data to quantify abundance of smaller, infaunal taxa, but recognise that our infaunal sampling (mesh = 25 mm) was not sufficient to estimate abundance of small, opportunistic species. Regardless, our analysis, conducted with both intrusive (benthic sled) and non-intrusive (video) sampling techniques provides indications of validity of using remote techniques to determine benthic community structure and potential responses to impacts across both large spatial scales, and in deep areas where remote sampling is the most cost-effective option.

The conceptual functional groups used here represent broadly applicable functional roles important for assessing impacts on and recovery from seafloor disturbances, and also ecosystem processes that affect integrity. As such they would be useful to apply to a variety of situations and the results should be easily transferable to different regions where functional trait data are available. The responses to disturbance we observed within functional groups across broad spatial scales and regional species pools suggests that we can populate a conceptual model of a generalised recovery trajectory in order to advance our understanding of the role of disturbance on marine community dynamics, as is indeed our intention (Lundquist et al., 2010, 2013). Results from this functional group approach can then be extrapolated to predict the impacts of varying disturbance rates on ecosystem function across an entire seascape based on changes in diversity and abundance of different functional groups at different rates of fishing disturbance.

Author Contributions

CL and JH developed the conceptual design and performed the statistical analyses. KC and FS assisted with data collation and preparation of tables and figures. DB, JH, and IT designed the initial OS 20/20 sampling strategy. CL drafted the manuscript, assisted by DB, FS, IT, and JH.

Funding

Funding to support this analysis was provided by the National Institute of Water and Atmospheric Research, Coasts & Oceans Research Programme, Project COME1601. This analysis utilises data compiled to validate a seafloor disturbance model funded by the Biodiversity Research Advisory Group of New Zealand’s Ministry for Primary Industries (Project #ZBD200925). Data collection and sample identification underlying the analyses within the Ocean Survey 20/20 Chatham-Challenger Post-Voyage Analysis Project were supported by a collaboration between the Ministry for Primary Industries, Land Information New Zealand, and the Department of Conservation.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Mary Livingston and the Ministry for Primary Industries, Biodiversity Research Advisory Group, for supporting this research, and for constructive advice throughout this project.

References

Anderson, M. J., and Robinson, J. (2003). Generalized discriminant analysis based on distances. Aust. N. Z. J. Stat. 45, 301–318. doi: 10.1111/1467-842X.00285