Predicting the current and future suitable-habitat distribution of tropical adult and juvenile targeted fishes in multi-sector fisheries of central Queensland, Australia

Keywords: catch per unit effort, climate change, coastal, conservation, estuarine, Great Barrier Reef, habitat-suitability model, tropical.

Coastal and estuarine ecosystems provide human societies with many ecosystem services, particularly coastal protection, as a nursery for important fish species and carbon sequestration. Although coastal wetlands are highly valued (second only to coral reefs in value per 10 000 m²) (Costanza et al. 2014; Barbier 2019), these habitats are some of the most threatened natural systems globally and managing them is complicated because of the ecological connectivity among multiple realms (marine, freshwater and terrestrial ecosystems; Edwards et al. 2010; Sheaves et al. 2014; Reuter et al. 2016). Many marine species distributions are dynamic and strongly linked to temperature preferences (Cheung et al. 2012; Sunday et al. 2012). Climate warming is changing the distribution of coastal and estuarine species, which alters the structure and functioning of these communities (Galaiduk et al. 2018). In marine systems, species redistributions are occurring approximately one order of magnitude faster than in terrestrial systems (Chen et al. 2011; Poloczanska et al. 2013), because marine ectotherms live in habitats closer to their thermal limits and have a narrower thermal tolerance than do terrestrial species (Pinsky et al. 2019).

The potential impacts of climate change on temperate aquatic faunal assemblages are well studied (Jerry et al. 2013), but most of the available information is for northern temperate regions (Poloczanska et al. 2016; Clarke et al. 2020). It is inherently difficult to disentangle factors determining coastal and estuarine ecosystem functioning and isolating possible effects of habitat change. The complication arises from multiple interactions among complex oceanographic and terrestrial processes that determine how many species need to shift or locate accessible thermal refuges (Lauchlan and Nagelkerken 2020). There is currently a real void in knowledge on how tropical coastal assemblages will respond to climate change and what adaptation strategies may be available for decision-makers in regional communities to minimise impacts (Jerry et al. 2013; Poloczanska et al. 2016); particularly for commercially valuable tropical-specific fisheries such as barramundi, snapper and mackerels. Determining suitable habitat and spatial distribution is necessary to improve the conservation and management of exploited tropical fish species (Champion et al. 2021). Predictive habitat modelling has become an increasingly useful technique for decision-makers to estimate the patterns of species distribution and to include these in conservation strategies (Syah et al. 2016).

The niche concept relative to species distribution modelling is considered in the Hutchinsonian manner as a hypervolume in multivariate space depicting a species environmental limitations congruent with a stable population (Elith and Franklin 2013). Environmental constraints delimit the species fundamental niche, determining the space and resources that an organism can exploit to maintain populations given its intrinsic physiological limits (Yates et al. 2018). Processes and drivers of range dynamics and community structure act at different spatial scales and patterns observable at one scale may be driven by processes at other scales (Chave 2013). Physiological limitations (and therefore variability in environmental conditions) will help define the occurrence of species at a given location and, consequently, the structure of interaction networks (Franklin 2010; Thuiller et al. 2013). Coastal and estuarine ecosystems are characterised by their high spatial and temporal variability, estuarine species occupy more of their niche than their marine counterparts; thus, they are more vulnerable to climate-change stressors. Empirical studies with fish in estuaries have shown that many species live close to their upper limits of thermal tolerance and may face thermal stress under climate warming (Madeira et al. 2012; Lauchlan and Nagelkerken 2020).

Estuary-dependent fisheries species are important because they contribute 75% of the total value of Australia’s commercial fisheries catches and 90% by numbers of Australia’s recreational fish catch (Creighton 2013; Creighton et al. 2015). Understanding the spatio-temporal variability of key environmental variables within commercial fishing grounds is crucial in resource management to elucidate species distributions and their habitat associations (Naimullah et al. 2020). Moreover, many finfishes utilise a mosaic of habitats (Nagelkerken et al. 2015) and only few, e.g. barramundi, are confined to a single juvenile habitat (Sheaves et al. 2015). Furthermore estuary-dependent fisheries species are at particular risk because they are affected by both changes in natural flow regimes driven by climate change and those generated as a result of coastal land-use and land-cover change (clearing, modifying coastal habitats and artificial barriers to flow; Sheaves et al. 2014). For example, one of the highest risks to the Great Barrier Reef that has been identified by the Australian Government is the degradation to coastal habitats and connectivity impairment as a result of land-use changes disrupting the ecosystems of the region (Commonwealth of Australia 2018). Understanding coastal and estuarine ecosystems is critical to meeting the UN sustainable development goal (SDG) 14, which aims to ‘conserve and sustainably use the oceans, seas and marine resources for sustainable development’ (United Nations General Assembly 2017, p. 18). Predictive habitat modelling supports SDG 14, to help reduce negative anthropogenic impacts through improved management and control of human activities.

Habitat suitability models or species distribution models (SDMs) have the advantage of being tractable, easy to interpret and permit the predictability of phenomena that depend on differences between components. Furthermore, SDMs can explore the challenges of planning at the appropriate scale and can implicitly capture many complex ecological responses. Our modelling approach involves two conceptually different modelling techniques that cope with non-linear empirical relationships between responses and environmental predictors (Bučas et al. 2013). We used two SDMs in this study: the maximum entropy model (MaxEnt) and the generalised additive model (GAM; Sahri et al. 2020). We built the models using environmental data relevant to coastal systems, e.g. habitat mapping with remote-sensing data, salinity, sea-surface temperature. We incorporate climate change at current conditions, near future (2015–2054) and distant future (2055–2100). The objectives of the study were (1) to quantify the effect of recent past and current (2004–2019) environmental conditions on estuary-dependent fisheries in a large tropical coastal region, (2) to examine how those distributions are likely to change as a consequence of climate change, and (3) to inform regional coastal planning, conservation efforts and policymakers. The overarching objective was to find options that maximise the chance of retaining a productive near-shore fishery while accommodating climate change and human responses to climate change. This work will inform the prioritisation of conservation actions for the long-term persistence of coastal and estuarine species.

Our study area is located within the north-eastern coast drainage division of central Queensland, encompassing the Mackay Whitsunday Natural Resource Region and part of the Isaac Region, central Queensland, Australia (Supplementary Fig. S1). Cape Palmerston National Park is positioned in the Ince Bay Receiving Waters adjacent to the World Heritage listed Great Barrier Reef. The Shoalwater and Corio Bays Area Ramsar site is located south of Cape Palmerston National park. A large proportion of the marine waters in the Ramsar site are included in marine parks (Commonwealth and Queensland), including the Great Barrier Reef Marine Park (Commonwealth) and Great Barrier Reef Coast Marine Park (State; Department of Agriculture, Water and the Environment 2019). The Ramsar wetland supports a broad range of natural values, including nationally or internationally threatened wetland species, significant species diversity and large populations of waterbirds, green turtles, dugong and fish that use the site for vital life-history functions such as roosting, nesting, feeding and breeding. Cape Palmerston National Park is listed as a Category II protected area on the International Union for Conservation of Nature (IUCN) World Database on Protected Areas (International Union for Conservation of Nature 2020) and covers 72.0 km². The Shoalwater and Corio Bays Area Ramsar site covers 2019.7 km².

The primary intensive land use in the region is the cultivation of sugar cane, making up 18% of the catchment area, with Mackay being the largest sugar-producing region in Australia (Folkers et al. 2014). Grazing is also an important land use, accounting for 42% (Reef Catchments 2014). The estuaries in the region directly support several commercial fisheries, e.g. East Coast Inshore Finfish Fishery, East Coast Otter Trawl Fishery, and Line Fishery (Reef) (Pascoe et al. 2016). Additionally, recreational fishing is a considerable activity in the region, with 25% of the population participating in fishing for recreation, which is a far greater percentage than the state average of 15% (Webley et al. 2015). Mangroves and associated communities cover 879.94 km² of tidal land in the region, with nine wetland areas recognised as nationally important (Reef Catchments 2013) and one Ramsar wetland recognised as internationally important. The study area is located between latitude 20°22′–22°58′S and longitude 147°19′–150°46′E (Fig. S1).

We obtained commercial fisheries logbook grid records at a resolution 0.50° from the Queensland Department of Agriculture and Fisheries for the years 2004–2019. Georeferenced presence points with catch site, catch weight (kg), operational days and number of licences was included in the data. The data contain 59 estuarine fish species. Species are categorised as estuarine if their ontogenetic (juveniles in our study) or adult stages live within estuaries of the central Queensland coast. We calculated monthly catch per unit effort (CPUE) for each species. For the MaxEnt models, occurrence data were also sourced from the Atlas of Living Australia (see https://fish.ala.org.au/, accessed 1 March 2021) and recreational fishing databases.

Landsat fractional cover data (four bands) were downloaded from the Australian Geoscience Data Cube (Geoscience Australia 2020) through an interface with the National Computational Infrastructure (NCI) at a 0.00025° resolution. The multi-temporal Landsat scenes are available as analysis ready, having been geometrically and radiometrically corrected. Quality assessments are applied to the data before conversion to surface reflectance (Lewis et al. 2017). The four bands comprised Band 1-bare (bare ground and rock), Band 2-green (photosynthetic vegetation), Band 3-non-green vegetation (indicative of drier habitats with less vegetative cover) and Band 4-model fitting error (Purss et al. 2015; Gill et al. 2017). We used images captured in April and May for the years 2004, 2006, 2013, 2015, 2017 and 2019. In tropical areas, images captured in autumn–winter are often cloud-free. As each image was a subset of the region of interest, we catalogued the images for each year in a mosaic dataset, specifying footprints and boundary lines masked to our entire study area. Our purpose for including vegetation fractional cover is to delineate the coastline with a time-series of satellite scenes characterising estuaries, particularly mangrove forests. Because mangroves are vegetated areas that are persistently green throughout the year, they can be distinguished from grasses, wetland vegetation (e.g. tussock grassland, forbland) and woodlands (e.g. estuarine wetlands of Eucalpytus spp. or Acacia spp.), which follow a more fluctuating and seasonal trend in greenness (Lymburner et al. 2020).

To constrain the analysis to the coastal region, we obtained mangrove and seagrass distribution data, which we masked to our study site. Mangrove forests are a critically important habitat providing nursery and breeding areas for many commercial and non-commercial finfish species, such as barramundi (Lates calcarifer Bloch) and dusky flathead (Platycephalus fuscus Cuvier; Sheaves et al. 2016; Great Barrier Reef Marine Park Authority 2019). The data consisted of the Australian subset of the global temporal maps of mangrove extent generated by the Global Mangrove Watch (GMW; Bunting et al. 2018; Lymburner et al. 2020). Although we acknowledge that there is no causal relationship between mangrove area and fisheries catches (Sheaves et al. 2020), it is well recognised that mangroves are a component of seascape nurseries and a breeding ground for several species of finfish, thereby enhancing the productivity of fisheries (Saenger et al. 2013; Nagelkerken et al. 2015; Sheaves et al. 2015). Similarly, many finfish species are dependent on seagrass meadows for food, shelter or some part of their life-cycle (Coles et al. 2015). Seagrass spatial data for the Great Barrier Reef World Heritage Area from 1984 to 2014 was sourced from the eAtlas Data Catalogue (Carter et al. 2016, 2017).

The purpose of the models is to determine which variables explain more of the model variance. Variables such as sea-surface temperature (SST) anomaly, mean temperature, mean sea-surface salinity (SSS) and bathymetry were selected (Supplementary Table S1) on the basis of their identified importance in influencing tropical estuarine finfish distribution from published literature (Semmens et al. 2010; Jerry et al. 2013; Olds et al. 2013; Sheaves et al. 2016). SST anomaly data for the same years as the satellite images were accessed from the Australian Bureau of Meteorology (BOM, see http://www.bom.gov.au/marine/sst.shtml). Hydrology data (salinity and dissolved oxygen) collected from the Commonwealth Scientific Industrial Research Organisation (CSIRO) Marine National Facility RV Investigator voyage IN2016 were acquired for the region. Hydrology samples were collected from Niskin bottles sampled at various depths during conductivity–temperature–depth deployments (Sherrin and Schwanger 2016). Dissolved oxygen (DO) is a key variable for species occurrences and behaviour and can oscillate with tidal fluctuations in estuarine and coastal areas (Sheaves 2017; Feitosa et al. 2020).

Mean monthly values of temperature and practical salinity of the water column were acquired from the Australian Shelf Salinity Data Atlas, a collection of data assembled from different organisations by the Integrated Marine Observing System (IMOS; Integrated Marine Observing System, CSIRO Oceans and Atmosphere, Australian Institute of Marine Science, Royal Australian Navy Hydrography and METOC Branch, Defence Science and Technology Organisation, Department of Defence, Australian Government 2015). We selected a time-series of 2004–2014. We used a geostatistical kriging interpolation method (Fletcher and Fortin 2018) to generate a density surface for the environmental covariates derived from the CSIRO and IMOS hydrology data, including dissolved oxygen and salinity for the MaxEnt analysis. We used a bathymetry model, which is the compilation of all available source bathymetry data within the Great Barrier Reef into a 100-m-resolution Digital Elevation Model (Bearman 2010).

The environmental predictors were assessed for co-correlation using a Pearson r correlation. Of the 25 variables, eight variables were co-correlated below r² of 0.9, a threshold indicative of excessive autocorrelation between pairs of predictors (Phillips et al. 2006; Feng et al. 2019; Jones et al. 2019) and were subsequently used for species distribution modelling (Table S1). Finally, all variables were clipped to the study area and rasterised at 0.001° resolution (~111 m).

First, to reduce the effect of sampling bias and the resultant overfit in the model, it is appropriate to correct for spatial autocorrelation (spatial clusters) of occurrence points (Boria et al. 2014). We used the spatially rarefy occurrence data tool in the SDM toolbox to truncate occurrence localities to a single point within the user-specified Euclidian distance, maximising the number of spatially independent localities and reducing spatial autocorrelation (Brown and Anderson 2014). Subsequently, 21 estuarine fisheries species were used in the MaxEnt analysis (Supplementary Table S2) and 10 estuarine finfish species in the GAM analysis (Supplementary Table S3). These species are most important for both recreational and commercial fisheries in terms of catch rates on the central Queensland coast (Department of Agriculture and Fisheries 2022). We chose the following five finfish species on which to focus the study in the MaxEnt analysis: barramundi, listed as Least Concern on the IUCN Red List (Pal and Morgan 2019); red-throat emperor (Lethrinus miniatus J.R. Forster) listed as Least Concern on the IUCN Red List (Carpenter et al. 2016); dusky flathead listed as Not Evaluated on the IUCN Red List; blue threadfin salmon (Eleutheronema tetradactylum Shaw) listed as Not Evaluated on the IUCN Red List, but Endangered in the Persian Gulf (Motomura et al. 2015); and the flowery rockcod (Epinephelus fuscoguttatus Forsskål) listed as Vulnerable on the IUCN Red List (Rhodes et al. 2018). The five species have very different life histories, even though they all use estuaries. Barramundi can remain in estuaries throughout its life cycle, dusky flathead lives in soft-bottom estuaries and seagrass meadows, blue threadfin salmon frequents estuaries as nursery habitats, red-throat emperor and flowery rockcod reside in estuaries and seagrass meadows as juveniles, later moving to marine coastal areas (Great Barrier Reef Marine Park Authority 2011; Jerry et al. 2013; Olds et al. 2013; Atlas of Living Australia, see https://fish.ala.org.au/). These five species were chosen because of their high value in multi-species, multi-sector fisheries and possible susceptibility to climate warming in both the juvenile ecosystems and the recruitment response to adult populations.

We constructed spatial and temporal models of estuarine finfish species distributions on the basis of the covariates described above, by using the following two frameworks: (1) MaxEnt and (2) generalised additive models (GAM). We chose these methods because they have been demonstrated to provide strong predictive performance but differ considerably in their modelling approach (Fiedler et al. 2018; Sahri et al. 2020). Climate scenarios were run only on MaxEnt models.

MaxEnt is a deterministic algorithm that converges to the optimal (maximum entropy) probability distribution of a species by estimating a probability surface subject to constraints imposed by the occurrence data (Phillips et al. 2006; Owens et al. 2013). Recent studies have demonstrated the high performance of MaxEnt in modelling the distribution of presence-only species data and it has been increasingly used and applied to coastal and marine species (Robinson et al. 2017; Shahparian et al. 2017; Wang et al. 2018; Syah et al. 2020). MaxEnt has previously been used to simulate future potential geographical distributions of species under future climate-change scenarios (Herrera Montiel et al. 2019; Sudo et al. 2020).

The models generated were constructed using a default value of 1 for the regulation multiplier to reduce over-fitting during model training (Yusop and Mustapha 2019). The run type was a cross-validation method that divides the original samples into a set of training and testing of the models, namely, 75% of the location data used for training and the remaining 25% were used to test the predictive ability of the model. Our models were developed with default values for the algorithm settings, including feature selection, convergence threshold (10⁻⁵) and 10 000 background samples from the geographical space, with the maximum number of iterations set to 1000, allowing enough time for model convergence (Basher et al. 2014; Deb et al. 2017). The program was set to run jackknife tests to measure variable importance, create response curves and write background predictions to be used for model assessment and validation (Gullage et al. 2017; Charrua et al. 2020). A calibrated model portrays species’ response to environmental gradients, so response curves offer a proficient method to visualise models (Owens et al. 2013; Qiao et al. 2019). Monthly records were input to the models from 2004 to 2019, by using the same catch survey data (presence or true absence data). A receiver operating characteristic (ROC) analysis where area under the curve (AUC) represents the capability of the model to adequately predict presence (sensitivity) and absence (specificity) was provided by MaxEnt to evaluate the model, with the metric varying from 0 to 1 (Fernández-Manso and Quintano 2020). AUC values of >0.7 indicate models with acceptable predictive performance (Feitosa et al. 2020). Three independent models were run with three time scenarios. Except for SST, environmental variables were kept the same as we investigated the predicted effect of SST on suitable finfish and elasmobranch habitat. A correlation matrix was generated as a measure of dependency between raster layers.

GAMs are a set of linear models that allow the linear functions of the predictors to be replaced by arbitrary non-parametric smooth functions to relate explanatory variables to the dependent variable. GAMs offer a flexible modelling framework that can handle a variety of complex data (Wood 2017). Moreover, GAMs can explore non-linear relationships between a dependent variable such as the local relative abundance index, catch per unit effort (CPUE) and multiple predictor variables (Tseng et al. 2013). GAM model construction involves a choice of smooth functions (Franklin and Miller 2010). We chose a Gaussian distribution with canonical (identity) link function, cubic smoothing spline and low-rank tensor product smooth to add modelling functionality (Wood 2006). GAMs model the monthly and interannual response of the species, which in our study is CPUE (Villarino et al. 2015; Erauskin-Extramiana et al. 2020), the number of positions at site level (0.001°) in each 0.5° (55.5 km) fishing grid square corresponding to the total catch (kg) by number of days fished (Silva et al. 2015) within the central Queensland coastal domain. Spatio-temporal factors were considered by fitting a three-dimensional smoother to the product of three variables, namely, latitude, longitude and month. Finally, these factors were taken as variables in the analysis using GAMs fit with given-presence data from which absences (zero values) were excluded (Mbaye et al. 2020).

The environmental data included temporal (month, year), spatial (longitude, latitude), thermal conditions (SST) and practical salinity of the water column (Psal mean; Supplementary Table S4; Silva et al. 2015). SST and salinity (SSS) were explanatory variables in the analysis because of their recognised influence on the spatial distribution of tropical estuarine finfish. We used temperature and salinity from the IMOS dataset because it provided a monthly time series, and we converged the data using georeferenced occurrence points. The models were trained by establishing a relationship between recent past and current oceanographic conditions (2004–2019) and the explanatory variable at the known occurrence points (Galaiduk et al. 2017). Two distinct models were built to avoid confounding issues with multiple collinearity (Zuur et al. 2010; Mourato et al. 2014), namely (1) spatio-temporal, including longitude, latitude, month, weight (kg, the total catch for each species), mean temperature and mean salinity, and (2) oceanographic, including month, year, weight (kg), mean temperature and mean salinity. A cubic smoothing spline method (França et al. 2012) and tensor products were chosen as the smoothers for the separate covariates and interaction terms respectively (Yi et al. 2016). Terms denoted by ‘ti()’ are non-parametric smoothed terms. Candidate GAMs were constructed to model estuarine finfish monthly and yearly distributions and for exploring factors potentially influencing finfish species in the coastal region of the central Great Barrier Reef. Preferred models were selected according to the highest R²_(adj) value, the lowest Akaike information criterion (AIC) value, generalised cross-validation (GCV) value and the diagnostic plots (residuals distribution plot, quantile-quantile plot and response v. fitted values plot) (Elepathage et al. 2019). We focused the analysis on the following three species of commercially valuable finfish on the central Queensland coast: barramundi, red-throat emperor and blue threadfin salmon. A correlation matrix was generated as a measure of dependency among variables.

We used outputs from two new state-of-the-art climate models under the Coupled Model Inter-comparison Project Phase 6 (CMIP6) to project SST to the future. The CMIP6 models include a coupled atmosphere–ocean general circulation model (AOGCM) with incremental improvements in the simulation of the climate in the Australian region (Eyring et al. 2016). Projections of Australian temperature and rainfall from the available CMIP6 ensemble broadly agree with those from CMIP5, except for a group of CMIP6 models with higher climate sensitivity and greater warming and increase in some extremes after 2050. We avoided these ‘hot models’ and chose two models from the medium range of climate sensitivity. CMIP6 projections are conducted under a new framework of socioeconomic as well as emission pathways, with some GCMs including climate processes and earth system elements that were not included in previous generations (Grose et al. 2020). We chose the climate experiment at the highest level of emissions because this provides the strongest climate-change signal, the SSP5-85, described as fossil-fuelled development, taking the highway that broadly resembles the previous CMIP5 RCP8.5 (‘business-as-usual’) simulations. Increased SST scenarios were downscaled to the same resolution as the BOM SST surfaces and established for the following time ranges: (1) Model GISS-E2-1-G, multidecadal time periods of 2015–2050 and 2051–2100; (2) Model GFDL-CM4, four 19-year time periods and one 5-year period of 2015–2034, 2035–2054, 2055–2074, 2075–2094, 2095–2100 (Supplementary Table S5). We describe the time periods of 2015–2050 and 2015–2034, 2035–2054 as being the early period of the 21st century, and the time periods of 2051–2100 and 2055–2074, 2075–2094, 2095–2100 as being the latter period of the 21st century. Geographical data were manipulated in the R statistical environment (ver. 4.1.3, R Foundation for Statistical Computing, Vienna, Austria, see https://www.R-project.org/) with the packages rgdal (ver. 1.5-25, R. Bivand, see https://cran.r-project.org/package=rgdal) and raster (ver. 3.4-13, R. J. Hijmans, see https://cran.r-project.org/package=raster). GAMs were built using the package mgcv (ver. 1.8-36, S. Wood, see https://cran.r-project.org/package=mgcv/; Wood 2011). MaxEnt model building and model prediction was undertaken with MaxEnt software (ver. 3.4.1, S. J. Phillips, M. Dudik and R. E. Shapire, see https://biodiversityinformatics.amnh.org/open_source/maxent/). Removal of highly correlated environment variables and graduated spatial rarefying of occurrence data were undertaken with a Python-based GIS toolkit, SDMtoolbox Pro (ver. 2.5, see http://www.sdmtoolbox.org/; Brown and Anderson 2014). GIS analyses were performed in ArcGIS (ver. 10.8, Esri) and ArcGIS Pro (ver. 2.6.1, Esri).

The aim of this research was two-fold, namely, to quantify the effect of recent past and current (2004–2019) environmental conditions on estuarine-dependent fisheries in a large tropical coastal region and to how those distributions are likely to change as a consequence of climate change.

As expected, the correlation matrix shows a positive correlation between variables salinity and rainfall and also the variables mangrove presence and the Landsat vegetation band (Supplementary Table S6). We predict that finfish species such as red-throat emperor, dusky flathead (Fig. 1, Supplementary Fig. S2–S4a, b) and blue threadfin salmon (Supplementary Fig. S5–S7a, b) will have substantially reduced suitable habitat across the region in the early period (2015–2054) and latter period (2055–2100) of the 21st century. Similarly, for the elasmobranch spot-tail shark (Carcharhinus sorrah). The impact of climate change on these species is significant. Surprisingly, some species are predicted to expand their distribution with increased temperature such as flowery rockcod (Fig. 2, Supplementary Fig. S8b) and barred javelin (Pomadasys kaakan). The expansion in distribution may be related to an increase in prey species with increased warming. Our models performed better than random (i.e. AUC > 0.5) and with significantly high values (AUC > 0.9) for all models (Table 1). We found consistency in time ranges across simulations with two global climate models in the latter half of the 21st century (2051–2100 and 2055–2100). Overall, bathymetry was the key explanatory variable explaining the current and future distribution of estuarine finfish of the central Queensland coast (Supplementary Fig. S7a, S9a), ranking highest (12 species, 57%) in variable importance in most simulations (Table S2). Jackknife tests also identified bathymetry as the variable with the highest model training, testing and AUC gain when used in isolation (i.e. most useful information by itself). SST for both global climate models in the early (5 species, 24%), latter (6 species 29%) and entire (5 species, 24%) periods of the 21st century had a strong contribution, defining the future niche space for 76% of the taxa. Percentage variable importance in the MaxEnt simulation of four commercially valuable central Queensland estuarine finfish species with the global climate models GISS-E2-1-G and GFDL-CM4 are shown in Table 2. Of the five species that had global climate models predict the entire period of the 21st century, dusky flathead shows prominent range shifts and movement away from the shallow coastal area of Shoalwater–Corio Bay Ramsar wetland to deeper water both north and south of the Ramsar wetland with an increase in SST (Fig. 1). North of the Ramsar wetland, a preferred location for the dusky flathead under climate warming, is Repulse Bay, a sheltered, deep-water basin (Fig 1, Fig. S1). As expected, salinity ranked highly for the euryhaline barramundi (Supplementary Fig. S9b–S10); however, the stenohaline characteristic of this species was also evident in the response curves (Fig. S9c, d).

Fig. 1.  Difference map between the spatio-temporal distribution in the current (2004–2019) climate and that under future climate conditions for dusky flathead (Platycephalus fuscus). The green colour indicates a positive shift and the brown colour indicates a negative shift in the study region, namely Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, central Queensland coast.

Fig. 2.  Difference map between the spatio-temporal distribution in the current (2004–2019) climate and that under future climate conditions for flowery rockcod (Epinephelus fuscoguttatus). The green colour indicates a positive shift and the brown colour indicates a negative shift in the study region, namely, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, central Queensland coast.

Table 1.  Climate models and training and test AUC for the five focal species.

Table 2.  Percentage variable importance in the MaxEnt simulation of four commercially valuable central Queensland estuarine finfish species with the global climate models GISS-E2-1-G and GFDL-CM4.

Response curves describe the relationship between probability of occurrence and explanatory variables (Jones et al. 2019). The curves show the mean response of the 10 replicate MaxEnt runs (red) and the mean ± s.d. (blue). The response curves for our focal species show that logistic probability, and thus habitat suitability, is concentrated in shallow depths along the coastline and in the estuaries of our study region (Figs S7a, S9a). The strength of the response curve and the location of the point at which depth no longer led to increasing occurrence differed among species. High probabilities of occurrence are at depths of <20 m for the flowery rockcod, dusky flathead (strong coastal habitation; Fig. S4c) and barramundi, for which the relationship between bathymetry and presence suitability peaked at a depth of 0–5 m in a distinct bell-shaped curve (Fig. S9a). The response curve for the blue threadfin salmon showed a steep depth-related reduction from <50 m (Fig. S7a).

Bathymetry was followed by SST in the earlier climatic periods (2015–2050 and 2035–2054) as a variable of importance for the red-throat emperor and painted sweetlip (Diagramma pictum). SST in the latter climatic period (2051–2100) and the entire period (2015–2100) was important for the distribution of dusky flathead, largetooth flounder, rockcods (Epinephelus spp.), spot-tail shark, blackspot snapper (Lutjanus fulviflamma), blue threadfin salmon, giant trevally (Caranx ignobilis), bigeye trevally (Caranx sexfasciatus) and purple tuskfish (Choerodon cephalotes) (Table S2). Salinity was important for the barramundi (lower range), red-throat emperor (upper range) and the flowery rockcod (lower range) in both climate simulations. Dissolved oxygen ranked highly for species such as grey mackerel (Scomberomorus semifasciatus), stripey snapper (Lutjanus carponotatus), red-throat emperor, whose preference is for 205–210 µM of DO in seawater (Supplementary Fig. S11b), and blue threadfin salmon, whose preference is for 205–209 µM of DO in seawater (Fig. S7c). Fractional cover mosaics 2004, 2006, 2009, 2017 and 2019 are categorised as important for many species of estuarine finfish, including dusky flathead (Fig. S4d), flowery rockcod (Fig. S7d), blacktip rockcod (Epinephelus fasciatus), goldspotted rockcod (Epinephelus coioides), moses snapper (Lutjanus russellii), king threadfin salmon (Polydactylus macrochir) and purple tuskfish (Table S2). We suggest that the photosynthetic (Band 2-green) component of the fractional cover mosaics correlates with mangrove presence within estuaries and along marine shorelines. Seagrass distribution had the lowest ranking for all species and simulations (Table 2).

Photosynthetic vegetation (~25–62%), mangroves and salinity were important explanatory variables for predicted habitat suitability of barramundi in the latter climatic period (Fig. S9d–f). Considering that barramundi is a ubiquitous estuarine resident, it is unsurprising that the probability of occurrence increases with the density of mangroves and fractional cover (2006), both variables showing a positive relationship in the response curves (Fig. S9e, f).

As expected, the correlation matrix shows a positive correlation between the variables CPUE and weight, and a negative correlation between the variables temperature and salinity (Supplementary Table S7). Overall, the candidate spatio-temporal models explained 91–95% of deviance in the data for 10 estuarine, commercially valuable finfish species on the central Queensland coast (Supplementary Table S8). The preferred models (Models 1 and 2) explained 95.6 and 95.3% of the deviance in the data respectively, and had the lowers AIC values (Table S8). Generally, the graphical diagnostics of the models (QQ-plot, histogram, residuals v. linear predictors and response v. fitted values) were normal, and we found no violation of independence (Supplementary Fig. S12). Model 1 included spatial and temporal terms (longitude, latitude, month, year) and a joint tensor product smooth, including the CPUE elements, days fished and catch weight. The relative magnitudes for the monthly and yearly factors with a smooth on fishing effort (i.e. monthly fluctuation v. long-term trend; Supplementary Fig. S13). All smooth terms were highly significant (R²_(adj) = 0.956, P < 0.0001; Supplementary Tables S8–S9). Similarly, Model 2 included spatio-temporal factors but with a bivariate smooth term allowed to change monthly, as well as an environmental factor, mean temperature. CPUE was highest at longitude 149.75 and latitude −20.85 (Supplementary Fig. S14). Additionally, the smooths on the bivariate term (longitude by month), latitude and weight in Model 2 were highly significant (R²_(adj) = 0.952, P < 0.0001). There was a clear positive response of mean temperature to catch rate and a negative response of mean salinity to catch rate (Fig. S14). The smooth on mean temperature was significant at the 0.05 level (P < 0.05; Supplementary Table S10). Further, maximum finfish abundance (CPUE) corresponded with the month of December at a longitude of 149.75, according to the three-dimensional graphic output from GAM constructed with a thin plate regression spline (Fig. 3). Finally, the bivariate smooth term on month and longitude produced a contour plot showing peak abundance in December (Fig. 4). Model 3 was the third-most highly ranked spatio-temporal model, with an interaction of month and longitude, and a joint smooth on environmental factors, mean temperature and salinity (Supplementary Fig. S15). The smooths on the interaction term, weight and mean temperature were highly significant (R²_(adj) = 0.919, P < 0.0001). The smooth on salinity was significant at the 0.05 level (P = 0.05; Supplementary Table S11). The highest-ranked model showed the relationship between the CPUE of 10 estuarine, finfish species and explanatory variables derived from the GAM (Supplementary Fig. S12).

Fig. 3.  Three-dimensional graphic output from Model 2 GAM constructed with a thin-plate regression spline, a bivariate smooth term allowed to change monthly, as well as an environmental factor, mean temperature, spatio-temporal model for 10 species of commercially valuable finfish, central Queensland coast.

Fig. 4.  Contour plot from Model 2 GAM constructed with a thin-plate regression spline, a bivariate smooth term allowed to change monthly, as well as an environmental factor, mean temperature for 10 species of commercially valuable finfish, central Queensland coast.

Oceanographic models explained 90–94% of deviance in the data for 10 estuarine, commercially valuable finfish species on the central Queensland coast (Table S8). Mean temperature had greater relevance and a lower AIC value for the distribution of estuarine finfish than did salinity. The two most highly ranked models, Models 7 and 8, included a smooth on all temporal terms, as well as the environmental factors mean temperature and salinity respectively. All smooth terms were highly significant in Model 7 (R²_(adj) = 0.94, P < 0.0001; Supplementary Table S12). Temporal and fishing effort terms in Model 8 were highly significant (R²_(adj) = 0.94, P < 0.0001); however, salinity was less significant (P < 0.05; Supplementary Table S13).

The most highly ranked models for our three focal finfish species (barramundi, red-throat emperor and blue threadfin salmon) were constructed with a temporal component and a joint smooth on catch weight and mean temperature. We chose a cubic smoothing spline method and degrees of freedom (knots) between 3 and 12. All smooth terms in the oceanographic models for our focal species were highly significant (barramundi R²_(adj) = 0.954, P < 0.0001; red-throat emperor R²_(adj) = 0.942, P < 0.0001; blue threadfin salmon R²_(adj) = 0.72, P < 0.0001; Supplementary Tables S14–S16).

We found that tropical estuary-dependent finfish and elasmobranch abundance and distribution were influenced by changes in the environmental parameters from the variations of SST, SSS, DO and photosynthetic cover of the shoreline across a large coastal region (2641 km²) adjacent to the Great Barrier Reef (GBR) in Australia. In particular, the results showed that the spatial habitat patterns were explained predominantly by bathymetry and SST (Lee et al. 2018). Changes in environmental temperature may lead to non-linear responses in abundance among marine species; therefore, small increases in temperature can have large impacts on predicted outcomes (Flanagan et al. 2019). A high greenhouse-gas emissions scenario has been predicted to generate large-scale and long-term changes on regional biodiversity (Cheung et al. 2016) and, thus, a climate-driven reduction in fisheries productivity in tropical coastal regions (Lam et al. 2020).

Our findings highlighted that a future under CMIP6 climate scenario SSP5-85, which broadly agrees with the previous CMIP5 climate scenario RCP 8.5 (‘business-as-usual’ greenhouse-gas emissions scenario) will lead to substantial shifts in the distribution of species away from coastal habitats to deeper water. These shifts may lead to a host of fisheries governance challenges, including movements in populations across stock management boundaries (Pinsky and Fogarty 2012), outdated or contested management policies and quota allocations, greater costs and hence displaced fisherman, and changes in the price of catches (Jones et al. 2015). Many studies have found that global warming trends vary regionally and through time (Pinsky and Mantua 2014). In addition, substantial variability exists in climate impacts among species, for example in our study the impact of climate change on the finfish species red-throat emperor, dusky flathead and blue threadfin salmon, and the elasmobranch spot-tail shark is substantially reduced suitable habitat across the region in the early period (2015–2054) and latter period (2055–2100) of the 21st century. The species mostly affected by warming have limited acclimation capacity or live close to their upper thermal limits (Madeira et al. 2012; Pinsky et al. 2019).

A study on the sensitivity to climate change in 2030 across Australian fisheries found that because pelagic species can move large distances in search of suitable environments, they show flexible responses and low sensitivity (Fulton et al. 2018). This is not the case however, for coastal demersal finfish which show high sensitivity to changing environmental conditions (Pethybridge et al. 2020). An assemblage of coastal fishes restricted to the GBR has little scope for southern range shift because of limited dispersal, unsuitable habitat, competition with cold water species or asymmetric thermal responses in consumer-resource dynamics (Booth et al. 2011; Dell et al. 2014; Bonebrake et al. 2018). Consequently, coastal fish assemblages may dramatically contract in range, experience declining population abundance or local extinction (Bustamante et al. 2012). As identified by Fulton et al. (2018), species with the highest sensitivity included red-throat emperor, red emperor, king threadfin salmon and the spot-tail shark (listed as near threatened on the IUCN Red List). We found finfish species such as dusky flathead, largetooth flounder (Pseudorhombus arsius), blackspot snapper, blue threadfin salmon, giant trevally, bigeye trevally, purple tuskfish and rockcods shifted away from marginal coastal environments in response to a changing climate. By contrast, we predict some species to expand their distribution with increased temperature such as flowery rockcod and barred javelin.

Our findings corroborate previous studies of climate change impacts in tropical Australia, for example Welch et al. (2014) determined that variations in barramundi CPUE were significantly correlated to variations in both river height and rainfall. Diadromous fishes such as barramundi are the most susceptible to climate change induced habitat loss and degradation (Arthington et al. 2016). Consequently, altered precipitation patterns and decreased rainfall linked to salinity increase, is the environmental driver likely to have the greatest influence on the abundance of this species. Lower rainfall and higher salinities will have negative consequences for barramundi populations in our study region. Similarly, the flowery rockcod, listed as vulnerable on the IUCN Red List, may not show a response to warming, but will be constrained in the future by increasing sea-surface salinity attributed to surface evaporation and reduction in river flows. Although tropical waters are characterised by low oxygen (Deutsch et al. 2015) an increase in water temperature may affect water quality in tropical seas and further reduce dissolved oxygen affecting pelagic-neritic species with high metabolic demand such as grey mackerel, blue threadfin salmon and spot-tail shark. Estuary-dwelling juveniles of neritic species such as red-throat emperor, barred javelin and stripey snapper may have their developmental and survival rates altered as nursery habitats change (Creighton et al. 2015; Sheaves et al. 2015).

Multiple data types coupled with mixed method approaches of habitat suitability are gaining traction in conservation studies. These methods combined with increased sample size and geographic coverage are likely to provide the best and most informative evaluation of important habitats particularly for data-poor populations and geographic regions (Moore et al. 2016). The results of the MaxEnt analysis showed that most species (90%) will likely experience a significant range contraction in the future. If demersal finfish and elasmobranchs of the central Queensland coast conserve their historical niche and habitat preferences, that habitat may be reduced at a regional scale leading to local extinction of many species (Petatán-Ramírez et al. 2019).

In contrast to the MaxEnt analysis we did not include future climate in the GAM analysis. However, the GAM model was used to describe the spatio-temporal patterns of abundance, given presence of finfish and elasmobranch species (Mourato et al. 2014). The spatio-temporal models that explained the greatest variance in the data (Models 1 and 2) included the effect of the month, which incorporates possible seasonal variations in behavioural patterns that are not attributable to temperature or salinity. The effect of spatial location (longitude, latitude), included as a joint tensor smoothing function, proved to be highly significant. Environmental factors such as SST and SSS had a strong influence on predicted CPUE, and the distribution of demersal finfish and elasmobranch species of the central Queensland coast.

Here, we have presented some of the first fine-grained projections for distribution changes that encompass many of the geographic ranges for tropical estuary-dependent species on Australia’s north-eastern coast. Furthermore, our environmental variables included photosynthetic activity from remote-sensing datasets. Seasonal changes in the distribution of juvenile and adult finfish and elasmobranchs defined in the gridded data were largely described by the models, in particular, longitudinal shifts in distributions were well predicted. This seasonality in longitudinal distributions within the coastal regions of the GBR highlighted the importance of inclusion of a month term in the models. Model 1 explained 95.6% of the deviance in the data followed by Model 2, with a bivariate smooth term allowed to change monthly, which explained 93.5% of the variance in the data. The GAM model produced robust estimates of relative abundance (CPUE) and distribution along the central coast of the GBR and has been especially effective in capturing the spatial and temporal variability of CPUE. The spatio-temporal models explained marginally more (average 93.1%) of the variability in catch rates than did the oceanographic models (average 92.5%). The variation of nominal CPUE during the period 2004–2019 indicated that the highest monthly nominal CPUEs and the main fishing ground were observed in December (summer) at longitude 149.75 and latitude −20.85, within the GBR lagoon. Substrate composition in this region is mud and shell and the depth is ~41 m. From our analysis, the high-value target species included red-throat emperor, red emperor and tropical snappers. Little is known of the spawning behaviour of demersal finfish on the GBR, despite their commercial and recreational importance; thus, we cannot determine if the model prediction is due to spawning aggregations. However, some finfish, such as snappers (e.g. stripey snapper), rockcods (e.g. flowery rockcod, goldspotted rockcod), javelins (e.g. barred javelin), jewfish (e.g. black jewfish) and mackerels (e.g. grey mackerel), form spawning aggregations in large numbers (Russell and Pears 2007; Semmens et al. 2010). In addition, subtidal areas of consolidated and mobile sands and muds are important habitats for fishery species across the spectrum of managed coastal fisheries (Sheaves et al. 2016; Brown et al. 2019). Consequently, the greater CPUE at the site predicted by the model may be due to species’ reproductive or connectivity-based strategies. Empirical studies have shown red-throat emperor to be affected by temperature increases because their distribution is restricted to south of 18°S, in the central and southern GBR, a region expected to experience a greater climate-change influence than the north (Lough 2007; Munday et al. 2007; Northrop and Campbell 2020). The southern GBR particularly is affected by the climate-induced strengthening of the East Australian Current (Booth et al. 2011; Bustamante et al. 2012). Furthermore, the red-throat emperor has less site-specificity and greater mobility than do other demersal fish on the GBR, but an apparent upper thermal limit of ~28°C (Munday et al. 2008). Clearly, the catch rate of the red-throat emperor will reduce as their distribution shifts south. To make informed decisions, policy makers require knowledge of the distribution of exploited species to develop new fisheries or restructure the management of existing ones because of the non-stationary nature of climate drivers (Gomez et al. 2015).

For models to be used in the management of fishery species, model uncertainty should be assessed, and model limitations evaluated to understand cumulative errors (Gomez et al. 2015). Many studies have asserted that predictions from species distribution models are sensitive to model assumptions and uncertainties (Cheung et al. 2009); however, the predictive capacity of MaxEnt has been consistently ranked as one of the highest-performing among methods, despite operating on presence-only (with random background) data (Elith et al. 2006; Moore et al. 2016; Mugo and Saitoh 2020). Similarly, Hijmans and Graham (2006) reported that MaxEnt and GAM performed as well under current climates as under past and future climates.

An underlying assumption of species distribution models is that the species is at equilibrium with its environment and that relevant environmental gradients have been adequately sampled (Elith and Graham 2009). However, other processes are important to consider, such as population and fisheries dynamics, dispersal and trophic interactions, including an amplified warming signal and declining productivity response across tropical oceans (Chust et al. 2014; Erauskin-Extramiana et al. 2019). The effect of warming on ecological interactions will determine the predictions of biomass change for species (Brown et al. 2010). Mechanistic models can be useful for this task because they establish a causal relationship between species distribution and mechanistic variables, independently of the species records, but generally require detailed biological knowledge for their parameterisation (Jørgensen et al. 2012).

The integrated land–sea management initiative requires tools and approaches that provide robust projections of future changes, particularly within regional seas and their coastal habitats (Reuter et al. 2016; Peck et al. 2018). Within GBR coastal provinces, ecosystem effects and cumulative impacts on fishery resources are poorly understood (Great Barrier Reef Marine Park Authority 2014). Moreover, disparate jurisdictional responsibilities hinder assessment efforts. Australian fisheries management requires coordination over large spatial areas, across sectors and jurisdictions (Fulton et al. 2018). This study reflects a primary motive for producing these species projections, namely, to contribute a spatio-temporal approach that would be available to policy makers and managers considering climate adaptation of marine fisheries management (Morley et al. 2018; Torrejón-Magallanes et al. 2019).

In this paper, we have used habitat-suitability models to understand spatial and temporal patterns of species distributions so as to assess climate-change threats to biodiversity, with the objective to inform fisheries management on the tropical eastern coast of Australia. We focused on delineating key areas of the seascape that are crucial for different life-history stages of multiple species and may benefit exploited species, particularly those species that exhibit relatively localised movement patterns.

Given the uncertainty around climate-change policy within the government domain and the true magnitude of future climate change and resulting species redistribution, managers (and the mechanisms they use) will need to remain flexible and adaptable. This process involves accommodating the changing nature of the broader socioecological system, providing recognition of the trade-offs between the objectives held for the different sectors and incorporating these in the resource-management system. Conservation strategies that initiate system modifications in support of climate-change adaption through extending the protected-area estate (e.g. the protection of estuarine nursery habitats, or re-establishment of hydrological regimes to restore ecosystem function, enhance fecundity and migration efficacy) may be more appropriate than are traditional reference levels such as maximum sustainable yield in fisheries management. Climate change has an unequivocal impact on the juvenile and adult stages of valuable finfish and elasmobranch species; we have captured that effect by identifying important environmental drivers and incorporating the most up-to-date climate simulations under future scenarios.

Supplementary material is available online.

Some of the data generated in this study are available in the Supplementary material.

The authors declare that they have no conflicts of interest.

This research did not receive any specific funding.