Using joint species distribution modelling to predict distributions of sea�oor taxa and identify vulnerable marine ecosystems in New Zealand waters

Effective ecosystem-based management of bottom-contacting �sheries requires understanding of how disturbances from �shing affect sea�oor fauna over a wide range of spatial and temporal scales. Using an extensive dataset of faunal abundances collected using a towed camera system, with spatially explicit predictor variables including bottom-trawl �shing effort, we developed spatial predictions of abundance for 67 taxa using Hierarchical Modelling of Species Communities. The model �t metrics varied by taxon: the mean ten-fold cross-validated AUC score was 0.70 ± 0.1 (standard deviation) for presence-absence and an R 2 of 0.11 ± 0.1 (standard deviation) for abundance models. Spatial predictions of probability of occurrence and abundance (individuals per km 2 ) varied by taxon, but there were key areas of overlap, with highest predicted taxon richness in areas of the continental shelf break and slope. The resulting joint predictions represent signi�cant advances on previous predictions because they are of abundance, allow the exploration of co-occurrence patterns and provide credible estimates of taxon richness (including for rare species that are often not included in community-level species distribution assessments). Habitat-forming taxa considered to be Vulnerable Marine Ecosystem (VME) indicators (those taxa that are physically or functionally fragile to anthropogenic impacts) were identi�ed in the dataset. Spatial estimates of likely VME distribution (as well as associated estimates of uncertainty) were predicted for the study area. Identifying areas most likely to represent a VME (rather than simply VME indicator taxa) provides much needed quantitative estimates of vulnerable habitats, and facilitates an evidence-based approach to managing potential impacts of bottom-trawling.


Introduction
Effective ecosystem-based management of bottom-contacting sheries requires understanding of how disturbances from shing affect sea oor fauna and habitats over a wide range of spatial and temporal scales (Clark et al., 2016;Pitcher et al., 2017).A key element needed to generate such understanding is knowledge of the spatial distributions of sea oor fauna.Reliable benthic spatial information is rare however, particularly for areas where deep-sea sheries occur and available data generally consist of point records of taxon presences assembled from disparate sources spanning many years or decades.Consequently, correlative modelling methods, in which statistical relationships between observed point occurrences of fauna and continuous environmental data layers are used to predict faunal occurrence across unsampled space (referred to as habitat-suitability, species-environment, or species-distribution models, e.g., Guisan and Zimmermann, 2000;Elith and Leathwick, 2009), are increasingly used to generate fullcoverage maps of either predicted habitat suitability or taxon presence for use in assessments of sea oor (benthic) impacts (e.g., Mazor et al., 2021).Spatially explicit models for sea oor fauna distributions are the cornerstone of marine conservation planning (Loiselle et al., 2003;Sundblad et al., 2011;Marshall et al., 2014;Por rio et al., 2014).
However, in most cases there are many species' distribution that are poorly described because there are too few records to generate robust species distribution models (SDMs) (Ellingsen et al., 2007).Consequently, the full complement of sea oor biodiversity is typically not represented in conservation planning, despite the important roles that less common species can play in the stability and functioning of marine ecosystems (Ellingsen et al., 2007;Zhang et al., 2018;2020a).
In the New Zealand region, SDMs have been used for over two decades to provide predictions of sea oor faunal Until recently, all these models (other than one small-scale study, Rowden et al., 2017) have been informed by faunal occurrence data compiled from research trawl bycatch and scienti c museum/institute records, which offer information about the presence of a taxon at any given site but not its abundance (density) or absence (Bowden et al., 2021b).Such models can only yield predictions of relative habitat suitability (the likely distribution of species), rather than predictions of expected abundance (Elith and Leathwick, 2009).These existing models ful l their purpose in that they provide the best estimates of the distributions of sea oor fauna in an environment that remains datalimited in terms of knowledge about both faunal distributions and the physical characteristics of their habitats.However, knowledge about spatial variations in species' abundances is crucial for understanding ecosystem functioning (Rullens et al., 2022); for instance, the presence of a single bryozoan colony at a site will not have the same ecological in uence, or conservation value, as a high density of bryozoan thickets (Wood et al., 2013).A number of studies have assessed whether presence-absence models can be used as surrogates for abundance distributions, with contrasting results; there is some evidence that the correlation between probability of occurrence and density is weaker for species with broader ecological niches than those with narrow niches (e.g., as in Rullens et al., 2021 and references therein).
While the value of single-taxon SDMs is certainly recognised, so too are their limitations (Loiselle et al., 2003;Bowden et al., 2021b; Lee-Yaw et al., 2022), which include their disconnect from most ecological theory (Zhang et al., 2020a;Stephenson et al., 2022a).Combining or stacking SDM predictions for different species is often used to examine community-level species distribution patterns (Guisan and Rahbek, 2011;Calabrese et al., 2014;Zhang et al., 2019).However, as these stacked or combination approaches model each taxon separately and combine spatial distributions post-hoc, species predictions are independent and do not allow for inter-species interactions and their combined relationships with the environment (e.g., Compton et al., 2013;Zhang et al., 2020a).In addition, rare species are often not included in community-level species distribution patterns because there are too few data to generate quantitative SDMs for these species (Zhang et al., 2020a) (but mechanistic models can be used for rarer species, e.g., Stephenson et al., 2020).
Joint species distribution models (JSDMs) offer a solution to these limitations (Warton et al., 2015;Ovaskainen et al., 2016a;Ovaskainen et al., 2016b;Zhang et al., 2018;Tikhonov et al., 2020).Where SDMs link abiotic (environmental) covariates to species occurrences, JSDMs make use of latent variable approaches as well (i.e., variables which are indirectly inferred within the model itself), describing occurrence as a function of environment but incorporating biotic associations such as competition, predation, or parasitism (Pichler and Hartig, 2021).The latent variable approach makes use of residual correlation between species, where positive residual correlation is assumed to represent co-occurrence, and negative residual correlation means the species co-occur less often than expected, given the environment (Zurell et al., 2018).While unobserved biotic interactions are the target in the spatial latent variable approach, it is important to consider that unknown or unmeasured abiotic contributions could also be having an effect on distribution patterns (Blanchet et al., 2020;Pichler and Hartig, 2021;Poggiato et al., 2021).
Since the early 2000s there has been a concerted effort in New Zealand to collect quantitative information of sea oor taxa abundances using non-destructive methods, i.e., using underwater cameras (e.g., Rowden et al., 2002;Clark and Rowden, 2009).Since 2006 imagery data has been collected using the Deep Towed Imaging System (DTIS) (Hill, 2009;Bowden and Jones, 2016).These data provide the opportunity to generate spatial estimates of taxon abundances (rather than simply occurrences) and, in this case, to also test the use of JSDMs which may better account for species-interactions and more easily incorporate rarer species in a quantitative manner (Zhang et al., 2020a).Furthermore, outputs from JDSMs can be used to predict the occurrence of ecosystems represented by a composite of species (Ovaskainen et al., 2017), including Vulnerable Marine Ecosystems (VMEs) (Gros et al., 2022).
VMEs are ecosystems (comprised of species groups, communities or habitats) that are physically or functionally vulnerable to anthropogenic impacts; the most vulnerable ecosystems are those that are both easily disturbed and very slow to recover or may never recover.Examples of VMEs include cold-water coral reefs, sponge beds, and deepsea hydrothermal vents (Roberts et al., 2009).VMEs are important because they may be unique, comprise structurally complex species that support and/or provide essential habitat for other species (including sh), and serve as natural carbon sinks (Cathalot et al., 2015;De Froe et al., 2019).However, VMEs are also threatened by human activities such as bottom trawling, oil and gas exploration, and deep-sea mining (Clark et al., 2016).Although VME is a term primarily associated with application in Areas Beyond National Jurisdiction (ABNJ), from an ecological and conceptual point of view, the concept of VME extends to all sea oor habitats, including areas within national jurisdictions.
Efforts to conserve VMEs are critical for maintaining healthy and resilient marine ecosystems (Van Dover, 2010).Management measures can include the designation of marine protected areas and implementing regulations to reduce the impacts of human activities on VMEs (e.g., Brodie and Clark, 2003;Australia and New Zealand, 2020).One approach to VME management involves identifying VME indicator taxa (e.g., Parker and Bowden, 2010), and then calculating vulnerability indices to assess the susceptibility of VMEs to various stressors, including shing activities (Clark et al., 2016;Morato et al., 2018;Gros et al., 2023).If benthic taxa data exist, the vulnerability of areas based on the occurrence or abundance of VME indicator taxa can be used to highlight particular areas which may represent a VME (Gros et al., 2023).This approach can help policymakers and stakeholders prioritize areas for protection and management and make more informed decisions about the sustainable use of marine resources.
Using an extensive sea oor imagery dataset providing taxon occurrence and abundance estimates across a wide spatial area (Bowden et

Study area
The study area spans a broad segment of New Zealand's marine environment, encompassing two major topographic features; Chatham Rise and Campbell Plateau, and the intervening areas of the Otago continental shelf and slope, and Bounty Trough (Fig. 1).Chatham Rise is a continental rise extending eastwards from the South Island of New Zealand for approximately 1000 km.The crest of the Rise is at depths of 300-400 m, with shallower features at its western end rising to approximately 50 m depth in places, and the Chatham Islands mark its eastern end.The Sub-Tropical Front, which is a major oceanographic feature in the area, coincides with and is partially constrained by Chatham Rise and, in consequence, intense phytoplankton blooms occur along its length (Nodder and Northcote, 2001;Stevens et al., 2019).The Chatham Rise is the most biologically productive large-scale, deep-sea sheries region in New Zealand's EEZ (Clark et al., 2000;Marchal et al., 2009).Campbell Plateau is a broad submarine plateau extending to the south and southeast of New Zealand's South Island.Much of the plateau lies at water depths of 500-1000 m but with shoal areas rising to the surface around islands in the west, and to less than 150 m on Pukaki Rise in the east (Fig. 1).On its western and southern boundaries, the plateau descends steeply to depths greater than 3000 m.The Sub-Antarctic Front, which forms at the northern boundary of the east-going Antarctic Circumpolar Current is constrained by the southern edge of the plateau.The Sub-Tropical Front lies to the north of the plateau, where warmer subtropical waters entrained down the west coast of the South Island recurve northwards, transporting nutrient rich, relatively colder waters towards the southern ank of Chatham Rise (Murphy et al., 2001;Nelson and Cooke, 2001;Hayward et al., 2007;Hurlburt et al., 2008;Mackay et al., 2014).Biological productivity, in terms of water-column primary production, is lower than on Chatham Rise but is elevated above levels in the surrounding ocean (Gutierrez-Rodriguez et al., 2020) Commercially important bottom trawl sheries exploiting populations of scampi (Metanephrops challengeri), hoki (Macruronus novaezelandiae), orange roughy (Hoplostethus atlanticus), oreos (Pseudocyttus maculatus, Neocyttus rhomboidalis), arrow squid (Nototodarus spp.), jack mackerels (Trachurus spp.), southern blue whiting (Micromesistius australis -although largely midwater trawl) and others have operated across the area in depths down to approximately 1200 m since the 1970s (Fisheries New Zealand, 2023).Recent summaries of bottomcontacting trawl history (Baird and Mules, 2021) show highest trawling intensity, primarily from the hoki shery, at 450-700 m depth west of Mernoo Bank and on the southern and northern central anks of Chatham Rise.

Faunal data
The faunal dataset used is an augmented version of that developed by Bowden et al. (2019) and Anderson et al. (2020) for application as an independent test of existing SDM predictions (Bowden et al., 2021b).In brief, these data were derived from analyses of sea oor video transects collected with the DTIS towed camera system during seven surveys spanning the period 2006 to 2020.Data consist of quantitative abundance estimates for invertebrate taxa at 467 sea oor sites: 358 from Chatham Rise and 109 from Campbell Plateau (Fig. 1).The survey data span a depth range from 49 to 1813 m but because there were few sites at the shallow and deep extremes of this range, model predictions here were restricted to depths from 100 to 1500 m (and to the spatial extent of the features, Fig. 1).This depth range encompassed 449 sites (96% of the total available).
More than 380 individual taxa were identi ed from these surveys but many of these were either operational taxonomic units (OTUs) or species-level names that were often not consistently recorded within and between surveys.In an extensive audit process (see Bowden et al., 2019), the full dataset was aggregated to a set of 139 taxa, spanning a range of taxonomic levels, which were reliably and consistently recorded across all studies and years (Table S1, supplementary materials).A subset of 67 taxa was then extracted from the aggregated taxon list by selecting only those that were observed at ve or more survey sites -a minimum number of samples required for the JSDM, see next section for details (Table 1).
Table 1.Seafloor invertebrate taxa (operational taxonomic unit [OTU] and common name) for which abundance data (as individuals per 1000 m 2 ) were used in development of Joint Species Distribution Models (JSDMs).Orange cells represent VME indicator taxa recognised by the South Pacific Regional Fisheries Management Organisation (SPRFMO) and used in this study (see Methods for why some taxa that are apparently VME indicator taxa were not included in this study; shown a grey cells).The number of individual OTU encompassed by each taxon (N) and the number of sites at which each taxon occurred (Site) are provided for context.The order Gorgonacea is no longer accepted as a taxonomic group, but is used to here represent the VME indicator taxon used by SPRFMO to identify "sea fans"; suborders Holaxonia, Calcaxonia and Scleraxonia of the Octocorallia.
Pennatulacea is used here to represent those taxa in the Order Pennatuloidea.

Explanatory variables
Explanatory variables available for the JSDMs were those selected in previous species distribution models of invertebrate fauna developed for the region and were available at 1 km grid resolution (Anderson et al., 2016;Georgian et al., 2019;Stephenson et al., 2022b) (Table 2).The 16 variables selected for consideration in this study describe sea oor characteristics, water chemistry at the sea oor, water physics, productivity, and bottom trawl shing distribution and effort (swept area per 5 x 5 km grid cell for the period 1989 to 2016) (Table 2).

Joint Species Distribution Modelling
The faunal data were analysed using the JSDM method Hierarchical Modelling of Species Communities (HMSC), which yields inference both at single-taxon and community levels from a single modelling process (Ovaskainen et al., 2017).HMSC partitions variation in species occurrences or abundances into components that can be interpreted in relation to environmental ltering, species interactions, and random processes.That is, in addition to modelling the relationship between biological data and explanatory variables, the HMSC framework can incorporate information on inter-species interactions and the spatial context of the sampling design.
With most SDM methods, the inclusion of many explanatory variables should be avoided because they generally only provide minimal improvement in predictive accuracy, complicate interpretation of model outcomes and can lead to over-tting and loss of generality (Elith et al., 2006).In order to produce a parsimonious and generally-applicable model (i.e., a model with appropriate explanatory or prediction power with as few predictor variables as possible), we rstly removed explanatory variables which had high co-linearity (i.e., ≥ 0.7, as per Dormann et al. ( 2013)); the decision of which of the two co-linear variables to retain was based on expert ecological knowledge of the likely in uence on faunal distributions.Following this process, variables standard deviation of depth, Bathymetric Position Index -broad, dissolved oxygen and particulate organic carbon export to seabed were removed.A further re nement of predictor variables was undertaken based on convergence scores (i.e., based on whether relationships between explanatory variables and faunal data were observed; a model may not converge if this is not the case) and the exploration of model ts (principally focusing on increasing model ts for taxa with low model t scores) in preliminary JSDM models where at least one variable from each of the 4 explanatory variable categories was used.
Seven explanatory variables were selected based on exploration of model ts and expert knowledge for use in nal JSDMs (Table 2): depth (bathy); tidal current speed (tidcurr); salinity (salinity); temperature residuals (tempres); pro le curvature (profcurv); Eppley Primary Productivity (epp); and trawling history (trawl).This set of variables included measures associated with food supply (epp, tidcurr), depth (bathy), water mass (salinity and tempres), sea oor topography (profcurv), and trawling history (trawl), all of which were considered likely to in uence distributions of fauna and community assembly processes.
HMSC from the hmsc R package (Tikhonov et al., 2020) was used to t JSDMs to the sea oor data simultaneously combining information on environmental covariates, random spatial effects and inter-taxa interactions in a single model.A two-part hurdle model was used to model distribution of abundance.In this procedure, a binomial model with a probit distribution was used initially to predict the probability of occurrence, followed by a separate model with a Gaussian distribution to estimate (log(x + 1) transformed) abundance for locations where presence was predicted (referred hereafter as abundance conditional on presence).The outcomes of the presence-absence and abundance conditional on presence models were multiplied together (hurdled) to create an overall estimate of abundance.We refer to these nal (hurdled) outputs as abundance.In addition to the explanatory variables (environmental data and trawling history) a random spatial effect (that also models co-occurrence among species) was included at the level of sampling station, using a latent factor approach (Ovaskainen et al., 2016a).The explained variation among the xed and random effects was partitioned using methods described by (Ovaskainen et al., 2017).The model was tted to the data with Bayesian inference, using the posterior sampling scheme of Tikhonov et al. (2020).Four chains, each with 250 samples and thinning parameter of 10 were used, which ensured that there was convergence of the Markov chain Monte Carlo simulations (i.e., following best practice recommendations from Ovaskainen and Abrego (2020)).
Model t metrics included AUC (area under the curve of the receiver operating characteristic) and R 2 (coe cient of determination) for each taxon for the presence-absence and abundance conditional on presence models, respectively.Model t metrics were calculated using withheld data (predictive power) by performing a ten-fold crossvalidation.AUC is an effective measure of model performance and a threshold-independent measure of accuracy (Allouche et al., 2006).AUC scores range from 0 to 1, with a score of 0.5 indicating model performance is equal to random chance.R 2 measures the proportion of the variance in the dependent variable that is predictable from the independent variables.Overall model performance was assessed by calculating the mean model ts across all taxa for both presence-absence and abundance models.
Finally, the parameter estimates were explored, and spatial predictions and associated uncertainty (standard deviation of the mean) were made for the study area using the 'constructGradient' and 'predict' functions in the hmsc R package.Taxon richness was estimated by summing the individual taxa occurrence predictions (Calabrese et al., 2014).

Mapping VME indices
Adapting the approach developed by Gros et al. (2023), the spatial predictions of abundance from the JSDMs were used to map areas most likely to contain VMEs.This approach consisted of: 1) assessing the vulnerability to bottom shing of all taxa within our dataset that were identi ed as VME indicator taxa; 2) calculating and mapping an "abundance-based VME index", using the cumulative abundance of VME indicator taxa weighted by a VME vulnerability score; 3) calculating and mapping a "richness-based VME index", using the richness of VME indicator taxa weighted by the VME vulnerability score; 4) mapping a "con dence index" which estimates the con dence associated with the spatial predictions of the VME indicator taxa; and 5) identifying the most likely areas in which VMEs might occur, based on the overlap of the highest abundance-based and richness-based VME index scores in areas with the highest model prediction certainty.All R code used for JSDM modelling and mapping VME indices is available on Github (see data availability).

Assigning vulnerability scores to VME indicator taxa
In the rst instance, VME indicator taxa were identi ed from the 67 taxa used in the JSDM models, following de nitions from the Conservation and Management Measure for the Management of Bottom Fishing in the High Seas of the South Paci c (SPRFMO, 2022  S1, in Supplementary Materials 1).
The set of VME indicator taxa identi ed in SPRFMO (2022) includes some high-level taxonomic groupings that encompass morphotypes in our image-based dataset that have different vulnerabilities to shing.Therefore, we reviewed the source observation data and excluded morphotaxa that were small in height, primarily associated with soft sediment, or were recorded at high densities in areas of intensive trawling (emphasising the "Vulnerable" element for this assessment).Five taxa recorded in the image data that would nominally be grouped as VME indicators were excluded from the VME analysis.(1) The grouping 'Anemones' included some large, erect morphotypes that were associated with seamounts and rocky reefs (high vulnerability) but also many small, lowpro le morphotypes that occurred in high densities on soft sediments (low vulnerability).( 2) Because identi cation of anemones in the dataset was generally at coarse taxonomic level it was not possible to separate the taxa reliably based on vulnerability, so the entire grouping Anemones, and the morphologically similar group Corallimorpharia, were excluded.(3) A distinctive pennatulacean taxon recorded as 'Kophobelemnon spp.' was only a few centimetres high and was recorded in very high densities on soft sediments in areas subject to intensive disturbance from trawl sheries (low vulnerability).( 4) Similarly, the solitary soft coral Taiaroa tauhou is small in size, can retract into the sediment, and occurred at high densities in areas of high shing disturbance.(5) Finally, a hexactinellid sponge morphotype recorded as Cladorhiszidae was excluded as this is also small in stature and occurred only on soft sediments.

VME indices
For the abundance-based VME index, spatial predictions of abundance for each VME indicator taxon were rst multiplied by their vulnerability scores, so that areas with high predicted densities would have higher vulnerability scores than those with low densities.Spatial estimates of taxon-speci c vulnerability scores were then summed across all VME indicator taxa to produce an overall abundance-based VME index (Gros et al., 2023).For the richnessbased VME index, spatial estimates of each taxon's probability of occurrence were multiplied by its vulnerability to bottom shing score and then summed, to produce a richness-based index of vulnerability (Gros et al., 2023).Despite the similarity in these two VME indices, abundance and richness represent distinct and important measures by which VMEs can be identi ed (Lockhart and Hocevar, 2021).For example, in the High Seas of the South Paci c, an encounter with a potential VME can be judged both on the weight of bycatch of certain key VME indicator taxa (i.e., presumably representing areas of high abundance on the sea oor) or based on a lower weight of bycatch of three or more VME indicator taxa (i.e., representing areas of high diversity on the sea oor) (Geange et al., 2020; SPRFMO, 2022).

Con dence index
Gros et al. (2023) developed a con dence index based on the quality of imagery used, but in our dataset all imagery was from the same camera system and thus of essentially the same quality (following aggregation of fauna to broader OTUs).We develop a different con dence index, here, by multiplying the spatially explicit uncertainty estimates (standard deviation of the mean abundance) for each VME indicator taxon by its associated VME vulnerability score.Uncertainty estimates for all VME indicator taxa were then summed to give an overall uncertainty index which was then divided by the abundance-based VME index.This approach is akin to calculating the coe cient of variation (CV, also sometimes known as relative standard deviation) of the abundance-based VME index and therefore represents a relative estimate of uncertainty.Here, we use a conservative threshold of 0.33, with a CV > 0.33 indicating moderate to high uncertainty, and therefore focus on areas of VME indices with CV < 0.33 (low uncertainty).

Identifying potential VME locations
In the absence of quantitative thresholds for de ning VMEs, areas with the highest abundance-based and richnessbased VME indices (subjectively de ned as the 95th percentile) were used to identify likely locations of VMEs.Similarly to Gros et al. (2023), the areas of high abundance-based VME index and high richness-based VME index were mapped, along with the con dence index.The subjectively de ned 95% percentile threshold aligns with examples of other studies predicting "primary habitat" or "core habitat", i.e., areas which are most likely to contain the study taxa (e.g., Tong et al., 2013;Anderson et al., 2022;Stephenson et al., 2023b).The overlapping of areas with highest VME index scores (and high model con dence) are more likely to represent locations of VME, or areas of priority for future investigation / con rmation of VME presence.

Model ts
Model ts for presence-absence and abundance models varied by taxon (Table 3).Model ts were reasonable, with a mean cross-validated AUC score of 0.70, but a lower mean cross-validated R 2 of 0.11 (these values were near identical when calculated for VME and non-VME indicator taxa).The family Gorgonocephalidae (basket stars) had the lowest cross-validated AUC score (0.45), suggesting predictions for this taxon were no better than would be expected by random chance, whereas the species Enypniastes eximia (swimming sea cucumber) had the highest cross validated AUC score (0.94).Cladorhizidae (family of demosponges) had the lowest cross-validated R 2 of -0.26, (no predictive power), whereas the subclass Echiura (spoon worms) had the highest R 2 of 0.49 (Table 3).Of the 67 taxa included in the JSDM models, 17 and 21 taxa were considered to have poor performance for presence-absence and abundance models, respectively (noting that these poorly performing models are predominantly for non-VME indicator taxa used in this study, Table 3).

Parameter estimates
There was strong evidence (95% posterior probability, or higher) of both positive and negative relationships between explanatory variables and taxon occurrence and abundance (see mean posterior values and support in Supplementary Materials 2).Relationships with environmental explanatory variables varied by taxon (see the variance explained by each explanatory variable and by the random spatial effect in Supplementary Materials 2) but there were general trends for both taxon occurrence and abundance to be negatively related to increasing depth (bathy) and decreasing salinity (salinity), and positively related to increasing pro le curvature of the seabed (profcurv) and primary productivity (epp_mean).Trawl history (trawl) was negatively associated with the occurrence and abundance of many VME indicator taxa (e.g., Goniocorella dumosa, Demospongiae, Antipatharia, Scleractinia, Hydrozoa) but in some cases positively associated with occurrence of non-VME indicator taxa (e.g., Echinothurioida, Gastropoda, Buccinidae) (Supplementary Materials 2).
Two broad groups of pairwise taxon associations were observed when using presence-absence data, primarily split between taxa found in soft-sediment and on hard substrate (Supplementary Materials 2).Primarily soft sediment taxa such as Holothuroidea, M. challengeri, Anemones, and Pennatulacea (among others) had a strong positive association with each other (i.e., co-occurrence when accounting for the environmental niche) and a negative association with other taxa associated with hard substrates such as Goniocorella dumosa, Scleractinia (Enallopsammia spp., Madrepora spp., Solenosmillia spp.), Antipatharia and Demospongiae (amongst others).There was also strong support for positive association among a large number of other taxa which may represent several distinct communities, but which are not obviously apparent in the taxon-to-taxon association matrices (for details see Supplementary Materials 2).

Spatial predictions of taxon richness and abundance
Spatial predictions of probability of occurrence, abundance (individuals per km 2 ), and associated uncertainty (standard deviation) varied by taxon (see Supplementary Materials 3 for all predictions) but there were key areas of overlap in the distribution of taxa, which are re ected in the predicted patterns of taxon richness.The highest predicted taxon richness occurs in depths from 400-800 m on the northern and southern anks on the western part of Chatham Rise and adjacent areas of the continental shelf break and slope (Fig. 2).There were also discrete (relatively small) areas of high richness areas of the continental shelf break and slope across Campbell Plateau and just off the coast of Fiordland (western part of the study area -Fig.2).
Predicted spatial patterns of abundance for Goniocorella dumosa (Fig. 3) and Hexactinellid sponges (Fig. 4) are presented here as example taxa because they are considered key habitat-forming VME indicator taxa in the study area (and based on observations from the DTIS imagery) (Tracey et al., 2011).Many of the habitat-forming taxa associated with VMEs had restricted distributions, with low to moderate probability of occurrence and abundance over most of the study area and concentrated areas of high probability of occurrence and abundance, primarily along steep slopes and on seamounts (e.g., Fig. 3 and Fig. 4, but also see Supplementary Materials 2).
The highest abundances for Goniocorella dumosa were predicted on seamounts and areas of steep slopes across the study area (e.g., North East Solander Trough, North West Bounty Trough, Graveyard Seamounts (centre-north Chatham Rise), South East Chatham Rise, South West Chatham Rise, Fig. 3), with moderate abundance predicted across central parts of Chatham Rise and on deeper parts of the continental shelf, especially in northern parts of the study area (Fig. 3).Abundances were predicted to be low in most parts of the study area (Fig. 3).Predictions of Hexactinellid sponge abundance were low across most of the study area and highest on seamounts and areas of steep slopes (Fig. 4).However, in contrast to Goniocorella dumosa, areas of moderate abundance were restricted to very deep parts of the study area (Fig. 4), noting these areas also have the highest model uncertainty (e.g., see Supplementary Materials 2).

Spatial predictions of VME indices and potential VME locations
Values of both the abundance-based and the richness-based VME indices were predicted to be low across most of the study area (approx.1.17 million km 2 , yellow in Fig. 5), with the threshold upper values (95th percentile) having restricted distributions that differed between the two indices (Fig. 5, and see detail of abundance-based and the richness-based VME indices in Figure S2, in Supplementary Materials 1).The abundance-based VME index was highest on the shelf break, shallower banks, and around offshore islands, as well as on seamount complexes throughout the study area (approx.59,000 km 2 , blue in Fig. 5).The richness-based VME index was highest on slopes, particularly in the deeper waters on the northern and southern anks of Chatham Rise (approx.59,000 km 2 , green in Fig. 5).
There were numerous areas of small spatial extent where the highest predicted values of both abundance-based and the richness-based VME indices coincided, and thus where VMEs are most likely to occur (approx.5500 km 2 , or 0.4% of the study area, red in Fig. 5, and see inset maps).Most of the areas where the highest predicted values of the two VME indices coincided were also in areas with low model uncertainty (i.e., no hashing in Fig. 5) providing greater con dence that these areas represent the most likely areas where VMEs mayoccur.

Discussion
Using a spatially extensive dataset of faunal abundance in areas around New Zealand we developed spatial predictions of distribution for 67 taxa using recently developed JSDMs.The resulting predictions are signi cant advances on previous predictions because they are joint predictions (i.e., each taxa is predicted simultaneously rather than sequentially, Wilkinson et al., 2021), therefore allowing the exploration of co-occurrence patterns and providing credible estimates of taxon richness (Ovaskainen et al., 2016a;Zhang et al., 2019).Most importantly, the spatial predictions of abundance for a large number of taxa are an advance on previous work in the region, which has either consisted of spatial predictions of abundance for a limited number of taxa across relatively small areas (e.g., Rowden et al., 2020) or are predictions of probability of occurrence, or habitat suitability (Stephenson et al., 2021;Stephenson et al., 2023a).Probability of occurrence or habitat suitability only indicates the likelihood of a taxon occurring in a particular area but does not necessarily provide information on relative abundance; in many cases, the highest abundances are found in smaller, more localised areas than the wider occurrence predictions (e.g., Rullens et al., 2021).This lack of direct information for abundance has important implications for conservation planning because the more extensive area identi ed from the occurrence models may not have the highest abundances and therefore may not be the best areas for conservation (Stephenson et al., 2022a).In addition, in the case of habitat formers (e.g., as several of the VME indicator taxa are in this study) if abundance information is not used it means there is no spatial information about the locations of what might be classi ed as actual VMEs.The presence of a coral VME indicator taxa does not indicate that there is a coral reef VME at that location (Howell et al., 2011;Rowden et al., 2017).Here we provide abundance estimates for 67 taxa, including for several habitat-forming taxa considered VME indicator taxa.We also, for the rst time, explore how these abundance data can be directly related to one or more of the FAO (2009) functional de nitions of a VME (e.g., as implemented by Gros et al., 2023) to highlight important areas which are most likely to contain VMEs and may be at most risk from the impact of bottom trawling.
Information of this nature is of critical importance for effective spatial management that aims to prevent or mitigate signi cant adverse impacts to VMEs (Gros et al., 2022).

Critical appraisal of JSDM
Given that JSDMs can explain spatial variation in community composition via inclusion of spatial latent variables, they show promise as an improvement to generic niche models for a variety of applications and have been shown to have some of the highest predictive power compared to other commonly used SDM approaches (Norberg et al., 2019;Zhang et al., 2020b).Inclusion of inter-species associations means that aspects of community ecological theory are implicit in JSDM predictions to inform spatial planning to identify priority areas for protection.For example, we would expect that species with mutualistic associations would have positive residual correlation, given their co-occurrence (Zurell et al., 2018).Therefore, priority areas for protection inferred from JSDM spatial predictions should better identify overlapping areas for community-level protection (via positive residual correlations), for both conservation targets and the interspecies relationships they depend on.In contrast, negative residual correlations between species will infer where differences lie.The ability to consider these interspecies relationships (differences and similarities) will empower conservation planners and reduce the likelihood of conservation con icts (Inoue et al., 2017).
JSDMs also provide the opportunity to perform conditional joint predictions (i.e., predicting species' occurrence or abundance given the known occurrence or abundance of other species, Ovaskainen and Abrego, 2020; Zhang et al., 2020a).This approach was not undertaken in the present study but may provide gains in predictive power, particularly for rare taxa, and may be especially useful in accounting for (assumed) biotic interactions from the cooccurrence matrix output from JSDM ( The effect of anthropogenic stressors on species' distributions are rarely accounted for in species distribution modelling (Elith and Leathwick, 2009), yet predictions that incorporate historic and current impacts are likely to represent more realistic (reduced or expanded) estimates of current distribution (Bowden et al., 2021a).For example, the erect branching coral Goniocorella dumosa was negatively associated with shing effort and therefore the predicted distribution or densities were reduced compared to a prediction using environmental estimates alone (the latter was not presented in this study).In contrast, the predatory scavenging whelk taxon Buccinidae was positively associated with shing effort distributions, with densities predicted to be higher compared to predictions using environmental estimates alone.Given the potential for altered taxa distributions resulting from the inclusion of anthropogenic impacts, it is important to use impact-adjusted species distribution or habitat suitability layers in spatial planning processes (Moilanen et al., 2011;Stephenson et al., 2023b).Taking this approach will avoid the possibility of designing ineffective conservation measures by protecting areas which previously had high abundances, but which may no longer be of high conservation value due to an anthropogenic impact (e.g., in the case of bottom trawling distribution and effort; Rowden et al., 2019).Nevertheless, previously high abundance but impacted areas could also be considered for their value as areas which may support recovery (Baco et al., 2019); these areas should especially be considered in cases where there are few pristine areas within a species range.
Bottom trawl sheries are sometimes focused where benthic taxa are particularly abundant or where VMEs are present (e.g., seamounts) and the impact of bottom sheries can be profound for these seabed ecosystems (FAO, 2009; Clark et al., 2016; e.g., Baco et al., 2020;Goode et al., 2020).In line with this issue, many nation states have recognised the threat posed to VMEs by bottom-contact trawling and have sought to identify and protect these ecosystems (Morato et al., 2010), including New Zealand.Here, we identify areas that are most likely to contain VMEs in part of the New Zealand EEZ where bottom trawling occurs (acknowledging that these areas are what is predicted to remain since bottom trawling distribution and effort are incorporated as a predictor in the models for VME indicator taxa).Given these already potentially reduced distributions, the estimates presented here of where and how much potential VME area occurs could be used as relevant and up-to-date information for assessing current marine protected areas, which were designed in part to protect VMEs.For example, the prediction maps could be used to inform modi cation of the boundaries of current Seamount Closure Areas and Benthic Protection Areas  2023) noted that the determination and categorisation of the VME indices into two relative (not absolute) levels (low versus high) was not ideal, and that future studies should "determine and validate how to reliably identify priority sites from the gradient of VME index values".Therefore, we recommend that future sampling and analyses of underwater camera imagery from the New Zealand EEZ undertake further exploration of the usefulness of the approach for identifying VMEs demonstrated here.In particular, to identify ecologically justi able thresholds for identifying VMEs (rather than the arbitrary threshold used here) for conservation and management purposes.
Following development of JSDMs for benthic taxa for the entire New Zealand EEZ, and on the completion of the further sampling and analyses suggested above, it may be possible to quantitatively incorporate robust predictions of VMEs into future spatial planning efforts in the region.Ideally, this planning would also include other anthropogenic threats to the biodiversity supported by VMEs, including climate change effects (e.g., as in Anderson et al., 2022).Achieving these steps would provide a wealth of additional information on sea oor communities and habitats for use in future spatial planning efforts.

Declarations Funding
distributions in the deep sea to inform research, environmental management, and prediction of climate change effects across a range of spatial scales (see examples and references in Tracey et al., 2011; Wood et al., 2013; Anderson et al., 2016; Rowden et al., 2017; Stephenson et al., 2021; Anderson et al., 2022; Stephenson et al., 2023a).
al., 2019; Anderson et al., 2020) and the recently developed JSDM approach (Ovaskainen & Soininen 2011, Warton et al. 2015, Ovaskainen et al. 2016), we predict the distribution and densities of 67 invertebrate taxa (including 26 VME indicator taxa) and estimates of taxonomic richness.We then use the recently developed concept of VME indices (Morato et al., 2018; Gros et al., 2023) to generate spatial estimates of likely VME distribution (as well as associated estimates of uncertainty), and their relative vulnerability, for a large region of the New Zealand marine environment.Identifying areas most likely to represent VME (rather than simply VME indicator taxa) provides much needed quantitative estimates of the most vulnerable habitats and facilitates an evidence-based approach to management (Gros et al., 2022).2Methods

(
Brodie and Clark, 2003;Helson et al., 2010).Some of the areas with high VME indices identi ed by this study, such as the Graveyard and Andes seamount complexes, have been well studied and VMEs in the form of coral reefs are known to occur on these seamounts (Bowden et al., 2019) supporting the contention that high value indices likely indicate the presence of VMEs.But other areas, such as the north-east Solander Trough and the north-west Bounty Trough have not been well studied, and thus present opportunities to independently assess the usefulness of the spatial predictions and VME indices method.Furthermore, Gros et al. (

Figure 1 Map
Figure 1

Table 2 .
Explanatory variables considered for use in Joint Species Distribution Models (JSDMs).The environmental variables are grouped into four categories: sea oor characteristics, water chemistry, water physics, and productivity.The symbol * represents those variables upscaled by Georgian et al. (2019) and the symbol ++ those variables upscaled by Stephenson et al. (2022).Highlighted cells represent those variables selected in the nal model (see section Joint Species Distribution Modelling for further information -spatial distributions of these explanatory variables are available in Figure S1, Supplementary Materials 1).
).A vulnerability to bottom shing score for each VME indicator taxon was assigned based on scores calculated for morphotaxa by Gros et al. (2023), which included consideration of six criteria: habitat forming; rare or unique; fragility; life history; larval dispersal and sessility (following Morato et al. (2018); Burgos et al. (2020)), which align with the FAO criteria for characterising VMEs (FAO 2009).Where our VME indicator taxa contained multiple morphotaxa, we used the mean vulnerability scores for each morphotaxa (morphotaxa and vulnerability scores are available in Table

Table 3 .
Overview of taxa records and JSDM model fits (ordered alphabetically by taxa, with VME indicator taxa used in this study shaded orange).Number of occurrences, prevalence (occurrences / number of samples), mean abundance (excluding absences), ten-fold cross-validated AUC scores of presence-absence model and R 2 of abundance conditional on presence model.Model fit metrics highlighted in red indicate models with poor model performance (subjectively defined as < 0.65 AUC and < 0.05 R 2 ).
(Roberts and Hirsh eld, 2004;Clark et al., 2016;Goode et al., 2020)conditional joint predictions shows promise when predicting distributions under future conditions due to climate change because it may allow environmental preferences of species and their interactions to be better captured than simply predicting each species separately under future conditions(Zhang et al., 2020b;Stephenson et al., 2022a).This capability will be particularly relevant for New Zealand because the region is predicted to be a hotspot of climate change-related impacts(Rickard et al., 2016;Law et al., 2018).Conditional joint predictions may also facilitate the identi cation of community assemblages.Further work to develop working de nitions of VMEs (consider the community assemblages) for the New Zealand regions would be of interest.However, care must be taken when interpreting the inter-taxa co-occurrence matrices from JSDM as they may not always represent biotic interactions but can instead re ect other spatial processes not accounted for in the model (i.e., missing important co-variates)(Ovaskainen et al., 2016a;Poggiato et al., 2021).4.2 Fishing impacts and identi cation of VMEsMany benthic taxa, particularly deep-sea taxa, are vulnerable to the direct and indirect impacts of shing, in particular, bottom-contact trawling (e.g.,Reed et al., 2007).Several taxa in our study, identi ed as VME indicator taxa, were found to have negative relationships between taxa occurrence and abundance with bottom trawl shing distribution and effort in the JSDM models.Similarly to other studies, we observed that these relationships varied based on biological traits (i.e., vulnerability), and the spatial overlap between taxa and the distribution of shing(Roberts and Hirsh eld, 2004;Clark et al., 2016;Goode et al., 2020).