Modelling current and future abundances of commercial benthic invertebrates using bathymetric LiDAR and oceanographic variables

A key requirement for managing commercial fisheries is understanding the geographic footprint of the resource, the level of exploitation and the potential impacts of changing climate or habitat conditions. The development of spatially explicit predictive models of species distributions combined with predictions of changing oceanographic conditions provides the opportunity to obtain new insights of species-habitat associations. Here, generalized linear models (GLMs) were used to model the abundance of two commercially important marine macro-invertebrates, blacklip abalone Haliotis rubra and long-spined sea urchin Centrostephanus rodgersii , along the coast of Victoria, Australia. We combined abundance data from fisheries independent diver surveys with environmental variables derived from bathymetric light detection and ranging (LiDAR) and oceanographic parameters derived from satellite imagery. The GLM was used to predict species responses to environmental gradients where reef complexity, sea surface temperature (SST) and depth were strongly associated with species distributions. The abundance of H. rubra declined with increasing summer SST. In comparison, the abundance of C. rodgersii increased with increasing winter SST. The GLM showed that the projected increase in ocean temperatures will likely lead to a decline in abundance across the H. rubra fishery. Conversely, a range expansion of C. rodgersii is likely due to the strengthening of the East Australian Current. For species that exhibit a high affinity to specific seascape features, this research demonstrated how recent advances in seabed mapping can allow the identification of areas with high conservation or fisheries value at a fine-scale relevant to resource exploitation across large geographic regions. This study demonstrates that both seafloor topography and SST explain the distribution of the commercial benthic species, H. rubra and C. rodgersii . Changes in the distributional range of these benthic species under future warming based on the fossil-intensive scenario were observed. Predicting the consequences of the interactive

of Victoria being about 28 t (~AU$123000) during 2009-10 (DSEWPC, 2011. Several studies in Australia (New South Wales and Tasmania) have shown that abalone are almost absent where urchin's dominate, although abalone often co-occur with urchins in fringe habitats (Wright et al., 2005;Strain & Johnson, 2009 It has been suggested that urchins and abalone compete for the same food source, which inhibit the settlement of abalone (Strain & Johnson, 2012). The range expansion of urchin barrens may also limit the resilience and productivity of the abalone industry in Australia and the effectiveness of management strategies aimed at arresting abalone stock declines due to disease, a changing climate and fishing pressure (Andrew & Underwood, 1992;Johnson et al., 2011).
The distribution of urchins, especially C. rodgersii, has historically been restricted to subtidal reefs along the coasts of New South Wales (NSW), eastern Victoria, and the Flinders Island Group in the Bass Strait (Andrew & O'Neill, 2000). In recent years, the range of C. rodgersii has expanded southwards to the east coast of Tasmania and westward along the coastline of Victoria. This phenomenon is likely due to strengthening of the East Australian Current, bringing warmer water containing urchin larvae further Ridgway, 2007 southwards ( ;Ling et al., 2009). The consequences of global warming are already being observed, with some local and regional environmental changes impacting species ranges (Parmesan et al., 1999;Thuiller, 2004;Lough & Hobday, 2011;Pecl er al., 2017). The Intergovernmental Panel on Climate Change (IPCC) assessment reports project that global temperature will rise between ~ 1.6°C and 6.4°IPCC, 2007 C by the end of the 21st century Betts et al., 2011). Oceanic waters in southeast Australia are experiencing rapid warming, three to four times the global average, making this region a hotspot for ocean temperature change (Ridgway, 2007). For instance, the mean sea surface temperature (SST) has already risen by 0.8°L ough, 2009 C on the east coast of Victoria and Tasmania since the 1960s ( ).
Shifts in species ranges depend on many factors, including life history traits, competition, and habitat requirements. Environmental factors, such as the physical characteristics of the seabed, influence the spatial patterns of species (Moore, Van Niel,

Author Manuscript
This article is protected by copyright. All rights reserved & Harvey, 2011; Rees, Jordan, Price, Coleman, & Davis, 2013). Seabed physical structure can be used to infer the complexity and connectivity of benthic habitats, especially with respect to the accuracy and detail that can be achieved as well as to the scale relevant to the resources exploited (Hedley et al., 2016;Valle et al., 2015). Recent advances in mapping technology, such as bathymetric light detection and ranging (LiDAR), allow the rapid collection of high resolution, full coverage bathymetry information over large geographic extents in nearshore waters. This technology has the potential to provide detailed information on the distribution and complexity of reefs that would, otherwise, be difficult to collect using other approaches (Pittman. Costa, & Wedding, 2013;Wedding, Friedlander, McGranaghan,Yost, & Monaco, 2008). For example, in comparison to the low resolution of bathymetric products derived from satellite altimetry (kilometres), the high resolution provided by LiDAR (metres) offers great potential for advancing our knowledge of the functional linkages between geophysical structures and ecological processes, especially in shallow coastal zones (Chust, Grande, Galparsoro, Uriarte, & Borja, 2010;Valle, Borja, Chust, Galparsoro, & Garmendia, 2011). Bathymetric LiDAR data has been effectively used to derive reef rugosity and complexity indices that correlate with fish abundance and species richness (Kuffner et al., 2007;Zawada & Brock, 2009;Zavalas, Ierodiaconou, Ryan, Rattray, & Monk, 2014). Such data is particularly useful for benthic species that might have certain requirements for recruitment success linked to substratum characteristics (Galparsoro, Borja, Bald, Liria, & Chust, 2009).
In this study, we used bathymetric LiDAR data and oceanic variable (i.e., SST) to model spatial variation in abalone and sea urchin abundance data collected by divers.
We used these models to define the importance of geophysical and climatic variables and predict species distributions over broad statewide scales. We then applied IPCC scenarios for SST to project future distributions of abalone and urchins, from which we quantified the potential impact of urchins on existing abalone fishing grounds. Our results are expected to inform the management of current abalone fishing grounds, as well as provide baseline information towards managing the stock into the future.

Author Manuscript
This article is protected by copyright. All rights reserved

Bathymetric and oceanographic variables
) was estimated by dividing the number of the individuals by the area of transects. Survey transects were randomly allocated between each pair of divers (three transects per diver per site), to overcome variability in estimates among observers.
Bathymetric LiDAR data were collected in 2007 using a LADS Mk II system coupled with a GEC-Marconi FIN3110 inertial motion sensing system and a dual frequency  (Table 1). These parameters might, in turn, influence the occurrence and distribution of abalone as well as the presence of the urchin and the propensity of overgrazing (Johnson et al., 2011;Ling & Johnson, 2012). In addition, the substrate was classified into three groups (reef; mixed; sediment) using a decision tree classification method integrating bathymetric
Mean Austral summer (December to February), winter (June to August) and annual SST (°F ranz, 2006 C) between 2002 and 2011 were derived from raw MODIS Aqua daily images using NASAs SeaDAS image processing software (version 6.1) at a spatial resolution of ~ 1 km. The SST algorithm used was the Ocean Biology Processing Group (OBPG) daytime long-wave SST algorithm ( ), which is based on MODIS band 31 (11 µm) and band 32 (12 µm). Daily images were aggregated to monthly images, and temporally averaged to obtain mean summer, winter and annual SST  Betts et al., 2011). Present-day SST values were extracted for the bathymetric LiDAR coverage at 1 km resolution (n = 3,000 pixels) using the nearest values. We then downscaled SST data to 5 m resolution to match the bathymetric LiDAR products using ordinary Kriging geo-statistical interpolation procedure in ArcGIS 10.2 (ESRI INC). Kriging has been shown to outperform other interpolation approaches (Murphy, Curriero, & Ball, 2010). The same procedure as for the present-day was applied to interpolate future summer and winter SST.
The interpolated present-day SST maps were validated using an independent dataset provided by Surf-Forecast (http://www.surf-forecast.com/) that reports daily

Author Manuscript
This article is protected by copyright. All rights reserved extracted for statistical analyses. Data from the survey sites located outside of the available bathymetric LiDAR coverage and those sites with inconsistent data that were added to the sampling design after 2004 were also excluded from the species distribution modelling explained in the next section. Therefore, abundance data from 122 survey sites with a consistent 10 year (2002-2011) series of survey records were retained in the analyses.

Analysis of species spatial distributions
Data exploration was conducted to identify outliers and determine colinearity between explanatory variables using the Pearson correlation coefficient and variance inflation factor (VIF) (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). Correlation coefficient thresholds of < 0.5 and VIF value of < 3 were used to retain only the least correlated variables. Strong correlations (> 0.7) were detected between slope, complexity and VRM. Considering higher variation explained by complexity and the important role of this variable as a potential determinant of benthic species habitats, it was retained to run the models whereas slope and VRM were excluded from the analysis. As high correlations (> 0.7) were also found between longitude, winter, summer and annual SST data, collinear variables were not included in the same model. Summer and winter SST data were respectively used in the abalone and urchin predictive models due to their high contributions to the model, as well as their relevance to the ecology of the species investigated.
A generalized linear modelling (GLM) approach was used to predict abundance of H. rubra and C. rodgersii as a function of bathymetric and oceanographic variables across the entire bathymetric LiDAR data coverage across the Victorian coastline. The GLM is a flexible generalization of ordinary linear regression for relating responses to linear combinations of predictor variables via a link function. This modelling method is suited for identifying ecological relationships and allows for non-linearity and nonconstant variance structures in the data (Guisan, Edwards, & Hastie, 2002). Due to the variations in H. rubra abundance, as well as the number of zero observations in urchin datasets, a GLM was developed using a negative binomial distribution with a log link function. Negative binomial GLMs are often used to address overdispersion in model residuals caused by high variation in abundance measurements (Hilbe, 2007;Zuur et al.,

Author Manuscript
This article is protected by copyright. All rights reserved 2009). Before fitting the final model, the data were randomly split into training and test datasets with 85% of the data for training and 15% for model validation.
Final fitted models were then used to predict the present and future abundance of H. rubra and C. rodgersii across the entire study area. To do this, raster datasets of the predictor variables were imported into the R program (R Core Team; version 2.15.3) using the "raster" package (Hijmans, 2014), and then a stack layer of the predictors was built to be used for the prediction of abundance for each of the pixels in the raster dataset. Predicted maps were then imported into ArcGIS for model visualization. Model selection was based on variable significance, percentage deviance explained and the Akaike information criterion (AIC). Goodness of fit of the GLMs was evaluated using a chi-square test based on the residual deviance and degrees of freedom with p > 0.05 indicating a good fit. Model diagnostics included examination of residual plots and the regression of fitted against observed values. Spatial independence was assessed by plotting the residuals of the best fitting model against their spatial coordinates. All GLMs were fitted using the "mgcv" (Wood, 2011) and "MASS" (Venables & Ripley, 2002) packages in R (R Core Team) version 2.15.3.

Potential impacts on abalone fishing grounds
The predicted H. rubra abundance maps were masked for each fishing zone in Victoria (Eastern, Central and Western fishing zones), and the percentage change in abundance between future and current predictions was estimated for each zone. To assess the potential impact of future urchin distributions on abalone reef habitats, the projected C.
rodgersii map for the year 2100 was used to evaluate how urchins impact productive abalone habitats, assuming reef geophysical characteristics will remain similar, and accounting for changes in SST.

RESULTS
H. rubra was detected along almost the entire coastline of Victoria where suitable areas of reef were present (Figure 5a). The highest and lowest H. rubra density (1.72 and 0.2 abalone m 2 , respectively) occurred in the eastern fishing zone. C. rodgersii was abundant at sites in eastern Victoria, near the New South Wales border, and was

Author Manuscript
This article is protected by copyright. All rights reserved generally absent from the western and central fishery zones during the survey period ( Figure 5b); however, isolated specimens (n = 3) were collected in the far west (20 km from Portland) of Victoria, outside of the survey period.
The negative binomial GLM predicting the abundance of H. rubra showed that substrate complexity and summer SST were the best predictors according to the model selection parameters (Table 2 and 3). Within this model, the abundance of H. rubra decreased with increasing summer SST ( Figure 6b, rodgersii was found in high abundance. The predicted present-day abundance of C. rodgersii covered known habitats occupied by the urchins, with the highest predicted abundance occurring in the far eastern part of Victoria. The GLM projection based on 2100 winter SST ranges (from the A1FI high emission scenario) showed that east coast of Victoria will likely support the expansion of areas with higher C. rodgersii abundance. There was also an increased probability of the distribution of C. rodgersii within other coastal fishing grounds in Victoria (Figures 7c, d).
°C and 6.4°C increase based on the IPCC A1FI high emission scenario) throughout the waters of Victoria by the end of the 21st century, with some variability across the state (Figure 8). The rate of warming was higher in the eastern zone; however, substantial increases in SST were also evident in central and western Victoria. Using the test data, it was found that the GLMs generally

Author Manuscript
This article is protected by copyright. All rights reserved predicted higher abundance with a positive Pearson correlation between observed and predicted values. There was no evidence of spatial dependence in the model residuals.
GLMs with the negative binomial distribution accounted for overdispersion (dispersion parameter close to 1), and none of the GLMs indicated patterns in residuals.
Differences in the abundance of H. rubra between current and future projections indicated 4% to 11% decline in abundance by the year 2100 across fishing zones, with the rate of loss being higher in the east ( Figure 9). However, the abundance of H. rubra might also increase on some reef habitats, such as specific sub-zones in the western fishing zone, even though as a whole losses were observed (Table 4). By quantifying the percentage of potentially impacted abalone reef areas, a future increase in the distribution of C. rodgersii by the year 2100 is expected to impact some important suitable reef habitats, especially within the abalone zone of eastern Victoria (Table 4).

DISCUSSION
Understanding the environmental parameters driving the distribution of species is key to understanding impacts of potential climate change scenarios, the implications of change in suitable habitat conditions and changes in species range. Our results confirm that both seafloor characteristics and climate-related factors determine current and future abundance distributions of commercial and ecologically valuable marine macroinvertebrates in southeast Australia.
Our predictive model showed that high seafloor complexity is the key abiotic predictor explaining the abundance patterns of both H. rubra and C. rodgersii as has been demonstrated at fine scales by SCUBA diving and towed underwater video (Johnson et al., 2011) and Autonomous Underwater Vehicles (Ling et al., 2016).
Complex benthic terrain supports a greater diversity and abundance of fish (Walker, Jordan, & Spieler, 2009;Moore et al., 2011) and invertebrate species (Holmes. Van Niel, Radford, Kendrick, & Grove, 2008;Galparsoro et al., 2009) than simple terrains, particularly in shallow areas. This is because complex terrains provide suitable cryptic habitat (i.e. shelter and protection from predators), as well as suitable feeding and

Author Manuscript
This article is protected by copyright. All rights reserved breeding grounds for benthic dwelling organisms (Kuffner et al., 2007;Alexander et al., 2014). In particular, gutters and surge channels with greater exposure help concentrate macro-algal drift enabling larger H. rubra to feed opportunistically whilst aggregated along vertical reef walls (Gorfine, 2002). Moreover, in a previous study by Ling and Johnson (2012) in assessing the relationship between habitat complexity and predation risk by lobsters, it was found that C. rodgersii had greater survival in crevice shelters compared to exposed urchins, indicating how reef complexity can impact predator-prey interactions.
Our study showed that the abundance of H. rubra was generally prevalent in shallow subtidal systems, as expected from previous studies, and across relatively similar depth ranges. Although there was no significant relationship between depth and the abundance of H. rubra, it is acknowledged that the depths observed were restricted to a diver depth range. In addition, the depth range of this study covers an important part of the ecological range of the species (2 m to 20 m), but the depth range was likely too narrow to find significant changes in abundance with depth. In contrast, the abundance of C. rodgersii was found to be greater on deeper substrate with high complexity with lower densities in shallow (< 5 m) inshore reefs (Johnson et al., 2005;Gorfine, Bell, Mills, & Lewis, 2012). In general, bathymetry and factors that co-vary with depth contribute towards structuring marine species by influencing light attenuation and exposure along with substrate complexity (Smith & Brown, 2002).
The results of this study supported direct observations of known habitats occupied by abalone and sea urchins. Our GLMs showed that SST was a significant predictor of the abundance for both species. C. rodgersii was highly correlated with warmer winter

Author Manuscript
This article is protected by copyright. All rights reserved The results support that both abalone and urchins are impacted by warming oceans. Temperature is the major factor affecting all biological processes, including those observed in marine ecosystems (Pörtner & Knust, 2007). Thus, it is important to understand how species thermo-tolerance affects breeding success, larval development and survival rate to infer how marine species will respond to climatic changes. Wild H.
rubra spawn between October and April, and exhibit restricted larval dispersal (Gorfine, 2002;Chick, 2010); however, recent studies have indicated gene flow across the abalone stock, with larval settlement probably being restricted by habitat rather than currents (Miller et al., 2016). The spawning of C. rodgersii around Australia occurs during winter, with fertilization peaking in August, and normal development especially being observed at temperatures over 12°C (Byrne & Andrew, 2013;Ling et al., 2009).
However, temperatures that exceed the tolerance of a given species suppress larval growth and survival, due to increased metabolic costs. For instance, experimental studies on the adaptive capacity of urchin (Heliocidaris erythrogramma) demonstrated deleterious effects when embryos and larvae were subjected to a +3°C-4°C increase above ambient temperatures (22°C-23°C) (Byrne, Selvakumaraswamy, Ho, Woolsey, E, & Nguyen, 2011); yet, those that normally developed at 26°C became juveniles, indicating the adaptive capacity of the species. Byrne et al. (2011) also suggested that H. erythrogramma may have an in-built plasticity to persist through ocean warming, and the species has also experienced gradual warming in local waters during recent decades, possibly resulting in evolutionary changes that have increased thermal tolerance. Furthermore, laboratory studies on the fertilization and development of C.

Author Manuscript
This article is protected by copyright. All rights reserved determining the resilience of echinoids to ocean warming, with existing studies suggesting that such resilience is high for C. rodgersii populations.
In contrast to C. rodgersii, H. rubra has less tolerance to higher temperatures, exhibiting different biological responses to temperature. H. rubra preferentially grows at ~ 17°C, with higher water temperature causing mortality to juveniles (Gilroy & Edwards, 1998). A recent study in South Australia suggested that abalone (H. laevigata and H. rubra) have a limited capacity to tolerate warmer waters during the summer months (Russell et al., 2012). This finding was supported by our results, which predicted that H. rubra would respond negatively to warmer summer temperatures, with a reduction in the fishery projected by 2100. In a recent study by Oliver et al. (2017) on the ecological impacts of marine heatwave in the Tasman Sea, unusual mortalities were observed among H. rubra populations in the east and southeast Tasmania. Oliver et al.
(2017) also suggested the potential to adapt to local thermal regimes, with mortalities occurring in cooler adapted abalone populations due to acute stress and elevated metabolic rates. Our study did not incorporate the thermo-tolerance limits and adaptation capacity of abalone and urchins to different temperature ranges in the species distribution model, due to data limitations and the inability of current models to account for physiological limits explicitly when predicting species geographic distributions (Sinclair, White, & Newell, 2010;Martínez, Arenas, Trilla, Viejo, & Carreño, 2015).
All models have some levels of uncertainty and subjectivity (strengths and weaknesses) in addressing how species relate to climate change (Beale & Lennon, 2012;Kerr & Dobrowski, 2013). Nevertheless, there tends to be a general lack of available habitat information at a scale relevant to exploitation. With the increasing availability of high resolution oceanographic and LiDAR/multibeam derived seafloor structure information, it is now possible to obtain information on how the ranges of species change at scales similar to those of terrestrial studies.
Our models indicated the potential for the distribution of C. rodgersii to expand its range under future climate change impacting on abalone reef habitats. However, in the assessment of the potential impact of urchins on abalone grounds, the effect of fishing was not included in the model due to the limitation in spatially and temporally consistent catch data for the survey sites. The role of fishing is important to consider in the future predictive models because if urchin abundance is low and/or kept low by

Author Manuscript
This article is protected by copyright. All rights reserved urchin harvesting then barrens and impacts of urchins on abalone will be minimal.
Nevertheless, some of the fishing sub-zones within the eastern fishing zone of Victoria have become heavily depleted of abalone, due to C. rodgersii range expansion.
Commercial abalone divers have observed that large areas of previously productive fishing grounds were impacted by formation of urchin barrens. In parallel, our findings indicate a high rate of occupancy by urchins in some of the fishing sub-zones. Culling urchins allows the recovery of canopy forming algae (Ling, Ibbott, & Sanderson, 2010;Strain & Johnson, 2013); thus, an urchin culling programme is currently being implemented in the eastern fishery zone of Victoria. Due to the costs of culling and the scale of urchin dispersal, decreasing the likelihood of urchin barren formation might be a more viable management option by protecting "beneficial predators", such as larger rock lobsters that feed on urchins but are not competitors of abalone (Pederson & Johnson, 2006;Babcock et al., 2010). The role of predators in regulating sea urchin abundance and distribution, and the impact of such prey-predator interactions on the resilience of marine ecosystems, is well known (Scheibling & Hamm, 1991;Hamilton & Caselle, 2015). However, implementing conservation strategies on the basis of predator activities requires comprehensive collaborations between managers and fisherman to prevent the overfishing of lobsters in abalone habitat that is also suitable for urchins. Recent fisheries assessments within the western fishing zone

Author Manuscript
This article is protected by copyright. All rights reserved

Author Manuscript
This article is protected by copyright. All rights reserved