Effects of human footprint and biophysical factors on the body-size structure of fished marine species.

Marine fisheries in coastal ecosystems in many areas of the world have historically removed large-bodied individuals, potentially impairing ecosystem functioning and the long-term sustainability of fish populations. Reporting on size-based indicators that link to food-web structure can contribute to ecosystem-based management, but the application of these indicators over large (cross-ecosystem) geographical scales has been limited to either fisheries-dependent catch data or diver-based methods restricted to shallow waters (<20 m) that can misrepresent the abundance of large-bodied fished species. We obtained data on the body-size structure of 82 recreationally or commercially targeted marine demersal teleosts from 2904 deployments of baited remote underwater stereo-video (stereo-BRUV). Sampling was at up to 50 m depth and covered approximately 10,000 km of the continental shelf of Australia. Seascape relief, water depth, and human gravity (i.e., a proxy of human impacts) were the strongest predictors of the probability of occurrence of large fishes and the abundance of fishes above the minimum legal size of capture. No-take marine reserves had a positive effect on the abundance of fishes above legal size, although the effect varied across species groups. In contrast, sublegal fishes were best predicted by gradients in sea surface temperature (mean and variance). In areas of low human impact, large fishes were about three times more likely to be encountered and fishes of legal size were approximately five times more abundant. For conspicuous species groups with contrasting habitat, environmental, and biogeographic affinities, abundance of legal-size fishes typically declined as human impact increased. Our large-scale quantitative analyses highlight the combined importance of seascape complexity, regions with low human footprint, and no-take marine reserves in protecting large-bodied fishes across a broad range of species and ecosystem configurations.

INTRODUCTION

Globally, the overexploitation of some coastal ecosystems has resulted in the widespread decline of larger‐bodied fish (Jackson et al. 2001; Olden et al. 2007) that likely play a key role in ecosystem function (Dulvy et al. 2004; Fisher et al. 2010) and population replenishment (Barneche et al. 2018). Regional examples of appropriate fisheries and conservation governance, including ecologically representative networks of no‐take marine reserves (NTMRs) and management tools, such as catch limits (MacNeil et al. 2020), can serve as a blueprint for achieving stock sustainability and marine biodiversity conservation targets (Duarte et al. 2020). In this context, quantitative syntheses that encompass multiple biogeographic regions, capturing broad‐scale variations in environmental conditions, ecosystem structure, and species composition are critical to disentangle anthropogenic, habitat, and environmental determinants of the abundance and size distribution of fishes (Cinner et al. 2018; Jouffray et al. 2019). Ultimately, this information can be used to optimize ecosystem functions and services for human well‐being under different social‐ecological contexts (Cinner et al. 2020).

In commercially exploited fishes, numerous size‐based indicators have been developed in recent decades to track the response of populations and communities to changes in management and environmental conditions (Bianchi et al. 2000; Shin et al. 2005). More recently, these indicators are being measured with nondestructive methods that can be broadly applicable to multispecies assemblages and ecosystem types, across both temperate (Fulton et al. 2005) and tropical regions (e.g., Wilson et al. 2010; Nash et al. 2016; Robinson et al. 2017). However, studies that test the robustness of these indicators over large geographical scales, capturing transitions in ecosystem structure from tropical to temperate waters, have been limited and relied on the use of underwater visual census techniques (UVCs) (e.g., Stuart‐Smith et al. 2017; Heather et al. 2020). Although global in extent, UVCs are limited in depth (usually <20 m) and can potentially underrepresent the diversity and abundance of large‐bodied fished species (Willis & Babcock 2000; Rojo et al. 2021). Further, some targeted fish display behavioral changes to the diver's presence along gradients of human footprints (Lindfield et al. 2014), including those associated with depth (Andradi‐Brown et al. 2018).

The advent of baited remote underwater stereo‐video (stereo‐BRUV) has proven valuable to overcome the limitations of diver‐based censuses in sampling large roving species across broad depth ranges (Langlois et al. 2020). Agonistic behavioral interactions around the bait may bias abundance and size‐frequency distribution estimates (Dunlop et al. 2015); however, the use of large horizontal fields of view provides a representative sample of the broader fish assemblage, including predatory and prey species (Cundy et al. 2017; Coghlan et al. 2017), and provides length‐frequency distributions comparable to fisheries‐dependent surveys (Langlois et al. 2015). Prior to the current study, the large‐scale representative syntheses of fish assemblage data from stereo‐BRUVs were limited to studies focusing on sharks and pelagic predatory species (Letessier et al. 2019; MacNeil et al. 2020). However, demersal marine teleosts encompass a diverse taxonomic group with high social, economic, and ecological importance, representing about one‐third of the global reconstructed marine fisheries catch (FAO 2016) and playing key roles in marine trophodynamics (Frank et al. 2015).

We investigated the potential drivers of body‐size structure of 82 species of marine demersal teleosts that are of importance to commercial, recreational, or both fisheries around Australia, a continent with a record of strong fisheries management (Alder et al. 2010) and one of the most complex and extensive networks of NTMRs worldwide (Boonzaier & Pauly 2016). We compiled length and abundance information from a unique, large‐scale (approximately 10,000 km of coastline) synthesis of stereo‐BRUV data encompassing a section of the continental shelf of Australia (1–50 m depth) (Figure 1). The diverse biogeographic distributions, contrasting life‐history strategies, and fisheries productivity of these species (Appendix S1) enabled us to test the role of human footprints in shaping body‐size structure of fished species. Information‐theoretic approaches were used to assess the relative influence of habitat, environmental, and anthropogenic factors for predicting the probability of occurrence of fishes above the 90th percentile of the observed body size distribution for a given species (i.e., large), as well as the relative abundance of fishes above (i.e., legal) and below (i.e., sublegal) the minimum legal size (MLS) of capture.

FIGURE 1

Locations of baited remote underwater stereo‐video (stereo‐BRUV) deployments (i.e., samples, n = 2904) across Australian neritic (1–50 m) continental shelf systems (green, samples in no‐take marine reserves [NTMRs] boundaries, n = 541; red, samples outside NTMR boundaries, n = 2363; light blue, temperate waters; orange, tropical waters)

METHODS

Data syntheses

We collected quantitative data on the body‐size structure of recreational and commercially fished species (marine teleosts) from a synthesis of stereo‐BRUV data across Australia (Figure 1; Appendix S2). The compilation of individual data sets was enabled through a central repository of image‐annotation data (GlobalArchive; https://globalarchive.org/). The full data set consisted of single and stereo‐BRUV deployments (2004–2017) and covering six administrative states and five marine ecoregions. Survey data were collected and analyzed with standard operating procedures (Langlois et al. 2020) and were mostly comparable in terms of bait type (pilchards), method of recording abundance (maximum number of individuals present in the field of view for each species at one time [MaxN] [Priede et al. 1994]), taxonomic resolution (species level), and soak time (60 min; only approximately 6% of samples had 30 min) (details in Harvey et al. [2021]). To improve the comparability of the data set, we excluded samples that did not include length information for all target species present (i.e., mono BRUVs), targeted offshore oceanic regions, targeted depths beyond the 50‐m isobath, and were not conducted during daylight hours (0900–1700). Small‐scale variations in the stereo‐BRUV setup (sample separation, bait quantity, and soak time) were accounted for by including CampaignID (i.e., a unique sampling event in time and space carried out by a particular research group) as a random factor in the statistical models (Appendix S3). Although MaxN has been criticized for its hyperstability (Schobernd et al. 2014) (i.e., counts remain stable despite changes in true population abundances), this typically occurs at relatively large MaxN values (>20), which are rare for targeted species that typically occur in low numbers (mean maximum MaxN 7.85 [SD 8.14] across species in our data set). Recent empirical analyses have also shown that MaxN scales linearly with metrics that quantify abundance across the entire length of the deployment (e.g., mean count), providing an unbiased estimate of abundance (MacNeil et al. 2020). Where temporal sampling had occurred, we retained the year closest to 2011 (Appendix S4), corresponding to the year the population density grid was available. The temporal coverage of samples we used spanned from 2006 to 2017, and samples were spread across seasons in temperate and tropical waters, although the distribution was unbalanced (Appendix S5). The inclusion of CampaignID as a random factor in the analyses enabled us to account for this spatial and temporal variation (Appendix S4).

Body‐size structure of fished species

We only considered fished species managed through state‐specific MLS regulations to better understand anthropogenic and environmental drivers of body‐size structure in fished species. These included 82 demersal marine teleosts in 18 families (Appendix S1). For each stereo‐BRUV deployment, we aggregated (i.e., summed) individuals of a species at the time of the species’ MaxN that were above the 90th percentile of the observed body distribution (i.e., large), as well as those that were above (i.e., legal) and below (i.e., sublegal) the state and region‐specific MLS. The aggregated abundance (summed MaxN across species) for each size class was then used as a size‐based indicator. Large fishes were rare in the samples, and we converted abundance data to presence or absence and used this as an indicator of fishing impacts in the probability of occurrence of the largest, and more vulnerable, size class. The large fishes were still included in the legal‐size fishes (i.e., not mutually exclusive) to retain abundance records of the largest individuals. Although this metric cannot adequately capture changes across the full spectrum of individual body sizes in a community (e.g., slope of the size spectra [Heather et al. 2020]), it provides a robust approach to assess abundance trends for size classes with different vulnerabilities to fisheries exploitation. Similar indicators have been proposed elsewhere to assess the ecological status of shallow reef ecosystems (Stuart‐Smith et al. 2017) and fisheries potential of global coral reefs (Cinner et al. 2020). Because some species are not managed through MLS regulations in all states, we used the median MLS across states as a conservative threshold for defining legal and sublegal fishes. This was a necessary trade‐off to capture interregional differences in life‐history parameters, as well as social considerations, that determine the legal‐size of capture, while preserving the majority of fished species contributing to the data across regions.

We tested whether assemblage‐level models were representative of community‐level responses to human footprints by partitioning our analysis for regional groups of conspicuous species. For each representative tropical and temperate family in the data set, we selected either the most conspicuous targeted species or grouped species with similar habitat preferences and suspected levels of fisheries exploitation at the genus level. The latter was a necessary trade‐off due to the large amount of 0s for some species across their biogeographic extent of occurrence. This resulted in seven regional species groups, including key targeted species for Australian commercial and recreational fisheries (Flood et al. 2014), but also iconic species for multispecies coral reef fisheries and temperate families with high ecological and socioeconomic importance in other regions (Appendix S6); thus, our results can be broadly applicable in a global context. We only included samples that were within the biogeographic extent of occurrence of each species group in the data set to down weight the importance of known biogeographic drivers (Appendix S7).

Biophysical and anthropogenic covariates

We extracted 19 biophysical and anthropogenic covariates that influence fish abundance and size structure at varying spatial scales. At macroecological scales (100s–1000s km), we sourced data on sea surface temperature (SST) and productivity (Fisher et al. 2010; Langlois et al. 2012). We initially sourced several climatological metrics (mean, minimum, maximum, and range) for each predictor: SST (degrees Celsius), nitrate (micromole per liter), phosphate (micromole per liter), and net primary productivity (NPP) (milligrams carbon per square meter per day). We retained the mean when they were highly collinear (r >0.8) (Appendix S8). We additionally sourced values of temperature at bottom depth from the CSIRO Atlas of Regional Seas (Huang et al. 2011) to test for potential mismatches with satellite‐derived measurements that might influence species distributions and body size throughout the water column (up to 50 m deep). Across Australia, these two metrics were strongly correlated (r = 0.99, p < 0.001) (Appendix S9), albeit some states displayed greater dispersion around the global relationship due to their higher vertical stratification (Langlois et al. 2012). We used geomorphological descriptors of seabed topography as proxies of within‐region (10s–100s of km) habitat availability and structure (Boström et al. 2011). Sample depth was measured in situ from the stereo‐BRUVs deployments. Covariates were sourced from the Australian marine physical environmental data set (Geoscience Australia, Huang et al. 2011) and the marine socioenvironmental covariates data set (Yeager et al. 2017).

We derived indirect potential proxies of human footprints from national and global socioeconomic products, including proximity to boat ramps, human population, human gravity (defined below), and management status (i.e., fished vs. NTMRs) (Appendix S10). Human population was derived from the LandScan 2011 population density grid by summing values around a 50‐km radius of each stereo‐BRUV deployment. Other human population metrics (i.e., census data from the Australian Bureau of Statistics) were considered, but the LandScan 2011 was the more accessible at the spatial scales considered that temporally matched the median year of the stereo‐BRUV data. We chose a 50‐km radius to capture the upper limit of small‐scale commercial and recreational fishing activities (Chuenpagdee et al. 2006). We acknowledge that there are other anthropogenic factors that may be associated with proximity to human population, such as water quality, habitat loss, vessel noise, and industrial development. Human gravity was computed, with a modification of the original metric proposed by Cinner et al. (2018), as the human population divided by the linear distance (kilometers) from the stereo‐BRUVs to the nearest coastal urban locality (from the Australia Bureau of Statistics significant urban areas data set; https://www.abs.gov.au/census). The prevailing type of management at the time of sampling, either fished (areas where local fisheries management regulations applied, e.g., size or catch limits) or NTMRs (areas closed to fishing and other extractive activities), was compiled from the Collaborative Australian Protected Area database (2020) and validated by regional managers.

Predictor variables were initially screened for the presence of potential outliers, collinearity, and relatively even spread of data points. Predictors with positively skewed distributions were log (x +1) or square‐root transformed, and we excluded predictors with correlations >0.8, resulting in a final set of nine predictor variables: relief, depth, SST mean, SST SD, nitrate, NPP mean, NPP SD, human gravity, and management status (Appendix S10).

Statistical modeling

We used generalized additive mixed models (Wood 2006) and a full‐subset model‐selection approach implemented in the R package FSSgam (Fisher et al. 2018) to test the relative influence of biophysical and anthropogenic covariates on the body‐size structure of fished species (assemblage‐level models, 82 species). Models were run in the gamm4 package (Wood & Scheipl 2014) with a binomial distribution with a logit link function for the probability of occurrence of large fishes and a negative binomial distribution with a log link function for the abundance of legal and sublegal fishes. The overdispersion parameter (tetha) of the negative binomial distribution was estimated via comparison of model fits in a range of potential values. The basis dimensions (k) of the cubic regression splines were limited to 4 to ensure unimodal relationships. The candidate set of models included combinations of all predictor variables, up to a maximum of 4 to avoid overfitting, and excluding any model with predictors correlated by >0.5 to avoid multicollinearity issues. Graham (2003) suggested a collinearity threshold of 0.28, but this is a rather conservative threshold that generally prevented SST and human gravity to occur in the same model, even when they were weakly correlated (Appendix S11). A null model containing the intercept and random effect terms was also included in the model set (Appendices S12 & S13). The random effect term contained the hierarchical nested structure of our data (1|state/clust.4 km + 1|Campaignid for assemblage‐level models, 1|ecoregion/clust.4 km + 1|Campaignid for regional species group models) (Appendices S12 & S13). This enabled us to control for large imbalances in sampling effort among states and small‐scale variations in stereo‐BRUV setup (e.g., sample separation, bait quantity, and soak time). We considered models in two units of the lowest AIC model as having substantial support (Burnham & Anderson 2003). When competing models were present, we selected the model with highest Akaike weights (ωAIC). The relative importance of each predictor was calculated by summing ωAIC of all models containing that variable (Fisher et al. 2018).

We used generalized linear mixed models (GLMMs) to explore the generality of the effect of indirect proxies of human footprints on the body‐size structure of regional fished species groups. Models were fitted with a binomial (i.e., probability of occurrence of large fishes) and negative binomial (i.e., abundance of legal and sublegal fishes) distribution in the glmmTMB R package (Brooks et al. 2017). Model‐selection was conducted in the MuMIn R package (Kamil 2018) and followed identical criteria as the full‐subset GAMMs for assemblage‐level models. Continuous predictors were standardized (mean 0 and SD 1) before analysis to account for differences in scaling. Model‐averaged coefficients were used to draw inferences about the magnitude and direction of the effect of proxies of human footprints, while accounting for uncertainty (weight of evidence, ωAIC) in the set of competing models. All models were visually inspected for violations of homoscedasticity by checking residuals against fitted values (Appendices S14–S16). We evaluated the suitability of the negative binomial distribution to model the large proportion of zeros in the count data by inspecting quantile‐quantile plots, running tests of zero inflation via the DHARMa package in R (Hartig 2017), and comparing the fits of models with a 0‐inflated negative binomial (intercept‐only model for the 0 inflation) (Appendix S17). We carried out the analyses in R 3.5.2 (R Core Team, 2019) with the ecocloud Platform (www.ecocloud.org.au), which is supported by the National Collaborative Research Infrastructure Strategy through the Australian Research Data Commons.

RESULTS

Assemblage‐level models

The most supported model to predict the probability of occurrence of large fishes contained the effects of human gravity, depth, and relief (i.e., a proxy of seabed complexity), although there was high uncertainty (ωAIC = 0.05–0.13) (Appendix S18), with models containing the effect of environmental predictors also receiving considerable support (Figure 2a). Similarly, the abundance of legal‐size fishes was best predicted by a model containing the effect of human gravity, depth, NTMRs, and relief; but, unlike models for large fishes, there was weak support for other competing models containing the effect of environmental predictors (ωAIC = 0.66) (Figure 2a; Appendix S18). For the abundance of sublegal fishes, the most supported model contained the effect of temperature (mean and variance) and NTMRs, but there was a high level of uncertainty: eight competing models had similar levels of support (ωAIC = 0.08–0.15) (Appendix S18). Mean SST was the strongest predictor for sublegal fishes; a model containing this term alone had support similar to other more complex models (ωAIC = 0.08 cf 0.11) (Appendix S18). Seascape relief had moderate levels of support for sublegal fishes, whereas human gravity had weak support (Figure 2a). The variance explained by the explanatory covariates in the most supported models was low (about 2–5%) (Appendix S18), but these models outperformed a null model containing the intercept and random effect structure.

FIGURE 2

Drivers and patterns of fished species occurrence and abundance across Australian neritic (1–50 m) continental shelf systems: (a) relative importance of explanatory variables to predict the probability of occurrence of large fishes (>90th percentile of body size distribution) and abundance of legal (≥ minimum legal size [MLS]) and sublegal (

Human gravity had a strong negative effect on the probability of occurrence of large fishes and the abundance of legal fishes. Large fishes were about three times more likely to be encountered in areas with low human gravity (Figure 2b), whereas areas with the lowest human gravity supported about five times the abundance of legal fishes (Figure 2f). Relief had a positive relationship with the probability of occurrence of large fishes and the abundance of legal fishes (Figures 2d & h); a ceiling effect was observed at high relief areas for the latter. Depth showed a unimodal (hump‐shaped) response with the probability of occurrence of large fishes and the abundance of legal fishes (peaks observed at approximately 30 and 20 m, respectively) before declining in deeper waters (Figures 2c & g). The abundance of sublegal fishes followed a biogeographic trend, increasing rapidly as mean SST increased and displaying a peak at intermediate levels of SST variance (Figures 2j & k). The presence of NTMRs had a positive effect on the abundance of fished species, but predicted the probability of occurrence of large fishes less well (Figure 2a). Predicted differences between NTMRs and fished areas were twice as large for the legal‐size fishes (0.44 MaxN vs. 0.25 MaxN for legal and sublegal, respectively), and there were larger uncertainties for the sublegal‐size fishes (Figures 2e & i).

Regional species group models

Models for regional species groups generally had low model weights and varying support for anthropogenic, environmental, and habitat predictors across large, legal, and sublegal fishes (Appendix S19). Despite the high uncertainty in the models and relatively greater importance of environmental and habitat predictors for most species groups (Appendix S20), the summed ωAIC for potential anthropogenic predictors (human gravity and NTMRs) was consistently greater for large and legal‐size fishes after controlling for habitat and environmental covariates (Figure 3a). Further, human gravity better explained the abundance of legal‐size fishes, but was of negligible importance for sublegal fishes in most species groups (Appendix S20). The variance explained by explanatory covariates in the most supported models varied strongly across species groups (0.04 cf 0.39) (Appendix S19), but was generally higher than the assemblage‐level models.

FIGURE 3

Size‐specific determinants of fished species occurrence and abundance across Australian neritic (1–50 m) continental shelf systems: relative importance (pooled Akaike weights [ωAIC ]) of (a) anthropogenic (human gravity [i.e., a proxy of human impacts] and management status), (b) environmental (nitrate, net primary productivity [NPP] mean, NPP standard deviation [NPP SD], sea surface temperature [SST] mean, and SST standard deviation [SST SD]), and (c) habitat (relief and depth) covariates to predict the probability of occurrence of large (>90th percentile of body size distribution) fishes and the abundance of legal (≥minimum legal size [MLS]) and sublegal (

Model‐averaged coefficients showed a generalizable, negative trend of human gravity on the probability of occurrence of large fishes and the abundance of legal‐size fishes. However, in the former, standard errors around model predictions overlapped 0 in three species groups considered (Notolabrus spp., Nemadactylus spp., and Lutjanus spp.) (Figure 4a). In contrast, there was no apparent trend in the abundance of sublegal fishes, and model‐averaged coefficients had standard errors overlapping 0. An exception to this was the positive relationship between human gravity and trends in the probability of occurrence of large fishes and abundance of both legal and sublegal pink snapper (Chrysophrys auratus). The NTMRs had a positive effect in four key fisheries species’ groups (C. auratus, Choerodon spp., Plectropomus spp., and Lethrinus spp.). Large fishes for these species groups generally displayed a strong effect of NTMRs in their probability of occurrence, although there was substantial variability around model‐averaged predictions and standard errors overlapped 0 in two species’ groups (C. auratus and Choerodon spp.). For abundance models, a consistent positive effect was found for the legal fishes, whereas the effect for sublegal fishes was variable (Figure 4b). Two species groups (Nemadactylus spp. and Lutjanus spp.) had a neutral relationship with NTMRs for both legal and sublegal fishes. The labrids Notolabrus spp. displayed negative associations with NTMRs.

FIGURE 4

Model‐averaged coefficients (SE) from generalized linear mixed models testing the effect of (a) human gravity (i.e., a proxy of human impacts) and (b) no‐take marine reserves (NTMRs) on the probability of occurrence of large fishes (dark squares), and the abundance of legal (gray triangles) and sublegal (light gray circles) fishes. Coefficients depict the magnitude of the effect while controlling for other environmental and habitat predictors and accounting for uncertainty of the models

DISCUSSION

Our large‐scale (approximately 10,000 km of coastline) quantitative analyses showed consistent support for the ubiquitous role of human footprints in shaping the body‐size structure of marine fished species by selectively removing large‐bodied individuals, an effect that was generally consistent across conspicuous species groups with contrasting habitat, environmental, and biogeographic affinities. For large‐bodied individuals (>90th percentile of the body‐size distribution) and individuals >MLS, probability of occurrence and relative abundance was best predicted by gradients of human gravity, depth, and seascape relief. The NTMRs were strong predictors of the abundance of legal‐size fishes, but had weak support for predicting the probability of occurrence of large fishes, highlighting potential limitations in NTMR coverage and design (Edgar et al. 2018; Goetze et al. 2021). In contrast, models for sublegal individuals were weak, but consistently featured positive relationships with environmental predictors, including mean and variance in SST. Combined, these results stress the importance of geographic areas subjected to low levels of human footprint in the surrounding seascape (Cinner et al. 2013, 2018) as spatial refuges for large‐bodied individuals across multiple ecosystem configurations. Our results also suggest that well‐designed NTMR networks (Knott et al. 2021, Goetze et al. 2021) and areas that remain relatively unaltered by humans provide useful reference points against which target levels of large‐bodied fish indicators can be assessed (Stuart‐Smith et al. 2017; McClanahan et al. 2019). This will boost the understanding of the potential benefits of protecting large‐bodied fishes for conserving trophic links (Dulvy et al. 2004; Mumby et al. 2006) and contribute to the recovery of local stocks through larval supply (Le Port et al. 2017; Andrello et al. 2017; Marshall et al. 2019).

Areas with the lowest human gravity were about three times more likely to support large‐bodied fishes and contained approximately five times the abundance of legal‐size fishes relative to the highest human gravity areas. Uncertainty around model predictions suggests that increasing distance from access points and population centers does not provide a panacea of protection (Devillers et al. 2015) and that there are other factors that can modulate the probability of occurrence and abundance of large‐bodied fishes across human‐gravity gradients. Because we used buffers of 50 km, designed to capture the footprint of small‐scale fishing activities (e.g., recreational fishing [Chuenpagdee et al. 2006]), as our metric of human gravity, other distal pressures, such as commercial fishing effort or stochastic and seasonal tourism to remote areas, may have influenced the probability of occurrence and abundance of large‐bodied fishes in the low human‐gravity regions (Amoroso et al. 2018).

There is a clear need for future research and methodological development of metrics that capture both the persistent and transient nature of human footprints (e.g., satellite‐derived measures of recreational and commercial boat traffic) across the seascape (Smallwood et al. 2011). Natural variations in the surrounding physical environment, such as local‐scale habitat features (Wilson et al. 2010; Nash et al. 2016) and regional‐scale productivity constraints (McClanahan et al. 2019), may also have contributed to the large degree of variation in the probability of occurrence and abundance of large‐bodied fishes. We set out to test the generality of patterns across tropical and temperate waters by synthesizing data from multiple independent studies, and the low variance explained by our models could partly be attributed to unbalanced sampling effort and stochasticity in sampling across independent studies synthesized. Besides the above factors, the inherent natural variability of fish abundance data, use of spatially static environmental predictors, and pooling of abundances across species with contrasting habitat preferences and life‐history strategies may have contributed to low variance explained.

In Australia, where the majority of the human population lives in the south of the continent, our human‐gravity metric could be confounded by gradients in SST and the natural biogeographic distributions of the endemic fished species (Langlois et al. 2012) (Appendix S11). Yet, we found that the probability of occurrence of large fishes and abundance of legal fishes from regional species groups with similar habitat affinities and life‐history strategies generally consistently declined as human gravity increased, after controlling for other habitat and environmental predictors. This result supports fishing as a primary driver beyond the species‐specific response to environmental and habitat factors (Cinner et al. 2018). A marked exception to this trend was the increasing abundance of both legal and sublegal pink snapper (C. auratus), a cosmopolitan fished species from temperate Australia that typically forms spawning aggregations in shallow embayments (Wakefield 2010). Around warm‐temperate Australia, human population centers typically co‐occur with these shallow embayments, and although data used in this analysis avoided pink snapper spawning seasons, their suspected tendency to aggregate near embayments may have confounded the effect of human gravity on their abundance distribution.

Networks of NTMRs are being expanded globally as a tool to slow marine biodiversity loss (Tittensor et al. 2014), but their potential fisheries benefits remain debated (Hilborn 2016; Di Lorenzo et al. 2020). Our assemblage‐level model demonstrated that the network of NTMRs around Australia can enhance ecological benefits by boosting the relative abundance of individuals above the MLS across a broad array of species and ecosystem configurations. The probability of occurrence of large fishes, however, displayed a weaker response to the effect of NTMR protection, a result that can be partly explained by the generally low number of individuals observed in the samples, as well as their potential larger home range size and mobility (Daly et al. 2021), which may not be captured adequately by small reserves (Claudet et al. 2010; Edgar et al. 2014). For instance, NTMRs had a strong effect on the probability of occurrence of large‐bodied individuals of coral trout (Plectropomus spp.), a species that displays high residency and territoriality and responds strongly to NTMR protection and direct catch or effort control measures (Hopf et al. 2016). The variable response across regional species groups highlights potential limitations in the design of NTMR networks around Australia where fisheries management measures are in place (Edgar et al. 2018; Goetze et al. 2021), such as gaps in NTMR coverage across the biogeographic extent of occurrence of species groups we analyzed (Roberts et al. 2021). Multiple intrinsic factors can also influence the biodiversity benefits of NTMRs, such as the size, age, level of compliance, connectivity of the zones, and the level of exploitation in the surrounding seascape (Claudet et al. 2008; Edgar et al. 2014; Cinner et al. 2018). For instance, in our synthesis, there was an uneven representation of the size of NTMRs across the human gravity gradient; small (<10 km2) NTMRs were overrepresented in areas of high human gravity, whereas the largest NTMRs (>100 km2) exclusively occurred in areas of low to intermediate human gravity. Gains in ecological and societal benefits may exist in establishing representative networks of NTMRs, including large NTMRs in areas with high human gravity (Cinner et al. 2018, 2020) that cover broad‐depth ranges (Goetze et al. 2021).

Habitat complexity (relief) was one of the strongest predictors of the probability of occurrence and abundance distribution of fished species. High topographic relief areas were more likely to contain large fishes and more legal and sublegal individuals. This highlights the importance of the local availability of high complexity substrates across large geographic areas to support fished species irrespective of their life history stage, likely through the provisioning of topographic refuges and availability of prey (Almany 2004). Our demonstration of the importance of habitat complexity, based on the available coarse resolution (250 m) national bathymetry data, suggests that existing data can usefully inform conservation and management of fished species at national scales (Costa et al. 2014). Yet, spatial management in regional contexts would be advanced by the acquisition of finer‐resolution seabed structural data, potentially improving model predictions (Pittman & Brown 2011).

Our results showed the widespread effect of human activities on the body size structure of marine fished species, an effect that was generally consistent across regional species groups with contrasting habitat, environmental, and biogeographic affinities. This demonstrates the importance of the use of consistent fisheries‐independent body‐size information at macroecological scales (Heather et al. 2020) to quantify human footprints on aspects of diversity that link to food‐web structure. The use of stereo‐BRUVs can overcome the depth limitations of diver‐based censuses, providing a cost‐effective means to expand the knowledge of species responses to anthropogenic impacts, including fisheries exploitation, across geographic space and depth profiles. In this context, our synthesis highlighted the importance of depth to explain patterns of occurrence and relative abundance of large‐bodied fishes, which is consistent with a longstanding paradigm of ontogenetic deepening (Audzijonyte & Pecl 2018), with the upper limit to depth distribution typically limited by the availability of habitat and prey (Galaiduk et al. 2018). Recently, this concept has been challenged by the finding that large‐bodied, older individuals of highly targeted species can preferentially occupy shallower depth strata following a fishing moratorium (Frank et al. 2018). The creation of comprehensive, adequate, and replicated networks of NTMRs (Goetze et al. 2021) and large‐scale stereo‐BRUV data sets covering broad depth ranges across different management scenarios (e.g., NTMRs vs. fished areas) and gradients of human pressure (e.g., human gravity here) will provide an experimental framework to help decouple the relative importance of these processes. This represents a key area for future research because the effect of ocean warming and fisheries could compress the three‐dimensional distribution of marine species and have large consequences for biodiversity and ecosystem productivity.

Appendix S1. List of species, by family, included in the Australian national synthesis of stereo‐BRUV length data.

Appendix S2. Summary of the Australian national synthesis of stereo‐BRUV length data, including sampling coverage across biogeographic “realms” and administrative states.

Appendix S3. Controlling for methodological differences in GAMMs and GLMMs models

Appendix S4. Retaining a temporal snapshot of the Australian synthesis of stereo‐BRUV data

Appendix S5. Temporal distribution of stereo‐BRUV deployments included in the analyses across (a, b) years and (c, d) seasons.

Appendix S6. Regional groups of conspicuous species around Australia included in the analyses.

Appendix S7. Biogeographic extent of occurrence of regional fished species groups investigated.

Appendix S8. Pearson's correlation coefficients between predictor variables initially considered to model the probability of occurrence and abundance of fished marine species across Australian neritic (1–50 m) continental shelf systems; color ramp indicating the strength and direction of the correlation (positive = blue; negative = red).

Appendix S9. Bivariate plot showing the correlation between values of sea surface temperature (°C) and temperature at bottom depth (°C).

Appendix S10. Hypothesis and associated predictor variables considered to model the probability of occurrence and abundance of fished marine species across Australian neritic (1–50 m) continental shelf systems.

Appendix S11. Pearson's correlation coefficients between human gravity and the other predictor variables considered to model the probability of occurrence and abundance of regional fished species groups across Australian neritic (1–50 m) continental shelf systems; color ramp indicating the strength and direction of the correlation (positive = blue; negative = red)

Appendix S12. Hierarchical nested structure of GAMMs and GLMMs models

Appendix S13. Comparison of model structures to control for spatial and temporal nonindependencies in the Australian synthesis of stereo‐BRUV length data.

Appendix S14. Residuals versus fitted values from “best” GAMMs for predicting the (a) probability of occurrence of large fishes, (b) abundance of legal fishes, and (c) abundance of sublegal fishes.

Appendix S15. Residuals versus fitted values from “best” GLMMs for predicting the (a, d, g, j) probability of occurrence of large fishes, (b, e, h, k) abundance of legal fishes, and (c, f, i, l) abundance of sublegal fishes across regional species groups.

Appendix S16. Residuals versus fitted values from “best” GLMMs for predicting the (a, d, g) probability of occurrence of large fishes, (b, e, h) abundance of legal fishes, and (c, f, i) abundance of sublegal fishes across regional species groups.

Appendix S17. Comparison of models using a negative binomial (Nb) distribution and a zero‐inflated negative binomial (ZinB).

Appendix S18. Best (within two units of the lowest AIC) GAMMs for predicting the probability of occurrence of large fishes, and the abundance of legal and sublegal fishes, across Australia neritic (1–50 m) continental shelf systems.

Appendix S19. Ranking of GLMMs for predicting the probability of occurrence of large fishes, and the abundance of legal and sublegal fishes, for regional groups of conspicuous species across Australia neritic (1–50 m) continental shelf systems.

Appendix S20. Relative importance of predictor variables to predict the probability of occurrence of large fishes, and abundance of legal and sublegal fishes, for regional groups of conspicuous species across Australian neritic (1–50 m) continental shelf systems.

Appendix S21. Global Moran's I results testing for autocorrelation in the raw data and model residuals of GAMMs and GLMMs for predicting the probability of occurrence of large fishes, and the abundance of legal and sublegal fishes, for assemblage‐level and regional species group models.

Click here for additional data file.