Identifying refugia from the synergistic threats of climate change and invasive species

Climate change may reduce available habitat for native species, while simultaneously increasing suitable habitat for invasive species, which then compete with or predate on native species. Thus, climate change and invasive species can interact synergistically to negatively affect native species. It is important to identify climate refugia that are likely to be both suitable for native species and unsuitable for invasive species, under both present and future climate conditions. We propose a refugia habitat metric (RHM) based on ecological niche modelling. We demonstrate the utility of the metric via a case study of an endemic freshwater crayfish, or kōura (Paranephrops planifrons), which is threatened by both climate change and predation from the invasive brown bullhead catfish (Ameiurus nebulosus) in Aotearoa / New Zealand. We used maximum entropy (MaxEnt) ecological niche models to predict current and future habitat suitability for the two species across Aotearoa / New Zealand. By the period 2080–2100, suitable habitat will increase across the northern and western North Island for catfish, while suitable habitat for kōura will decrease overall and shift southwards and towards more mountainous regions. Using the refugia habitat metric, we identified areas of habitat within the current range of kōura, with significant potential refugia habitat outside the species current range. Using the refugia prioritization metric will allow conservation managers to identify habitat for protection and potentially translocation target sites for vulnerable native species.

1 Introduction

Climate change and invasive species can alter biodiversity and ecosystem processes. The impacts of climate change and invasive species are often considered independently, yet these two drivers of global environmental change interact in complex ways, such as enhancing the impacts of invasive species already present by skewing predatory and competitive interactions with natives, if the invader is better adapted to the changing conditions (Rahel et al., 2008). The presence of invasive competitors or predators may result in further declines in native species struggling to adapt to altered climate regimes. Thus, climate change can directly negatively affect a native species and indirectly affect it through facilitating invasive species (Rahel and Olden, 2008).

Figure 1 (a) Distribution of brown bullhead catfish (Ameiurus nebulosus) records in Aotearoa. (b) Distribution of kōura (Paranephrops planifrons) records in Aotearoa. Data from the New Zealand Freshwater Fish Database (Stoffels, 2022).

To identify conditions under which native species are more likely to persist given concurrent climatic change and interactions with introduced species, refugia habitat needs to be identified, ideally within the current range of native species but also further beyond if habitat within the current range becomes unsuitable (Gallardo et al., 2017). Attempts to identify climate change refugia provide an initial filter of environmental suitability (Ashcroft, 2010; Barrows et al., 2020; Keppel et al., 2012). Biotic interactions provide a second level of filtering; for example, the presence of predators or competitors (native and invasive) can reduce suitability, while the presence of prey or host species may be prerequisites for potential refugia (Van Der Putten et al., 2010; Wisz et al., 2013). Identifying all biotic interactions and how they may respond to climate change is effectively impossible, but the presence of invasive species (as a novel competitor, predator, or pathogen) known to threaten native species is likely to be a strong limiting factor to the quality of refugia habitat (McCarthy et al., 2021). Invasive predators in particular can impact lower trophic levels by intensified predator-controlled trophic cascades (Mpanza et al., 2024), which can be the result if refugia of prey are lost, in conjunction with the habitat of predators being sustained or increased. In addition, the presence of invasive predators can have non-consumptive effects on resident prey species, by negatively impacting on feeding, habitat use, and reproduction (MacNeil and Briffa, 2019).

Ecological niche models (ENMs) are widely used to predict the current and potential distribution of species and quantify habitat suitability. Ecological niche models typically use correlations between environmental predictors and species presences and absences (or pseudo-absences or background points) to understand where conditions are most suitable for a species. There are many algorithms available, such as Random Forests, GLM, and MaxEnt, that all make slightly different assumptions but produce relative probability of occurrence or habitat suitability scores (Elith and Leathwick, 2009). Extrapolating ENMs in time or space can be problematic because the assumption of environmental equilibrium may be violated, and novel species interactions and adaptation are not accounted for (Austin, 2007). Nevertheless, extrapolation is essential for conservation planning, and tools exist to quantify issues such as the extent to which extrapolation beyond training data has occurred (Velazco et al., 2024).

Here we address the need to identify refugia from both climate change and invasive species using a novel refugia habitat metric (RHM). We use EMNs to describe current and future habitat suitability for a native and an invasive species. We then use the habitat suitability of the invasive species to down-weight native habitat suitability. Our refugia habitat metric can be used to identify habitats that are both suitable for native species and unsuitable for invasive species. The metric can be applied using current or future environmental and habitat data and can include one or many invasive species. The metric is also agnostic to the modelling approach used to generate habitat suitability scores. We demonstrate the utility of the metric via a case study of an endemic freshwater crayfish, northern kōura (Paranephrops planifrons, hereafter kōura) in Aotearoa / New Zealand (hereafter Aotearoa).

Kōura may be exposed to warming-related stressors and predation by the invasive brown bullhead catfish (Ameiurus nebulosus, hereafter bullhead catfish). Kōura are considered a taonga (treasured) species by Māori (the indigenous people of Aotearoa), and they support important mahinga kai practices (customary food gathering; Kusabs and Quinn, 2009). Kōura occur across much of the North Island of Aotearoa and west of the Southern Alps in the South Island (Fig. 1); they prefer cool, slow-moving water, woody debris cover, and coarse substrates that provide a variety of food items (e.g. invertebrates, plant roots, and leaves) and shelter (Jowett et al., 2008; Kusabs et al., 2015; Parkyn et al., 2002; Usio and Townsend, 2000). Water temperature is an important factor dictating metabolic activity, growth, and feeding behaviour for crayfish (Verhoef and Austin, 1999). In the North Island, kōura have been observed where water temperature ranges from 6–18 °C in streams in native forest cover and from 5–25 °C in pasture streams (Parkyn et al., 2002). Jones (1981) reported the optimum mean daily temperature for kōura to be 19 °C, with a mean critical upper limit of 31.9 °C (Simons, 1984). Freshwater crayfish are vulnerable to climate change via increasing water temperatures and their limited adaptive capacity (Hossain et al., 2018). Kōura are currently declining (their official conservation status is “At Risk: Declining”).

Although population estimates for both kōura and bullhead catfish are uncertain, bullhead catfish are implicated in the decline in kōura in the locations where they have invaded and become established (Clearwater et al., 2014; Dedual, 2019). Although there are few studies in Aotearoa detailing stomach contents of bullhead catfish where kōura are present, Barnes and Hicks (2003) reported that 64 % of large catfish (<250 mm fork length) from rocky shoreline areas of Lake Taupō contained kōura, indicating the catfish is a major predator. Bullhead catfish are native to eastern North America, where they occur in a variety of habitats but particularly in still or slow-moving waterbodies with muddy substrates (Dedual, 2019). The range expansion of bullhead catfish is facilitated by its ability to survive in water temperatures up to 37.5 °C, spawn over a wide temperature range (14–29 °C; Blumer, 1985), and tolerate low oxygen concentrations and poor water quality for prolonged periods and by its broad omnivorous diet (Scott and Crossman, 1973). The bullhead catfish has been present in Aotearoa since 1877, when it was deliberately introduced for sport-fishing (McDowall, 1994). The current distribution of bullhead catfish in Aotearoa is largely confined to the upper North Island (Fig. 1), with only two restricted observations in the South Island (NIWA, 2020); however, in recent decades, the bullhead catfish has expanded its range (Hicks and Allan, 2018). Until 1985, the bullhead catfish was restricted to the lower parts of the Waikato River in the North Island and to Lake Mahinapua in the South Island (Fig. 1). In 1985, it was detected in Lake Taupō (Barnes and Hicks, 2003). It has since spread throughout the Waikato River and has reached lakes nearly 100 km to the northeast in the Rotorua region (lakes Rotorua and Rotoiti in 2016 and 2018, respectively) and multiple catchments in the north of the North Island (Hicks and Allan, 2018).

In Aotearoa, many freshwater native fish could face extinction or near-extinction from loss of viable habitat due to climate change (Canning et al., 2025), which will be exacerbated by interactions with other stressors, such as declining water quality and non-climate-related habitat loss (Ling, 2010). Given the broader temperature tolerances of bullhead catfish than kōura, and projected increases in temperature for Aotearoa under climate change, kōura face compounding pressures from climate change itself (NIWA 2020) and from a highly invasive predator favoured by such change. Climate change is likely to result in a net decrease in suitable habitat for kōura and a net increase in suitable habitat for bullhead catfish. Here, we describe the RHM, illustrate how habitat suitability for kōura and bullhead catfish is likely to change by the end of the century, and apply the RHM to identify a potential refugia habitat for kōura.

2 Materials and methods

Firstly, we describe the RHM mathematically, along with its key properties. We then use ecological niche models to identify suitable habitat across Aotearoa for bullhead catfish and kōura for the present day and under predicted climate conditions for 2080–2100. We then calculate RHM scores for Aotearoa to identify potential current and future refugia for kōura.

2.1 Refugia habitat metric

Refugia habitat metric (RHM) scores were calculated using habitat suitability maps produced via ecological niche modelling (explained below), where each raster pixel has a habitat suitability score ranging from 0–1 (0 unsuitable – 1 suitable). Habitat suitability maps are required for the native species and the n invasive species of interest. The RHM can only be calculated for pixels that have a score for all species of interest. The fundamental concept behind the RHM is that habitat suitability for the native species should be down-weighted based on the habitat suitability score(s) of invasive species. We chose to use a negative exponential weighting, as we assumed the presence of an increasing number of invasive species would have a diminishing effect on the native species. Firstly, a weighting parameter (ω) was calculated based on the habitat suitability score(s) of invasive species:

$$\begin{array}{}\text{(1)}& \mathit{\omega }={\sum }_{i=\mathrm{1}}^n{\text{InvHSS}}_i\cdot {\mathit{\alpha }}_i,\end{array}$$

where InvHSS_i is the habitat suitability score of the ith invasive species in each pixel, n is the total number of invasive species, and α_i describes the strength of interactions between the ith invasive species and the native species. Higher values of α_i indicate that high habitat suitability scores for the invasive species should more strongly negatively affect the habitat suitability scores of the native species. An α_i of 0.0 results in an invasive species having no effect on the RHM (RHM = NatHSS). An α_i of 1.0 can be used by default. The RHM was then calculated by reweighting the native species' habitat suitability score:

$$\begin{array}{}\text{(2)}& \text{RHM}=\text{NatHSS}\cdot {\mathit{\beta }}^{-\mathit{\omega }},\end{array}$$

Figure 2 Example of how the RHM responds to the invasive species habitat suitability score (InvHSS) for a range of invasive species habitat suitability score weighting parameter (α) values (a: 0; b: 0.5; c: 1; d: 2) when a single invasive species is present. β fixed at 2.0.

Figure 3 Example of how the RHM responds to the invasive species habitat suitability score (InvHSS) for a range of invasive species habitat suitability score weighting parameter (α) values and shape parameter (β) values (a: 1; b: 2; c: 4; d: 8) when a single invasive species is present. All values of NatHSS were fixed at 1.0. See Fig. A1 for the equivalent figure but with multiple invasive species.

where NatHSS is the habitat suitability score of the native species in a given pixel and β is a global weighting parameter. Values of β must be 1.0 (all invasive species have no effect on native species) or greater. The RHM metric has the following properties:

• It ranges between 0.0 and 1.0.

• If the invasive species has a habitat suitability score of 0.0, the NatHSS is equal to the RHM.

• If α is 0.0, the NatHSS is equal to the RHM.

• If β is 1.0, the NatHSS is equal to the RHM.

• If α = 1, β = 2, and ω = 1, the RHM = NatHSS/2. We recommend α = 1 and β = 2 as the default values.

• There is no limit to the number of invasive species that can be included.

• Down-weighting follows a negative exponential relationship; that is, as ω increases, the relative effect on the RHM deceases. (Figs. 3 and A1).

2.2 Case study: endemic kōura and invasive bullhead catfish

2.2.1 Species records

Occurrence records for bullhead catfish and kōura were downloaded from the New Zealand Freshwater Fish Database (NZFFD; Stoffels, 2022). Occurrence data were filtered to include records from 1986 onwards, to align with environmental data (see Sect. 2.2.2 for description of environmental data). Data were cleaned by removing duplicates and the following records: those within 1 km of country centroids, those within 1 km of the Global Biodiversity Information Facility (GBIF) headquarters in Copenhagen, those within 100 m of known biodiversity institutions (these coordinates emerge via automation when just the country is listed in the primary source; Zizka et al., 2019), and those located in the sea or at 0,0 longitude/latitude (Zizka et al., 2019). We only included “human observation” records (i.e. we discarded museum specimens, individuals at zoos, and other non-observation records). After these filters were applied, we were left with 7344 occurrence points (native range = 5766; introduced = 1578; Fig. 1). The NZFFD does not distinguish between Northern (Paranephrops planifrons) and Southern kōura (P. zealandicus; the only other crayfish species in Aotearoa); however, the records can be distinguished due to the two species maintaining distinct geographic ranges (Hopkins, 1970). After filtering the occurrence records, we were left with 255 and 3654 presences for bullhead and kōura, respectively (Fig. 1).

2.2.2 Environmental variables

We focused on environmental and habitat variables that are ecologically important to freshwater species (Table A1). We included river-specific variables from the Freshwater Environments of New Zealand (FENZ) database. FENZ consists of a network representing all river reaches in Aotearoa, and it contains over 500 000 individual reaches (Leathwick et al., 2010). Each reach is described in terms of a series of modelled physiochemical variables such as upstream forest cover and reach slope. From the FENZ database we included the following variables: catchment area, Strahler order, median, downstream distance to the coast, upstream calcium concentration, proportion upstream catchment covered by lakes, upstream indigenous forest cover, and reach habitat type (1 = still; 2 = backwater; 3 = pool; 4 = run; 5 = riffle; 6 = rapid; 7 = cascade). Additionally, we included dissolved reactive phosphorus (DRP; Whitehead et al., 2022) and total suspended solids (TSS; Unwin and Larned, 2013) from the River Maps database (Whitehead and Booker, 2019). There are no future predictions available for the FENZ variables; therefore we assumed they would remain relatively static through time. We included ecologically important climate variables (mean daily temperature, total annual rainfall, mean solar radiation, and potential evapotranspiration deficit) from the New Zealand Climate Projection Dataset for the period 1986–2005 (Gibson et al., 2024). Air temperature variables were used as surrogates for instream temperature, as instream temperature data were not available and air temperature is correlated with instream temperature (McGarvey et al., 2018). The data have been downsampled at 5 km resolution specifically for Aotearoa from global models. Since the FENZ data were in vector format and the climate data were gridded raster data, we extracted climate data along each FENZ reach; if a reach intersected more than one grid pixel, we took the mean value across intersecting pixels. All absolute pairwise correlations among environmental variables were less than 0.7. Predictors were not transformed prior to analysis.

To understand how habitat suitability is likely to change later this century, we used projections for mean daily temperature, total annual rainfall, mean solar radiation, and potential evapotranspiration deficit for the period 2080–2100 from the New Zealand Climate Projection Dataset, which consists of climate projections created using six global models downscaled to provide New Zealand-specific projections at a 5 km grid. Given the uncertainty in future climate projections, we used mean estimates derived from the six models (NorESM2-MM, AWI-CM-1-1-MR, GFDL-ESM4, ACCESS-CM2, EC-Earth3, CNRM-CM6-1) for the Shared Socioeconomic Pathway (SSP) SSP3-7.0 scenario (Gibson et al., 2024). We have assumed the non-climate predictors will not change into the future, which is likely a best-case scenario for kōura; the National Policy Statement for Freshwater Management 2020, which is the current overriding freshwater environmental protection legislation in Aotearoa / New Zealand, sets out the requirements that ecosystem health of water bodies must be either maintained or improved, which provides some justification for this assumption.

2.2.3 Ecological niche modelling

To predict habitat suitability across Aotearoa for bullhead catfish and kōura, we implemented a maximum entropy (MaxEnt) model (Phillips et al., 2006). MaxEnt is a presence–background approach based on an inhomogeneous Poisson process specifically designed for presence-only data (Phillips et al., 2017). We used MaxEnt due to its high predictive performance when using presence-only data and when looking to extrapolate in space or time (Ahmadi et al., 2023; Heikkinen et al., 2012; Valavi et al., 2022; Velazco et al., 2024).

When creating a calibration area for presence-only models, it is important to delineate areas where species have had the opportunity to colonize (Barve et al., 2011); therefore we constrained the calibration area to catchments that contained at least one occurrence record. We sampled 10 000 background reaches at random from within calibration areas (Hill et al., 2017; Merow et al., 2013).

We systematically varied the hyper-parameter settings to tune the MaxEnt model and tested regularization multiplier settings of 0.5, 1, 2, and 3 and feature class settings of linear (“l”), hinge (“h”), quadratic (“q”), product (“p”), and threshold (“t”). The hyper-parameter combination that maximized the true skill statistic (TSS) was selected, and the final model had the regularization multiplier set to 2.0 and 0.5 for bullhead catfish and kōura, respectively, and feature class combination “lhp” for both species.

Model discrimination was assessed using K-fold cross-validation (N=5). Firstly, we used the area under the receiver operator curve (AUC), where values range from 0 to 1, with a value of 0.5 indicating no better than random predictive ability and a value of 1 indicating perfect prediction (Fielding and Bell, 1997). Secondly, we used the Continuous Boyce Index (CBI), which is favoured in studies that use presence-only data and where values range from −1 to 1, with higher values indicating better performance. Thirdly, we calculated sensitivity, specificity, and the true skill statistic (TSS), where values closer to 1.0 indicate better model discrimination. Because we extrapolated habitat suitability to areas and times outside the training data, we quantified the degree of extrapolation in environmental space between the training and projection data, using the shape metric (Velazco et al., 2024). Shape is calculated using the multivariate distance from each extrapolated point to the nearest training data point, in environmental space. If a given extrapolated point has a shape value of 100, it is 100 times further away from the nearest training data point than the averaged Mahalanobis distance between training data points and the centroid of the training data (Velazco et al., 2024). Non-analogue environments do not necessarily equate to invalid predictions, since invasive species have often shown their ability to colonize new environments, but identify areas where predictions may be uncertain.

To understand how the amount of habitat available to kōura and bullhead catfish is likely to change between now and 2080–2100 and to avoid having to pick a single binarizing threshold, we calculated total available habitat for the two time periods across the full range of binarizing thresholds (0–1). Finally, we calculated the RHM for current and future time periods. All analysis was carried out in R v.4.3.2 (R Core Team, 2020). We used the tidyverse v2.0.0 (Wickham et al., 2019), terra v1.8.54 (Hijmans, 2024), and sf v1.0.21 (Pebesma, 2018) packages for data manipulation; the flexsdm v1.3.6 (Velazco et al., 2022) package for model fitting; and the tidyterra v0.7.2 (Hernangómez, 2023) and patchwork v1.3.1 (Pedersen, 2021) packages for visualization.

Figure 4 Predicted habitat suitability of Aotearoa for (a–c) brown bullhead catfish (Ameiurus nebulosus) and (d–f) kōura (Paranephrops planifrons) under current climatic conditions (a, d) and in the period 2080–2100 (b, e), with the change in suitability between current and future conditions (c, f).

Figure 5 Proportion of Aotearoa that is suitable habitat for kōura (Paranephrops planifrons) and brown bullhead catfish (Ameiurus nebulosus) using a range of threshold values. Current range refers to habitat within the current range of kōura.

3 Results

3.1 Current and future habitat suitability

The accuracy of the MaxEnt models varied for bullhead catfish and kōura (AUC = 0.92, CBI = 0.92, sensitivity = 0.92, specificity = 0.78, true skill statistic = 0.70 for bullhead catfish; AUC = 0.81, CBI = 0.99, sensitivity = 0.73, specificity = 0.73, true skill statistic = 0.46 for kōura). Currently, from known species records, bullhead catfish are largely confined to the upper North Island, and our habitat suitability modelling confirmed the northern and to a lesser extent the western North Island are highly suitable for bullhead catfish, while in the South Island there are small areas of suitability, particularly along the east and south of the island (Fig. 4). For kōura, known records show this currently occurring across much of the North Island and west of the Southern Alps in the South Island, and our habitat suitability modelling confirmed much of coastal Aotearoa had good to high (greater than 0.5) habitat suitability. For both bullhead catfish and kōura, the mountainous central South Island and eastern North Island had low suitability. Under future climate scenarios, bullhead catfish habitat suitability either remained similar to current conditions or increased in suitability for the period 2080–2100 (Fig. 4). Furthermore, there were significant increases in habitat suitability, particularly in the western North Island. In contrast, for kōura, habitat suitability declined across most of the current range in the North Island, with the South Island range maintaining high suitability. Outside of the current range, the South Island largely increased in suitability, as did the central and eastern North Island (Fig. 4). Overall, under climate change, the amount of Aotearoa favourable for kōura is projected to stay relatively stable but with shifts in where is suitable, while for catfish there is an overall increase in habitat suitability (Figs. 4 and 5).

Environmental conditions in the South Island were less similar to model training data than conditions in the North Island, as indicated by higher shape values (measure of environmental dissimilarity between training and extrapolated data) for bullhead catfish, indicating a greater level of uncertainty in model projections in the South Island. In contrast, for kōura, most areas had low (< 50) shape values, which is unsurprising given that the projection and training areas overlapped significantly (Fig. A3).

Figure 6 Partial dependence plots for the MaxEnt model predicting habitat suitability of Aotearoa for brown bullhead catfish (Ameiurus nebulosus) and kōura (Paranephrops planifrons).

For bullhead catfish, the most important predictors of habitat suitability were distance to coast, Strahler order, annual rainfall, and reach habitat (Fig. A4); for kōura, they were distance to coast, mean annual temperature, catchment area, and Strahler order (Fig. A5). On average, higher habitat suitability scores for bullhead catfish were associated with low distance to coast, high Strahler stream order (larger rivers), low annual rainfall, low habitat scores (still water/backwater habitats), and higher mean annual temperature (Fig. 6). Kōura habitat suitability was positively associated with low distance to coast, intermediate temperatures, larger catchments, and low Strahler order (smaller streams; Fig. 6).

Figure 7 Potential refugia (0 low priority, 1 high priority) for kōura (Paranephrops planifrons) based on current and future habitat suitability for kōura and brown bullhead catfish (Ameiurus nebulosus), for current (a) and future (b) conditions. Black polygons indicate kōura occurrences with a 20 km buffer, approximating the current distribution.

3.2 Refugia habitat metric

The refugia habitat metric indicates that under current conditions there are significant areas of high-suitability refugia in the current range of kōura, particularly in the western and southern North Island and across the north and west of the South Island (Fig. 7a). Additionally, there is significant refugia habitat outside the species current range, across the eastern South Island. Under future conditions, potential refugia habitat contracts, particularly within the current range in the northern North Island, although significant refugia still remain in the western North Island and South Island (Fig. 7b). There is also significant refugia habitat present outside the current range in the eastern South Island.

4 Discussion

The RHM introduced in this study represents an advancement in ecological niche modelling by integrating the influence of invasive species into projections of native species' habitat suitability under climate change. While most ENMs focus on abiotic predictors of habitat suitability, such as temperature and precipitation (Elith and Leathwick, 2009), they often overlook the role of biotic interactions (but see Pollock et al., 2014, for a description for joint species distribution models). By down-weighting the native species' habitat suitability in areas with high invasive species habitat suitability, the RHM accounts for a key ecological constraint on native persistence. Our approach explicitly acknowledges that habitat suitability is not solely a function of environmental drivers but also biotic interactions (Araújo and Luoto, 2007). The RHM offers a flexible and scalable approach that directly modifies the native species' suitability surface based on invasion potential, producing maps that simultaneously reflect abiotic and biotic constraints. The RHM is not intended to replace existing tools but rather to complement and enhance them. If used alongside mechanistic models, physiological threshold data, or population viability analyses, the metric could help provide a more complete picture of species persistence potential. Furthermore, integrating it with landscape connectivity or habitat fragmentation models could identify corridors that facilitate native species' movement to less invaded and climatically stable areas. Taking a multi-pronged approach would allow a more comprehensive conservation strategy that considers both environmental and ecological dimensions of habitat suitability.

Although we used brown bullhead catfish and northern kōura in our case study, the RHM is designed to be broadly applicable across taxonomic groups, ecosystem types, and geographic scales. It can be used in the context of different invasive pressures, from aggressive plant species to aquatic invaders disrupting native fish assemblages. As long as ENMs can be developed for both native and invasive species, the metric can identify refugia in terrestrial, freshwater, or marine systems. Moreover, the approach is scalable: it can be applied locally for site-level management or extended to regional or global biodiversity assessments (Thuiller et al., 2005; Tingley et al., 2014). This flexibility makes the RHM applicable to a wide range of ecological and conservation contexts, from targeted reserve design to continental-scale invasion risk mapping.

From a practical standpoint, the metric provides a valuable tool for conservation managers and policymakers tasked with navigating the complex realities of global change (Rillig, 2025). By identifying areas that are likely to offer both climatic stability and low invasion pressure, the RHM can help focus efforts where native species are most likely to persist. This has direct applications in habitat restoration, reserve design, and invasive species management and could guide assisted migration efforts, especially where funding or capacity is limited (Heller and Zavaleta, 2009). The metric could also inform spatial prioritization in national biodiversity strategies or contribute to global assessments of climate refugia (Morelli et al., 2016). Moreover, it aligns with emerging conservation frameworks that advocate for proactive, integrated approaches that address multiple drivers of biodiversity loss simultaneously (IPBES, 2019). With the recent (2023) invasion of the freshwaters of Aotearoa / New Zealand by the gold/Asian clam (Corbicula fluminea), an ecosystem engineer with the potential to outcompete and replace culturally valued native bivalves (Somerville et al., 2025), conservation translocations involving moving native species out of the pathway of the ongoing range expansion of such potentially damaging invaders is a current subject of debate within government proception agencies (MPI, 2023). The RHM could provide a highly useful and practical conservation tool in these endeavours.

While the RHM is a useful tool for identifying areas where native species may face less intense invasion pressure, it is important to consider the limitations of ecological assumptions underpinning its design. Firstly, the metric assumes that areas of high suitability for invasive species are associated with greater ecological risk to natives, primarily through competitive displacement or predation. However, this relationship is likely context-dependent, varying with species traits, resource overlap, and environmental conditions (Lee-Yaw et al., 2022; MacDougall et al., 2009). For example, the effects of catfish on kōura are likely to be greatest where eels (Anguilla dieffenbachia and A. australis) are not present or are rare (e.g. Rotorua and Taupō lakes); where eels are abundant, they generally dominate catfish, reducing their effect on kōura (Kusabs pers. obs.). Some introduced species may have minimal impact on certain native taxa despite co-occurrence, while others may be highly disruptive even at low densities (Ricciardi et al., 2013), highlighting the need for input from experts when deciding which invasive species to include in the modelling or how to weight different species. Secondly, the relative effect of invasive species may be context-dependent. For instance, in one of the few Aotearoa / New Zealand studies detailing stomach contents of bullhead catfish in invaded areas harbouring significant kōura populations, the stomach contents of 6247 catfish from Lake Taupō revealed that 64 % of large catfish (<250 mm fork length) from rocky shoreline areas contained kōura, while only 15 % of the same-size catfish from weedy areas of shoreline contained kōura (Barnes and Hicks, 2003). This indicates the importance of the habitat template in which the interaction is taking place in determining if the invader is “weak” or “strong” in terms of predatory impact (Kumschick et al., 2011). A strength of the RHM is the ability to incorporate variability in the effect of different invasive species via the α parameter (although determining α values may be challenging if one wants to move beyond qualitative rankings of invasive species impacts). Thirdly, the reliability of the RHM is contingent on the quality of the underlying ENMs. The ENM approach is biologically naive and ignores dispersal limitation and multi-species interactions (Wisz et al., 2013). For example, the area invaded by catfish may be less than that predicted by the model due to dispersal limitation and human interventions (e.g. eradication efforts). Ultimately, if all these mechanisms are to be included, more process-based models are required (Pilowsky et al., 2022). Furthermore, model outputs are influenced by the resolution and quality of environmental predictors, the choice of algorithms, and the completeness of occurrence data (Elith and Graham, 2009). For example, the river reaches in our case study are approximately 700 m in length, which will not capture microhabitat refugia, such as small headwater and spring-fed streams, and which will allow kōura to persist within their current range (Clearwater et al., 2014). Fourthly, the RHM is static and does not account for temporal dynamics such as lagged responses to climate change or adaptive evolution in native species. These limitations highlight the importance of using the metric as a component of an integrative framework, rather than a stand-alone predictive tool. Finally, uncertainties inherent in climate projections, land-use scenarios, or invasive species distributions can propagate non-additively (Buisson et al., 2010). Sensitivity analyses and ensemble modelling could help quantify and mitigate some of these uncertainties, but users of the RHM should remain cautious and transparent about potential sources of error.

Future work could integrate dynamic modelling frameworks (e.g. process-based or agent-based models) or couple the metric with trait-based or mechanistic approaches that better capture species-specific vulnerabilities and adaptive capacities (Kearney and Porter, 2009; Pilowsky et al., 2022). Incorporating temporal dynamics, such as simulating invasive species' spread rates or native species' dispersal abilities, would allow more realistic projections under changing environmental conditions. Validation of the RHM represents a next step for advancing its application and credibility. Empirical testing could involve comparing adjusted suitability predictions to known patterns of native persistence or decline in the presence of invasive species across environmental gradients. Long-term monitoring data, experimental removals, and invasion chronosequences offer valuable opportunities to test whether areas identified as dual refugia correspond with higher native survival or reproductive success. These directions represent promising avenues for future research and development.

Finally, in the bullhead catfish – kōura case study, we showed that significant parts of the North Island and parts of the eastern South Island of Aotearoa are predicted to become suitable for bullhead catfish under climate change. The predictions across the North Island are similar to those predicted in Canning et al. (2025), although they did not predict any suitable habitat in the South Island. For kōura, climate change is predicted to have mixed effects; some parts of Aotearoa will become more suitable, while others will become less suitable. Areas of increasing habitat suitability are concentrated in the cooler southern parts of the country, whereas habitat suitability decreased in northern lowland coastal areas, which, while currently suitable, may become too warm in the future. Bullhead catfish pressure is likely to be greatest (i) in shallow lakes/rivers where the bed substrate is composed mainly of sand/mud and (ii) in lakes that suffer from hypolimnetic deoxygenation; in these waterbodies, kōura have very limited, if any, rocky habitat or depth refugia from bullhead catfish (Francis, 2019). As discussed earlier, there are likely to be refugia within the current range of kōura, in the form of small headwater and spring-fed streams which were not necessarily identified here given the fine-resolution data that are needed to tease apart microhabitat variability. Many of these systems have large kōura populations and could easily be engineered to exclude predators, for example, using weirs (Ian A. K. Kusabs, personal observation). It is unlikely that kōura will be able to track refugia outside their current range given dispersal constraints. For example, the area extending northeast from Rotorua–Taupō (eastern North Island) has large areas of potentially suitable habitat, yet kōura are absent. McDowall (2005) hypothesized that this gap in the distribution is due to dispersal limitation and a failure to recolonize the area after volcanic eruptions locally extirpated kōura populations (McDowall, 2005). Translocations are increasing being considered for freshwater conservation in Aotearoa / New Zealand (Rayne et al., 2025) and provide one approach to facilitate climate tracking. For kōura there is an existing long history of translocating for fisheries purposes (Rayne et al., 2020).

We have developed and demonstrated the use of a novel metric to identify habitat refugia from the combined threats of climate change and invasive species. The RHM has the potential to be used in a wide variety of conservation settings and should be beneficial to species managers and restoration practitioners.

Appendix A

Figure A1 Example of how the RHM responds to multiple invasive species (ω) for a range of invasive species habitat suitability score weighting parameter (α) values and shape parameter (β) values (a: 1; b: 2; c: 4; d: 8). All values of NatHSS were fixed at 1.0.

Table A1 Description of the environmental and climate variables used in the MaxEnt models of brown bullhead catfish (Ameiurus nebulosus) and kōura (Paranephrops planifrons) habitat suitability.

Figure A2 Pairwise correlations among environmental predictors used in MaxEnt models for brown bullhead catfish (Ameiurus nebulosus) and kōura (Paranephrops planifrons) habitat suitability. See Table A1 for a description of variables.

Figure A3 Shape metric scores across Aotearoa for kōura (Paranephrops planifrons) and brown bullhead catfish (Ameiurus nebulosus). “2080–2100” refers to the SSP3-7.0 Shared Socioeconomic Pathway projection for the period 2080–2100.

Figure A4 Variable importance scores for the brown bullhead catfish (Ameiurus nebulosus) MaxEnt model.

Figure A5 Variable importance scores for the kōura (Paranephrops planifrons) MaxEnt model.